diff --git a/examples/fluids/README.md b/examples/fluids/README.md index 16df8a11e2..0a7bce7e9f 100644 --- a/examples/fluids/README.md +++ b/examples/fluids/README.md @@ -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 ::: diff --git a/examples/fluids/problems/blasius.c b/examples/fluids/problems/blasius.c index 4243d66b74..b990bc2f77 100644 --- a/examples/fluids/problems/blasius.c +++ b/examples/fluids/problems/blasius.c @@ -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]]; diff --git a/examples/fluids/problems/newtonian.c b/examples/fluids/problems/newtonian.c index bd07858d76..e1bee405fa 100644 --- a/examples/fluids/problems/newtonian.c +++ b/examples/fluids/problems/newtonian.c @@ -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; @@ -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); diff --git a/examples/fluids/problems/stg_shur14.c b/examples/fluids/problems/stg_shur14.c index e267d6ba79..8cae0a5f3d 100644 --- a/examples/fluids/problems/stg_shur14.c +++ b/examples/fluids/problems/stg_shur14.c @@ -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) { diff --git a/examples/fluids/problems/stg_shur14.h b/examples/fluids/problems/stg_shur14.h index af7fb58874..e7be94741a 100644 --- a/examples/fluids/problems/stg_shur14.h +++ b/examples/fluids/problems/stg_shur14.h @@ -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, diff --git a/examples/fluids/qfunctions/channel.h b/examples/fluids/qfunctions/channel.h index 96306663b4..b877d98841 100644 --- a/examples/fluids/qfunctions/channel.h +++ b/examples/fluids/qfunctions/channel.h @@ -83,7 +83,7 @@ CEED_QFUNCTION(ICsChannel)(void *ctx, CeedInt Q, for (CeedInt i=0; inewtonian_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]; diff --git a/examples/fluids/qfunctions/densitycurrent.h b/examples/fluids/qfunctions/densitycurrent.h index 8f503a8692..22a10d99ae 100644 --- a/examples/fluids/qfunctions/densitycurrent.h +++ b/examples/fluids/qfunctions/densitycurrent.h @@ -152,7 +152,7 @@ CEED_QFUNCTION(ICsDC)(void *ctx, CeedInt Q, for (CeedInt i=0; inewtonian_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]; diff --git a/examples/fluids/qfunctions/newtonian.h b/examples/fluids/qfunctions/newtonian.h index ffd03d39b8..9fddd27555 100644 --- a/examples/fluids/qfunctions/newtonian.h +++ b/examples/fluids/qfunctions/newtonian.h @@ -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; iis_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; @@ -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 @@ -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]; @@ -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; iis_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; @@ -1026,8 +1057,7 @@ CEED_QFUNCTION(PressureOutflow)(void *ctx, CeedInt Q, // 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], @@ -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; iP0; 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]; diff --git a/examples/fluids/qfunctions/newtonian_types.h b/examples/fluids/qfunctions/newtonian_types.h index 5b5d3af457..51d43f5877 100644 --- a/examples/fluids/qfunctions/newtonian_types.h +++ b/examples/fluids/qfunctions/newtonian_types.h @@ -50,7 +50,7 @@ struct NewtonianIdealGasContext_ { CeedScalar ijacobian_time_shift; CeedScalar P0; bool is_implicit; - bool primitive; + bool is_primitive; StabilizationType stabilization; }; diff --git a/examples/fluids/qfunctions/stg_shur14.h b/examples/fluids/qfunctions/stg_shur14.h index c92710e2a4..6201f43afd 100644 --- a/examples/fluids/qfunctions/stg_shur14.h +++ b/examples/fluids/qfunctions/stg_shur14.h @@ -321,11 +321,19 @@ CEED_QFUNCTION(ICsSTG)(void *ctx, CeedInt Q, for(CeedInt i=0; inewtonian_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; } @@ -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; } diff --git a/examples/fluids/src/setupdm.c b/examples/fluids/src/setupdm.c index 6659686f5b..26fa2c349a 100644 --- a/examples/fluids/src/setupdm.c +++ b/examples/fluids/src/setupdm.c @@ -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); } }