Skip to content

Commit

Permalink
Merge pull request #1029 from CEED/jrwrigh/stg-primitive
Browse files Browse the repository at this point in the history 
fluids: Make Strong STG compatible with primitive variables.
  • Loading branch information
@jrwrigh
jrwrigh authored Jul 20, 2022
2 parents dc805cc + efe9d85 commit 39d18f8
Show file tree
Hide file tree
Showing 11 changed files with 143 additions and 66 deletions.
5 changes: 5 additions & 0 deletions examples/fluids/README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -534,6 +534,11 @@ For the Density Current, Channel, and Blasius problems, the following common com
- Developer option to test properties
- `false`
- boolean

* - `-primitive`
- Use primitive variables (pressure, velocity, temperature) instead of conservative variables (density, momentum, total energy)
- `false`
- boolean
:::


Expand Down
2 changes: 1 addition & 1 deletion examples/fluids/problems/blasius.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,7 +126,7 @@ static PetscErrorCode ModifyMesh(MPI_Comm comm, DM dm, PetscInt dim,
faces[1]+1, *num_node_locs);
if (*num_node_locs > faces[1] +1) {
ierr = PetscPrintf(comm, "WARNING: y_node_locs_path has more locations (%d) "
"than the mesh has nodes (%d). This maybe unintended.",
"than the mesh has nodes (%d). This maybe unintended.\n",
*num_node_locs, faces[1]+1); CHKERRQ(ierr);
}
PetscScalar max_y = (*node_locs)[faces[1]];
Expand Down
10 changes: 9 additions & 1 deletion examples/fluids/problems/newtonian.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -145,6 +145,14 @@ PetscErrorCode NS_NEWTONIAN_IG(ProblemData *problem, DM dm, void *ctx) {
problem->apply_vol_ifunction.qfunction_loc = IFunction_Newtonian_Prim_loc;
problem->apply_vol_ijacobian.qfunction = IJacobian_Newtonian_Prim;
problem->apply_vol_ijacobian.qfunction_loc = IJacobian_Newtonian_Prim_loc;
problem->apply_inflow.qfunction = BoundaryIntegral;
problem->apply_inflow.qfunction_loc = BoundaryIntegral_loc;
problem->apply_inflow_jacobian.qfunction = BoundaryIntegral_Jacobian;
problem->apply_inflow_jacobian.qfunction_loc = BoundaryIntegral_Jacobian_loc;
problem->apply_outflow.qfunction = PressureOutflow;
problem->apply_outflow.qfunction_loc = PressureOutflow_loc;
problem->apply_outflow_jacobian.qfunction = PressureOutflow_Jacobian;
problem->apply_outflow_jacobian.qfunction_loc = PressureOutflow_Jacobian_loc;
} else {
problem->ics.qfunction = ICsNewtonianIG;
problem->ics.qfunction_loc = ICsNewtonianIG_loc;
Expand Down Expand Up @@ -288,7 +296,7 @@ PetscErrorCode NS_NEWTONIAN_IG(ProblemData *problem, DM dm, void *ctx) {
newtonian_ig_ctx->Ctau_E = Ctau_E;
newtonian_ig_ctx->stabilization = stab;
newtonian_ig_ctx->is_implicit = implicit;
newtonian_ig_ctx->primitive = prim_var;
newtonian_ig_ctx->is_primitive = prim_var;
ierr = PetscArraycpy(newtonian_ig_ctx->g, g, 3); CHKERRQ(ierr);

CeedQFunctionContextCreate(user->ceed, &problem->ics.qfunction_context);
Expand Down
14 changes: 11 additions & 3 deletions examples/fluids/problems/stg_shur14.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -470,13 +470,21 @@ PetscErrorCode StrongSTGbcFunc(PetscInt dim, PetscReal time,
PetscFunctionReturn(0);
}

PetscErrorCode SetupStrongSTG(DM dm, SimpleBC bc, ProblemData *problem) {
PetscErrorCode SetupStrongSTG(DM dm, SimpleBC bc, ProblemData *problem,
Physics phys) {
PetscErrorCode ierr;
DMLabel label;
const PetscInt comps[] = {0, 1, 2, 3};
const PetscInt num_comps = 4;
PetscFunctionBeginUser;

PetscInt comps[5], num_comps=4;
if (phys->primitive) {
// {1,2,3,4} for u, v, w, T
for(int i=0; i<4; i++) comps[i] = i+1;
} else {
// {0,1,2,3} for rho, rho*u, rho*v, rho*w
for(int i=0; i<4; i++) comps[i] = i;
}

ierr = DMGetLabel(dm, "Face Sets", &label); CHKERRQ(ierr);
// Set wall BCs
if (bc->num_inflow > 0) {
Expand Down
3 changes: 2 additions & 1 deletion examples/fluids/problems/stg_shur14.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,8 @@ extern PetscErrorCode SetupSTG(const MPI_Comm comm, const DM dm,
const CeedScalar P0,
const CeedScalar ynodes[], const CeedInt nynodes);

extern PetscErrorCode SetupStrongSTG(DM dm, SimpleBC bc, ProblemData *problem);
extern PetscErrorCode SetupStrongSTG(DM dm, SimpleBC bc, ProblemData *problem,
Physics phys);

extern PetscErrorCode SetupStrongSTG_QF(Ceed ceed, ProblemData *problem,
CeedInt num_comp_x, CeedInt num_comp_q, CeedInt stg_data_size,
Expand Down
2 changes: 1 addition & 1 deletion examples/fluids/qfunctions/channel.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,7 @@ CEED_QFUNCTION(ICsChannel)(void *ctx, CeedInt Q,
for (CeedInt i=0; i<Q; i++) {
const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]};
State s = Exact_Channel(3, 0., x, 5, ctx);
if (context->newtonian_ctx.primitive) {
if (context->newtonian_ctx.is_primitive) {
q0[0][i] = s.Y.pressure;
for (CeedInt j=0; j<3; j++)
q0[j+1][i] = s.Y.velocity[j];
Expand Down
2 changes: 1 addition & 1 deletion examples/fluids/qfunctions/densitycurrent.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,7 +152,7 @@ CEED_QFUNCTION(ICsDC)(void *ctx, CeedInt Q,
for (CeedInt i=0; i<Q; i++) {
const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]};
State s = Exact_DC(3, 0., x, 5, ctx);
if (context->newtonian_ctx.primitive) {
if (context->newtonian_ctx.is_primitive) {
q0[0][i] = s.Y.pressure;
for (CeedInt j=0; j<3; j++)
q0[j+1][i] = s.Y.velocity[j];
Expand Down
130 changes: 84 additions & 46 deletions examples/fluids/qfunctions/newtonian.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -790,12 +790,20 @@ CEED_QFUNCTION(BoundaryIntegral)(void *ctx, CeedInt Q,

const NewtonianIdealGasContext context = (NewtonianIdealGasContext) ctx;
const bool is_implicit = context->is_implicit;
State (*StateFromQi)(NewtonianIdealGasContext gas,
const CeedScalar qi[5], const CeedScalar x[3]);
State (*StateFromQi_fwd)(NewtonianIdealGasContext gas,
State s, const CeedScalar dqi[5],
const CeedScalar x[3], const CeedScalar dx[3]);
StateFromQi = context->is_primitive ? &StateFromY : &StateFromU;
StateFromQi_fwd = context->is_primitive ? &StateFromY_fwd : &StateFromU_fwd;


CeedPragmaSIMD
for(CeedInt i=0; i<Q; i++) {
const CeedScalar U[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]};
const State s = StateFromU(context, U, x_i);
const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
State s = StateFromQi(context, qi, x_i);

const CeedScalar wdetJb = (is_implicit ? -1. : 1.) * q_data_sur[0][i];
// ---- Normal vect
Expand All @@ -811,12 +819,12 @@ CEED_QFUNCTION(BoundaryIntegral)(void *ctx, CeedInt Q,

State grad_s[3];
for (CeedInt j=0; j<3; j++) {
CeedScalar dx_i[3] = {0}, dU[5];
CeedScalar dx_i[3] = {0}, dqi[5];
for (CeedInt k=0; k<5; k++)
dU[k] = Grad_q[0][k][i] * dXdx[0][j] +
Grad_q[1][k][i] * dXdx[1][j];
dqi[k] = Grad_q[0][k][i] * dXdx[0][j] +
Grad_q[1][k][i] * dXdx[1][j];
dx_i[j] = 1.;
grad_s[j] = StateFromU_fwd(context, s, dU, x_i, dx_i);
grad_s[j] = StateFromQi_fwd(context, s, dqi, x_i, dx_i);
}

CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3];
Expand Down Expand Up @@ -846,11 +854,15 @@ CEED_QFUNCTION(BoundaryIntegral)(void *ctx, CeedInt Q,
// -- Total Energy Density
v[4][i] = -wdetJb * Flux[4];

jac_data_sur[0][i] = s.U.density;
jac_data_sur[1][i] = s.Y.velocity[0];
jac_data_sur[2][i] = s.Y.velocity[1];
jac_data_sur[3][i] = s.Y.velocity[2];
jac_data_sur[4][i] = s.U.E_total;
if (context->is_primitive) {
jac_data_sur[0][i] = s.Y.pressure;
for (int j=0; j<3; j++) jac_data_sur[j+1][i] = s.Y.velocity[j];
jac_data_sur[4][i] = s.Y.temperature;
} else {
jac_data_sur[0][i] = s.U.density;
for (int j=0; j<3; j++) jac_data_sur[j+1][i] = s.U.momentum[j];
jac_data_sur[4][i] = s.U.E_total;
}
for (int j=0; j<6; j++) jac_data_sur[5+j][i] = kmstress[j];
}
return 0;
Expand All @@ -873,6 +885,13 @@ CEED_QFUNCTION(BoundaryIntegral_Jacobian)(void *ctx, CeedInt Q,

const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx;
const bool implicit = context->is_implicit;
State (*StateFromQi)(NewtonianIdealGasContext gas,
const CeedScalar qi[5], const CeedScalar x[3]);
State (*StateFromQi_fwd)(NewtonianIdealGasContext gas,
State s, const CeedScalar dqi[5],
const CeedScalar x[3], const CeedScalar dx[3]);
StateFromQi = context->is_primitive ? &StateFromY : &StateFromU;
StateFromQi_fwd = context->is_primitive ? &StateFromY_fwd : &StateFromU_fwd;

CeedPragmaSIMD
// Quadrature Point Loop
Expand All @@ -888,22 +907,22 @@ CEED_QFUNCTION(BoundaryIntegral_Jacobian)(void *ctx, CeedInt Q,
{q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]}
};

CeedScalar U[5], kmstress[6], dU[5], dx_i[3] = {0.};
for (int j=0; j<5; j++) U[j] = jac_data_sur[j][i];
for (int j=0; j<6; j++) kmstress[j] = jac_data_sur[5+j][i];
for (int j=0; j<3; j++) U[j+1] *= U[0];
for (int j=0; j<5; j++) dU[j] = dq[j][i];
State s = StateFromU(context, U, x_i);
State ds = StateFromU_fwd(context, s, dU, x_i, dx_i);
CeedScalar qi[5], kmstress[6], dqi[5], dx_i[3] = {0.};
for (int j=0; j<5; j++) qi[j] = jac_data_sur[j][i];
for (int j=0; j<6; j++) kmstress[j] = jac_data_sur[5+j][i];
for (int j=0; j<5; j++) dqi[j] = dq[j][i];

State s = StateFromQi(context, qi, x_i);
State ds = StateFromQi_fwd(context, s, dqi, x_i, dx_i);

State grad_ds[3];
for (CeedInt j=0; j<3; j++) {
CeedScalar dx_i[3] = {0}, dUj[5];
CeedScalar dx_i[3] = {0}, dqi_j[5];
for (CeedInt k=0; k<5; k++)
dUj[k] = Grad_dq[0][k][i] * dXdx[0][j] +
Grad_dq[1][k][i] * dXdx[1][j];
dqi_j[k] = Grad_dq[0][k][i] * dXdx[0][j] +
Grad_dq[1][k][i] * dXdx[1][j];
dx_i[j] = 1.;
grad_ds[j] = StateFromU_fwd(context, s, dUj, x_i, dx_i);
grad_ds[j] = StateFromQi_fwd(context, s, dqi_j, x_i, dx_i);
}

CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3];
Expand Down Expand Up @@ -949,14 +968,22 @@ CEED_QFUNCTION(PressureOutflow)(void *ctx, CeedInt Q,
const bool implicit = context->is_implicit;
const CeedScalar P0 = context->P0;

State (*StateFromQi)(NewtonianIdealGasContext gas,
const CeedScalar qi[5], const CeedScalar x[3]);
State (*StateFromQi_fwd)(NewtonianIdealGasContext gas,
State s, const CeedScalar dqi[5],
const CeedScalar x[3], const CeedScalar dx[3]);
StateFromQi = context->is_primitive ? &StateFromY : &StateFromU;
StateFromQi_fwd = context->is_primitive ? &StateFromY_fwd : &StateFromU_fwd;

CeedPragmaSIMD
// Quadrature Point Loop
for (CeedInt i=0; i<Q; i++) {
// Setup
// -- Interp in
const CeedScalar U[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]};
State s = StateFromU(context, U, x_i);
const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
State s = StateFromQi(context, qi, x_i);
s.Y.pressure = P0;

// -- Interp-to-Interp q_data
Expand All @@ -978,12 +1005,12 @@ CEED_QFUNCTION(PressureOutflow)(void *ctx, CeedInt Q,

State grad_s[3];
for (CeedInt j=0; j<3; j++) {
CeedScalar dx_i[3] = {0}, dU[5];
CeedScalar dx_i[3] = {0}, dqi[5];
for (CeedInt k=0; k<5; k++)
dU[k] = Grad_q[0][k][i] * dXdx[0][j] +
Grad_q[1][k][i] * dXdx[1][j];
dqi[k] = Grad_q[0][k][i] * dXdx[0][j] +
Grad_q[1][k][i] * dXdx[1][j];
dx_i[j] = 1.;
grad_s[j] = StateFromU_fwd(context, s, dU, x_i, dx_i);
grad_s[j] = StateFromQi_fwd(context, s, dqi, x_i, dx_i);
}

CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3];
Expand Down Expand Up @@ -1014,20 +1041,23 @@ CEED_QFUNCTION(PressureOutflow)(void *ctx, CeedInt Q,
v[4][i] = -wdetJb * Flux[4];

// Save values for Jacobian
jac_data_sur[0][i] = s.U.density;
jac_data_sur[1][i] = s.Y.velocity[0];
jac_data_sur[2][i] = s.Y.velocity[1];
jac_data_sur[3][i] = s.Y.velocity[2];
jac_data_sur[4][i] = s.U.E_total;
if (context->is_primitive) {
jac_data_sur[0][i] = s.Y.pressure;
for (int j=0; j<3; j++) jac_data_sur[j+1][i] = s.Y.velocity[j];
jac_data_sur[4][i] = s.Y.temperature;
} else {
jac_data_sur[0][i] = s.U.density;
for (int j=0; j<3; j++) jac_data_sur[j+1][i] = s.U.momentum[j];
jac_data_sur[4][i] = s.U.E_total;
}
for (int j=0; j<6; j++) jac_data_sur[5+j][i] = kmstress[j];
} // End Quadrature Point Loop
return 0;
}

// Jacobian for weak-pressure outflow boundary condition
CEED_QFUNCTION(PressureOutflow_Jacobian)(void *ctx, CeedInt Q,
const CeedScalar *const *in,
CeedScalar *const *out) {
const CeedScalar *const *in, CeedScalar *const *out) {
// *INDENT-OFF*
// Inputs
const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0],
Expand All @@ -1042,6 +1072,14 @@ CEED_QFUNCTION(PressureOutflow_Jacobian)(void *ctx, CeedInt Q,
const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx;
const bool implicit = context->is_implicit;

State (*StateFromQi)(NewtonianIdealGasContext gas,
const CeedScalar qi[5], const CeedScalar x[3]);
State (*StateFromQi_fwd)(NewtonianIdealGasContext gas,
State s, const CeedScalar dQi[5],
const CeedScalar x[3], const CeedScalar dx[3]);
StateFromQi = context->is_primitive ? &StateFromY : &StateFromU;
StateFromQi_fwd = context->is_primitive ? &StateFromY_fwd : &StateFromU_fwd;

CeedPragmaSIMD
// Quadrature Point Loop
for (CeedInt i=0; i<Q; i++) {
Expand All @@ -1056,24 +1094,24 @@ CEED_QFUNCTION(PressureOutflow_Jacobian)(void *ctx, CeedInt Q,
{q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]}
};

CeedScalar U[5], kmstress[6], dU[5], dx_i[3] = {0.};
for (int j=0; j<5; j++) U[j] = jac_data_sur[j][i];
for (int j=0; j<6; j++) kmstress[j] = jac_data_sur[5+j][i];
for (int j=0; j<3; j++) U[j+1] *= U[0];
for (int j=0; j<5; j++) dU[j] = dq[j][i];
State s = StateFromU(context, U, x_i);
State ds = StateFromU_fwd(context, s, dU, x_i, dx_i);
CeedScalar qi[5], kmstress[6], dqi[5], dx_i[3] = {0.};
for (int j=0; j<5; j++) qi[j] = jac_data_sur[j][i];
for (int j=0; j<6; j++) kmstress[j] = jac_data_sur[5+j][i];
for (int j=0; j<5; j++) dqi[j] = dq[j][i];

State s = StateFromQi(context, qi, x_i);
State ds = StateFromQi_fwd(context, s, dqi, x_i, dx_i);
s.Y.pressure = context->P0;
ds.Y.pressure = 0.;

State grad_ds[3];
for (CeedInt j=0; j<3; j++) {
CeedScalar dx_i[3] = {0}, dUj[5];
CeedScalar dx_i[3] = {0}, dqi_j[5];
for (CeedInt k=0; k<5; k++)
dUj[k] = Grad_dq[0][k][i] * dXdx[0][j] +
Grad_dq[1][k][i] * dXdx[1][j];
dqi_j[k] = Grad_dq[0][k][i] * dXdx[0][j] +
Grad_dq[1][k][i] * dXdx[1][j];
dx_i[j] = 1.;
grad_ds[j] = StateFromU_fwd(context, s, dUj, x_i, dx_i);
grad_ds[j] = StateFromQi_fwd(context, s, dqi_j, x_i, dx_i);
}

CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3];
Expand Down
2 changes: 1 addition & 1 deletion examples/fluids/qfunctions/newtonian_types.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ struct NewtonianIdealGasContext_ {
CeedScalar ijacobian_time_shift;
CeedScalar P0;
bool is_implicit;
bool primitive;
bool is_primitive;
StabilizationType stabilization;
};

Expand Down
37 changes: 27 additions & 10 deletions examples/fluids/qfunctions/stg_shur14.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -321,11 +321,19 @@ CEED_QFUNCTION(ICsSTG)(void *ctx, CeedInt Q,
for(CeedInt i=0; i<Q; i++) {
InterpolateProfile(X[1][i], u, cij, &eps, &lt, stg_ctx);

q0[0][i] = rho;
q0[1][i] = u[0] * rho;
q0[2][i] = u[1] * rho;
q0[3][i] = u[2] * rho;
q0[4][i] = rho * (0.5 * Dot3(u, u) + cv * theta0);
if (stg_ctx->newtonian_ctx.is_primitive) {
q0[0][i] = P0;
q0[1][i] = u[0];
q0[2][i] = u[1];
q0[3][i] = u[2];
q0[4][i] = theta0;
} else {
q0[0][i] = rho;
q0[1][i] = u[0] * rho;
q0[2][i] = u[1] * rho;
q0[3][i] = u[2] * rho;
q0[4][i] = rho * (0.5 * Dot3(u, u) + cv * theta0);
}
} // End of Quadrature Point Loop
return 0;
}
Expand Down Expand Up @@ -560,11 +568,20 @@ CEED_QFUNCTION(STGShur14_Inflow_StrongQF)(void *ctx, CeedInt Q,
for (CeedInt j=0; j<3; j++) u[j] = ubar[j];
}

bcval[0][i] = scale[i] * rho;
bcval[1][i] = scale[i] * rho * u[0];
bcval[2][i] = scale[i] * rho * u[1];
bcval[3][i] = scale[i] * rho * u[2];
bcval[4][i] = 0.;
if (stg_ctx->newtonian_ctx.is_primitive) {
bcval[0][i] = 0;
bcval[1][i] = scale[i] * u[0];
bcval[2][i] = scale[i] * u[1];
bcval[3][i] = scale[i] * u[2];
bcval[4][i] = scale[i] * theta0;

} else {
bcval[0][i] = scale[i] * rho;
bcval[1][i] = scale[i] * rho * u[0];
bcval[2][i] = scale[i] * rho * u[1];
bcval[3][i] = scale[i] * rho * u[2];
bcval[4][i] = 0.;
}
}
return 0;
}
Expand Down
2 changes: 1 addition & 1 deletion examples/fluids/src/setupdm.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,7 @@ PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree,
CHKERRQ(ierr);

if (use_strongstg) {
ierr = SetupStrongSTG(dm, bc, problem); CHKERRQ(ierr);
ierr = SetupStrongSTG(dm, bc, problem, phys); CHKERRQ(ierr);
}
}

Expand Down

0 comments on commit 39d18f8

Please sign in to comment.