Skip to content

Commit

Permalink
Define a new qfunction for converting the final solution from
Browse files Browse the repository at this point in the history 
conservative to primitive variables.
It specifically adds support for the blasius problem.
  • Loading branch information
@LeilaGhaffari
LeilaGhaffari committed Jan 3, 2023
1 parent 60d0891 commit 447e215
Show file tree
Hide file tree
Showing 7 changed files with 166 additions and 112 deletions.
4 changes: 2 additions & 2 deletions examples/fluids/blasius.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,13 +11,13 @@ output_freq: 10
## Linear Settings:
degree: 1
dm_plex_box_faces: 40,60,1
platemesh_nDelta: 45
platemesh_Ndelta: 45

# # Quadratic Settings:
# degree: 2
# dm_plex_box_faces: 20,30,1
# platemesh:
# nDelta: 22
# Ndelta: 22
# growth: 1.1664 # 1.08^2

stab: 'supg'
Expand Down
11 changes: 7 additions & 4 deletions examples/fluids/navierstokes.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,9 +110,9 @@ struct CeedData_private {
CeedVector x_coord, q_data;
CeedBasis basis_x, basis_xc, basis_q, basis_x_sur, basis_q_sur, basis_xc_sur;
CeedElemRestriction elem_restr_x, elem_restr_q, elem_restr_qd_i;
CeedOperator op_setup_vol, op_ics;
CeedOperator op_setup_vol, op_ics, op_to_prim;
CeedQFunction qf_setup_vol, qf_ics, qf_rhs_vol, qf_ifunction_vol, qf_setup_sur, qf_apply_inflow, qf_apply_inflow_jacobian, qf_apply_outflow,
qf_apply_outflow_jacobian, qf_apply_freestream, qf_apply_freestream_jacobian;
qf_apply_outflow_jacobian, qf_apply_freestream, qf_apply_freestream_jacobian, qf_to_prim;
};

// PETSc user data
Expand Down Expand Up @@ -188,7 +188,7 @@ struct ProblemData_private {
CeedInt dim, q_data_size_vol, q_data_size_sur, jac_data_size_sur;
CeedScalar dm_scale;
ProblemQFunctionSpec setup_vol, setup_sur, ics, apply_vol_rhs, apply_vol_ifunction, apply_vol_ijacobian, apply_inflow, apply_outflow,
apply_freestream, apply_inflow_jacobian, apply_outflow_jacobian, apply_freestream_jacobian;
apply_freestream, apply_inflow_jacobian, apply_outflow_jacobian, apply_freestream_jacobian, to_prim;
bool non_zero_time;
PetscErrorCode (*bc)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
void *bc_ctx;
Expand Down Expand Up @@ -294,7 +294,10 @@ PetscErrorCode DMPlexInsertBoundaryValues_NS(DM dm, PetscBool insert_essential,
PetscErrorCode RegressionTests_NS(AppCtx app_ctx, Vec Q);

// Get error for problems with exact solutions
PetscErrorCode GetError_NS(CeedData ceed_data, DM dm, User user, Vec Q, PetscScalar final_time);
PetscErrorCode ConvertStateVar_NS(DM dm, CeedData ceed_data, User user, Vec Qsource_soln_loc, Vec Qtarget_loc, Vec Qtarget);
PetscErrorCode ComputeErrorComponent_NS(Vec Q_true_loc, Vec Q_soln_loc, NormType norm_type, MPI_Op norm_reduce, PetscReal norm_true[5],
PetscReal norm_error[5], PetscReal rel_error[5]);
PetscErrorCode GetError_NS(CeedData ceed_data, DM dm, User user, ProblemData *problem, Vec Qsource_soln, PetscScalar final_time);

// Post-processing
PetscErrorCode PostProcess_NS(TS ts, CeedData ceed_data, DM dm, ProblemData *problem, User user, Vec Q, PetscScalar final_time);
Expand Down
1 change: 1 addition & 0 deletions examples/fluids/problems/blasius.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -247,6 +247,7 @@ PetscErrorCode NS_BLASIUS(ProblemData *problem, DM dm, void *ctx, SimpleBC bc) {
// ------------------------------------------------------
problem->ics.qfunction = ICsBlasius;
problem->ics.qfunction_loc = ICsBlasius_loc;
problem->non_zero_time = PETSC_TRUE;

CeedScalar U_inf = 40; // m/s
CeedScalar T_inf = 288.; // K
Expand Down
3 changes: 3 additions & 0 deletions examples/fluids/problems/newtonian.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,6 +88,8 @@ PetscErrorCode NS_NEWTONIAN_IG(ProblemData *problem, DM dm, void *ctx, SimpleBC
problem->setup_vol.qfunction_loc = Setup_loc;
problem->setup_sur.qfunction = SetupBoundary;
problem->setup_sur.qfunction_loc = SetupBoundary_loc;
problem->to_prim.qfunction = Newtonian_ToPrimitive;
problem->to_prim.qfunction_loc = Newtonian_ToPrimitive_loc;
problem->bc = NULL;
problem->bc_ctx = NULL;
problem->non_zero_time = PETSC_FALSE;
Expand Down Expand Up @@ -298,6 +300,7 @@ PetscErrorCode NS_NEWTONIAN_IG(ProblemData *problem, DM dm, void *ctx, SimpleBC
CeedQFunctionContextReferenceCopy(newtonian_ig_context, &problem->apply_inflow_jacobian.qfunction_context);
CeedQFunctionContextReferenceCopy(newtonian_ig_context, &problem->apply_outflow.qfunction_context);
CeedQFunctionContextReferenceCopy(newtonian_ig_context, &problem->apply_outflow_jacobian.qfunction_context);
CeedQFunctionContextReferenceCopy(newtonian_ig_context, &problem->to_prim.qfunction_context);

if (unit_tests) {
PetscCall(UnitTests_Newtonian(user, newtonian_ig_ctx));
Expand Down
24 changes: 24 additions & 0 deletions examples/fluids/qfunctions/newtonian.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,30 @@ CEED_QFUNCTION(ICsNewtonianIG_Conserv)(void *ctx, CeedInt Q, const CeedScalar *c
return ICsNewtonianIG(ctx, Q, in, out, StateToU);
}

// *****************************************************************************
// This QFunction converts the conservative variables to primitive variables.
// *****************************************************************************
CEED_QFUNCTION(Newtonian_ToPrimitive)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
// Inputs
const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
const CeedScalar(*q_cons)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1];

// Outputs
CeedScalar(*q_prim)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];

NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx;
// Quadrature Point Loop
CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
const CeedScalar x[3] = {X[0][i], X[1][i], X[2][i]};
CeedScalar q_c[5] = {0}, q_p[5] = {0};
for (CeedInt j = 0; j < 5; j++) q_c[j] = q_cons[j][i];
State s = StateFromU(context, q_c, x);
StateToY(context, s, q_p);
for (CeedInt j = 0; j < 5; j++) q_prim[j][i] = q_p[j];
} // End of Quadrature Point Loop
return 0;
}

// *****************************************************************************
// This QFunction implements the following formulation of Navier-Stokes with explicit time stepping method
//
Expand Down
217 changes: 111 additions & 106 deletions examples/fluids/src/misc.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -131,123 +131,127 @@ PetscErrorCode RegressionTests_NS(AppCtx app_ctx, Vec Q) {

PetscFunctionReturn(0);
}
PetscErrorCode GetConservativeFromPrimitive_NS(Vec Q_prim, Vec Q_cons) {
PetscInt Q_size;
PetscInt num_comp = 5;
PetscScalar cv = 2.5;
const PetscScalar *Y;
PetscScalar *U;

PetscFunctionBegin;
PetscCall(VecGetSize(Q_prim, &Q_size));
PetscCall(VecGetArrayRead(Q_prim, &Y));
PetscCall(VecGetArrayWrite(Q_cons, &U));
for (PetscInt i = 0; i < Q_size; i += num_comp) {
// Primitive variables
PetscScalar P = Y[i + 0], u[3] = {Y[i + 1], Y[i + 2], Y[i + 3]}, T = Y[i + 4];
PetscScalar rho = P / T;
U[i + 0] = rho;
for (int j = 0; j < 3; j++) U[i + 1 + j] = rho * u[j];
U[i + 4] = rho * (cv * T + .5 * (u[0] * u[0] + u[1] * u[1] + u[2] * u[2]));
}
PetscCall(VecRestoreArrayRead(Q_prim, &Y));
PetscCall(VecRestoreArrayWrite(Q_cons, &U));
PetscErrorCode ConvertStateVar_NS(DM dm, CeedData ceed_data, User user, Vec Qsource_loc, Vec Qtarget_loc, Vec Qtarget) {
PetscFunctionBeginUser;
// CEED Restriction
CeedVector q_source_ceed, q_target_ceed;
CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &q_source_ceed, NULL);
CeedElemRestrictionCreateVector(ceed_data->elem_restr_q, &q_target_ceed, NULL);

// Place PETSc vectors in CEED vectors
PetscScalar *q_source, *q_target;
PetscMemType q_source_mem_type, q_target_mem_type;
PetscCall(VecGetArrayReadAndMemType(Qsource_loc, (const PetscScalar **)&q_source, &q_source_mem_type));
CeedVectorSetArray(q_source_ceed, MemTypeP2C(q_source_mem_type), CEED_USE_POINTER, q_source);
PetscCall(VecZeroEntries(Qtarget_loc));
PetscCall(VecGetArrayAndMemType(Qtarget_loc, (PetscScalar **)&q_target, &q_target_mem_type));
CeedVectorSetArray(q_target_ceed, MemTypeP2C(q_target_mem_type), CEED_USE_POINTER, q_target);

// Apply CEED Operator
CeedOperatorApply(ceed_data->op_to_prim, q_source_ceed, q_target_ceed, CEED_REQUEST_IMMEDIATE);

// Restore vectors
CeedVectorTakeArray(q_source_ceed, MemTypeP2C(q_source_mem_type), NULL);
PetscCall(VecRestoreArrayReadAndMemType(Qsource_loc, (const PetscScalar **)&q_source));
CeedVectorTakeArray(q_target_ceed, MemTypeP2C(q_target_mem_type), NULL);
PetscCall(VecRestoreArrayReadAndMemType(Qtarget_loc, (const PetscScalar **)&q_target));

// Local-to-Global
PetscCall(VecZeroEntries(Qtarget));
PetscCall(DMLocalToGlobal(dm, Qtarget_loc, ADD_VALUES, Qtarget));

// Cleanup
CeedVectorDestroy(&q_source_ceed);
CeedVectorDestroy(&q_target_ceed);

PetscFunctionReturn(0);
}

PetscErrorCode GetPrimitiveFromConservative_NS(Vec Q_cons, Vec Q_prim) {
PetscInt Q_size;
PetscInt num_comp = 5;
PetscScalar cv = 2.5;
const PetscScalar *U;
PetscScalar *Y;

PetscFunctionBegin;
PetscCall(VecGetSize(Q_cons, &Q_size));
PetscCall(VecGetArrayRead(Q_cons, &U));
PetscCall(VecGetArrayWrite(Q_prim, &Y));
for (PetscInt i = 0; i < Q_size; i += num_comp) {
// Conservative variables
PetscScalar rho = U[i + 0], M[3] = {U[i + 1], U[i + 2], U[i + 3]}, E = U[i + 4];
// Primitive variables
for (int j = 0; j < 3; j++) Y[i + 1 + j] = M[j] / rho;
Y[i + 4] = (E - 0.5 * (M[0] * M[0] + M[1] * M[1] + M[2] * M[2]) / rho) / (cv * rho);
Y[i + 0] = (E - 0.5 * (M[0] * M[0] + M[1] * M[1] + M[2] * M[2]) / rho) / cv;
PetscErrorCode ComputeErrorComponent_NS(Vec Q_true_loc, Vec Q_soln_loc, NormType norm_type, MPI_Op norm_reduce, PetscReal norm_true[5],
PetscReal norm_error[5], PetscReal rel_error[5]) {
// Get the norm of each component
PetscCall(VecStrideNormAll(Q_true_loc, norm_type, norm_true));
PetscCallMPI(MPI_Allreduce(MPI_IN_PLACE, norm_true, 5, MPIU_REAL, norm_reduce, PETSC_COMM_WORLD));
// Compute the error of each component
PetscCall(VecAXPY(Q_soln_loc, -1.0, Q_true_loc));
PetscCall(VecStrideNormAll(Q_soln_loc, norm_type, norm_error));
PetscCallMPI(MPI_Allreduce(MPI_IN_PLACE, norm_error, 5, MPIU_REAL, norm_reduce, PETSC_COMM_WORLD));
for (int i = 0; i < 5; i++) {
rel_error[i] = norm_error[i] / norm_true[i];
}
PetscCall(VecRestoreArrayRead(Q_cons, &U));
PetscCall(VecRestoreArrayWrite(Q_prim, &Y));
PetscFunctionReturn(0);
}

PetscErrorCode GetError_NS(CeedData ceed_data, DM dm, User user, Vec Q, PetscScalar final_time) {
PetscInt loc_nodes;
Vec Q_exact, Q_exact_loc;
Vec Q_loc = Q;
PetscReal rel_error, norm_error, norm_exact;
// Get error for problems with true solutions
PetscErrorCode GetError_NS(CeedData ceed_data, DM dm, User user, ProblemData *problem, Vec Qsource_soln, PetscScalar final_time) {
Vec Qsource_soln_loc, Qsource_true, Qsource_true_loc;
PetscFunctionBegin;

// Get exact solution at final time
PetscCall(DMCreateGlobalVector(dm, &Q_exact));
PetscCall(DMGetLocalVector(dm, &Q_exact_loc));
PetscCall(VecGetSize(Q_exact_loc, &loc_nodes));
PetscCall(ICs_FixMultiplicity(dm, ceed_data, user, Q_exact_loc, Q_exact, final_time));

if (user->phys->state_var == STATEVAR_PRIMITIVE) { // If primitive, get conservative
Vec Q_cons, Q_cons_exact;
PetscReal rel_error_cons, norm_error_cons, norm_exact_cons;

PetscCall(VecDuplicate(Q_loc, &Q_cons));
PetscCall(GetConservativeFromPrimitive_NS(Q_loc, Q_cons));
PetscCall(VecDuplicate(Q_exact, &Q_cons_exact));
PetscCall(GetConservativeFromPrimitive_NS(Q_exact, Q_cons_exact));

// Get |exact solution - obtained solution|
PetscCall(VecNorm(Q_exact, NORM_1, &norm_exact));
PetscCall(VecAXPY(Q_loc, -1.0, Q_exact));
PetscCall(VecNorm(Q_loc, NORM_1, &norm_error));

PetscCall(VecNorm(Q_cons_exact, NORM_1, &norm_exact_cons));
PetscCall(VecAXPY(Q_cons, -1.0, Q_cons_exact));
PetscCall(VecNorm(Q_cons, NORM_1, &norm_error_cons));

// Compute relative error
rel_error = norm_error / norm_exact;
rel_error_cons = norm_error_cons / norm_exact_cons;

// Output relative error
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Relative Error in primitive variables: %g\n", (double)rel_error));
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Relative Error in conservative variables:: %g\n", (double)rel_error_cons));

} else {
Vec Q_prim, Q_prim_exact;
PetscReal rel_error_prim, norm_error_prim, norm_exact_prim;

PetscCall(VecDuplicate(Q_loc, &Q_prim));
PetscCall(GetPrimitiveFromConservative_NS(Q_loc, Q_prim));
PetscCall(VecDuplicate(Q_exact, &Q_prim_exact));
PetscCall(GetPrimitiveFromConservative_NS(Q_exact, Q_prim_exact));

// Get |exact solution - obtained solution|
PetscCall(VecNorm(Q_exact, NORM_1, &norm_exact));
PetscCall(VecAXPY(Q_loc, -1.0, Q_exact));
PetscCall(VecNorm(Q_loc, NORM_1, &norm_error));

PetscCall(VecNorm(Q_prim_exact, NORM_1, &norm_exact_prim));
PetscCall(VecAXPY(Q_prim, -1.0, Q_prim_exact));
PetscCall(VecNorm(Q_prim, NORM_1, &norm_error_prim));

// Compute relative error
rel_error = norm_error / norm_exact;
rel_error_prim = norm_error_prim / norm_exact_prim;

// Output relative error
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Relative Error in primitive variables: %g\n", (double)rel_error_prim));
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Relative Error in conservative variables: %g\n", (double)rel_error));
// Update time for evaluation
if (user->phys->ics_time_label) CeedOperatorContextSetDouble(ceed_data->op_ics, user->phys->ics_time_label, &final_time);

// Define norm types
NormType norm_type = NORM_MAX;
MPI_Op norm_reduce = MPI_MAX;

// Compute the error in the source state variable
// -- Get the local values of the final solution
PetscCall(DMGetLocalVector(dm, &Qsource_soln_loc));
PetscCall(DMGlobalToLocal(dm, Qsource_soln, INSERT_VALUES, Qsource_soln_loc));
PetscCall(DMPlexInsertBoundaryValues(dm, PETSC_TRUE, Qsource_soln_loc, final_time, NULL, NULL, NULL));
// -- Get true solution at final time
PetscCall(DMCreateGlobalVector(dm, &Qsource_true));
PetscCall(DMGetLocalVector(dm, &Qsource_true_loc));
PetscCall(ICs_FixMultiplicity(dm, ceed_data, user, Qsource_true_loc, Qsource_true, final_time));
PetscCall(DMPlexInsertBoundaryValues(dm, PETSC_TRUE, Qsource_true_loc, final_time, NULL, NULL, NULL));
// -- Compute error
PetscReal norm_true_source[5], norm_error_source[5], rel_error_source[5];
PetscCall(
ComputeErrorComponent_NS(Qsource_true_loc, Qsource_soln_loc, norm_type, norm_reduce, norm_true_source, norm_error_source, rel_error_source));
// -- Report error
const char *state_var_source = "Conservative";
if (user->phys->state_var == STATEVAR_PRIMITIVE) {
state_var_source = "Primitive";
}
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nRelative Error (absolute error norm, true solution norm):\n"));
PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %s variables:\n", state_var_source));
for (int i = 0; i < 5; i++) {
PetscCall(
PetscPrintf(PETSC_COMM_WORLD, "\tComponent %d: %g (%g, %g)\n", i, (double)rel_error_source[i], norm_error_source[i], norm_true_source[i]));
}

// Convert solution to primitive variables, if possible, and compute the error
if (user->phys->state_var != STATEVAR_PRIMITIVE && problem->to_prim.qfunction) {
Vec Qprim_soln, Qprim_soln_loc, Qprim_true, Qprim_true_loc;
// -- Convert the final solution from the source state variable set to the prim state variable set
PetscCall(DMCreateGlobalVector(dm, &Qprim_soln));
PetscCall(DMGetLocalVector(dm, &Qprim_soln_loc));
PetscCall(ConvertStateVar_NS(dm, ceed_data, user, Qsource_soln_loc, Qprim_soln_loc, Qprim_soln));
// -- Convert the true values from the source state variable set to the prim state variable set
PetscCall(DMCreateGlobalVector(dm, &Qprim_true));
PetscCall(DMGetLocalVector(dm, &Qprim_true_loc));
PetscCall(ConvertStateVar_NS(dm, ceed_data, user, Qsource_true_loc, Qprim_true_loc, Qprim_true));
// -- Compute error
PetscReal norm_true_prim[5], norm_error_prim[5], rel_error_prim[5];
PetscCall(ComputeErrorComponent_NS(Qprim_true_loc, Qprim_soln_loc, norm_type, norm_reduce, norm_true_prim, norm_error_prim, rel_error_prim));
// -- Report error
const char *state_var_prim = "Primitive";
PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %s variables:\n", state_var_prim));
for (int i = 0; i < 5; i++) {
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\tComponent %d: %g (%g, %g)\n", i, (double)rel_error_prim[i], norm_error_prim[i], norm_true_prim[i]));
}
// Cleanup
PetscCall(DMRestoreLocalVector(dm, &Qprim_true_loc));
PetscCall(DMRestoreLocalVector(dm, &Qprim_soln_loc));
PetscCall(VecDestroy(&Qprim_true));
PetscCall(VecDestroy(&Qprim_soln));
}

// Cleanup
PetscCall(DMRestoreLocalVector(dm, &Q_exact_loc));
PetscCall(VecDestroy(&Q_exact));
PetscCall(DMRestoreLocalVector(dm, &Qsource_true_loc));
PetscCall(DMRestoreLocalVector(dm, &Qsource_soln_loc));
PetscCall(VecDestroy(&Qsource_true));

PetscFunctionReturn(0);
}
Expand All @@ -259,13 +263,14 @@ PetscErrorCode PostProcess_NS(TS ts, CeedData ceed_data, DM dm, ProblemData *pro

// Print relative error
if (problem->non_zero_time && !user->app_ctx->test_mode) {
PetscCall(GetError_NS(ceed_data, dm, user, Q, final_time));
PetscCall(GetError_NS(ceed_data, dm, user, problem, Q, final_time));
}

// Print final time and number of steps
PetscCall(TSGetStepNumber(ts, &steps));
if (!user->app_ctx->test_mode) {
PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Time integrator took %" PetscInt_FMT " time steps to reach final time %g\n", steps, (double)final_time));
PetscCall(
PetscPrintf(PETSC_COMM_WORLD, "\nTime integrator took %" PetscInt_FMT " time steps to reach final time %g\n", steps, (double)final_time));
}

// Output numerical values from command line
Expand Down
18 changes: 18 additions & 0 deletions examples/fluids/src/setuplibceed.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -260,6 +260,16 @@ PetscErrorCode SetupLibceed(Ceed ceed, CeedData ceed_data, DM dm, User user, App
CeedQFunctionAddInput(ceed_data->qf_ics, "qdata", q_data_size_vol, CEED_EVAL_NONE);
CeedQFunctionAddOutput(ceed_data->qf_ics, "q0", num_comp_q, CEED_EVAL_NONE);

// -- Create QFunction for converting the state variable set
if (problem->non_zero_time && !user->app_ctx->test_mode && problem->to_prim.qfunction && user->phys->state_var != STATEVAR_PRIMITIVE) {
CeedQFunctionCreateInterior(ceed, 1, problem->to_prim.qfunction, problem->to_prim.qfunction_loc, &ceed_data->qf_to_prim);
CeedQFunctionSetContext(ceed_data->qf_to_prim, problem->to_prim.qfunction_context);
CeedQFunctionContextDestroy(&problem->to_prim.qfunction_context);
CeedQFunctionAddInput(ceed_data->qf_to_prim, "x", num_comp_x, CEED_EVAL_INTERP);
CeedQFunctionAddInput(ceed_data->qf_to_prim, "q_cons", num_comp_q, CEED_EVAL_INTERP);
CeedQFunctionAddOutput(ceed_data->qf_to_prim, "q_prim", num_comp_q, CEED_EVAL_NONE);
}

// -- Create QFunction for RHS
if (problem->apply_vol_rhs.qfunction) {
CeedQFunctionCreateInterior(ceed, 1, problem->apply_vol_rhs.qfunction, problem->apply_vol_rhs.qfunction_loc, &ceed_data->qf_rhs_vol);
Expand Down Expand Up @@ -353,6 +363,14 @@ PetscErrorCode SetupLibceed(Ceed ceed, CeedData ceed_data, DM dm, User user, App
CeedOperatorSetField(ceed_data->op_ics, "q0", ceed_data->elem_restr_q, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
CeedOperatorContextGetFieldLabel(ceed_data->op_ics, "evaluation time", &user->phys->ics_time_label);

// -- Create CEED operator for converting the state variable set
if (ceed_data->qf_to_prim) {
CeedOperatorCreate(ceed, ceed_data->qf_to_prim, NULL, NULL, &ceed_data->op_to_prim);
CeedOperatorSetField(ceed_data->op_to_prim, "x", ceed_data->elem_restr_x, ceed_data->basis_x, ceed_data->x_coord);
CeedOperatorSetField(ceed_data->op_to_prim, "q_cons", ceed_data->elem_restr_q, ceed_data->basis_q, CEED_VECTOR_ACTIVE);
CeedOperatorSetField(ceed_data->op_to_prim, "q_prim", ceed_data->elem_restr_q, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
}

// Create CEED operator for RHS
if (ceed_data->qf_rhs_vol) {
CeedOperator op;
Expand Down

0 comments on commit 447e215

Please sign in to comment.