Define a new qfunction for converting the final solution from conserv…

…ative to primitive variables.

examples/fluids/blasius.yaml

@@ -13,13 +13,13 @@ checkpoint_interval: 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'

examples/fluids/navierstokes.h

@@ -112,9 +112,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
@@ -190,7 +190,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;
@@ -296,7 +296,9 @@ 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, 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);

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

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));

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
 //

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, PetscReal norm_true[5], PetscReal norm_error[5], PetscReal rel_error[5]) {
+  // Define norm types
+  NormType norm_type   = NORM_MAX;
+  MPI_Op   norm_reduce = MPI_MAX;
+
+  // 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);
+
+  // 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_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_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));
+    CeedQFunctionDestroy(&ceed_data->qf_to_prim);
+    CeedOperatorDestroy(&ceed_data->op_to_prim);
   }
 
   // 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

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;