diff --git a/examples/fluids/blasius.yaml b/examples/fluids/blasius.yaml index db6c4ed8a2..624174fca0 100644 --- a/examples/fluids/blasius.yaml +++ b/examples/fluids/blasius.yaml @@ -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' diff --git a/examples/fluids/navierstokes.h b/examples/fluids/navierstokes.h index 28587e98b2..2706fa726d 100644 --- a/examples/fluids/navierstokes.h +++ b/examples/fluids/navierstokes.h @@ -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 @@ -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; @@ -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); diff --git a/examples/fluids/problems/blasius.c b/examples/fluids/problems/blasius.c index 8da08cda02..86e95a2ee7 100644 --- a/examples/fluids/problems/blasius.c +++ b/examples/fluids/problems/blasius.c @@ -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 diff --git a/examples/fluids/problems/newtonian.c b/examples/fluids/problems/newtonian.c index 65e557503e..3309768c06 100644 --- a/examples/fluids/problems/newtonian.c +++ b/examples/fluids/problems/newtonian.c @@ -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; @@ -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)); diff --git a/examples/fluids/qfunctions/newtonian.h b/examples/fluids/qfunctions/newtonian.h index d8b99f3ea1..3f1f13a348 100644 --- a/examples/fluids/qfunctions/newtonian.h +++ b/examples/fluids/qfunctions/newtonian.h @@ -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 // diff --git a/examples/fluids/src/misc.c b/examples/fluids/src/misc.c index a0ed7a8553..3ac06a6e64 100644 --- a/examples/fluids/src/misc.c +++ b/examples/fluids/src/misc.c @@ -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); } @@ -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 diff --git a/examples/fluids/src/setuplibceed.c b/examples/fluids/src/setuplibceed.c index 1cfe6ddb7b..2972a2fad6 100644 --- a/examples/fluids/src/setuplibceed.c +++ b/examples/fluids/src/setuplibceed.c @@ -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); @@ -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;