diff --git a/examples/fluids/qfunctions/newtonian.h b/examples/fluids/qfunctions/newtonian.h index f23b65966b..da362cf79e 100644 --- a/examples/fluids/qfunctions/newtonian.h +++ b/examples/fluids/qfunctions/newtonian.h @@ -88,10 +88,13 @@ CEED_QFUNCTION(Newtonian_L2Error)(void *ctx, CeedInt Q, const CeedScalar *const } // ***************************************************************************** -// This QFunction computes the L2 error in primitive variables if -// the simulation is run in conservative variables. +// This QFunction computes the L2 error for each component in the source +// state variable and convert them to the target state variable. +// +// Source state variable: STATEVAR_CONSERVATIVE or STATEVAR_PRIMITIVE +// Target state variable: STATEVAR_PRIMITIVE or STATEVAR_CONSERVATIVE // ***************************************************************************** -CEED_QFUNCTION(Newtonian_L2ErrorAsPrimitive)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { +CEED_QFUNCTION_HELPER int Newtonian_ConvertL2Error(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_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; @@ -105,50 +108,33 @@ CEED_QFUNCTION(Newtonian_L2ErrorAsPrimitive)(void *ctx, CeedInt Q, const CeedSca // 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_t[5] = {0}, q_c_s[5] = {0}, q_p_t[5] = {0}, q_p_s[5] = {0}; + CeedScalar q_s_t[5] = {0}, q_s_s[5] = {0}, q_t_t[5] = {0}, q_t_s[5] = {0}; + State s_s, s_t; for (CeedInt j = 0; j < 5; j++) { - q_c_s[j] = q_source_soln[j][i]; - q_c_t[j] = q_source_true[j][i]; // Automatically provided in the same variables as the solution. + q_s_s[j] = q_source_soln[j][i]; + q_s_t[j] = q_source_true[j][i]; // Automatically provided in the same variables as the solution. + } + if (context->state_var == STATEVAR_CONSERVATIVE) { + s_s = StateFromU(context, q_s_s, x); + StateToY(context, s_s, q_t_s); + s_t = StateFromU(context, q_s_t, x); + StateToY(context, s_t, q_t_t); + } else { + s_s = StateFromY(context, q_s_s, x); + StateToU(context, s_s, q_t_s); + s_t = StateFromY(context, q_s_t, x); + StateToU(context, s_t, q_t_t); } - State s_t = StateFromU(context, q_c_t, x); - StateToY(context, s_t, q_p_t); - State s_s = StateFromU(context, q_c_s, x); - StateToY(context, s_s, q_p_s); - for (CeedInt j = 0; j < 5; j++) q_target_error[j][i] = q_data[0][i] * Square(q_p_t[j] - q_p_s[j]); + for (CeedInt j = 0; j < 5; j++) q_target_error[j][i] = q_data[0][i] * Square(q_t_t[j] - q_t_s[j]); } // End of Quadrature Point Loop return 0; } -// ***************************************************************************** -// This QFunction computes the L2 error in conservative variables if -// the simulation is run in primitive variables. -// ***************************************************************************** +CEED_QFUNCTION(Newtonian_L2ErrorAsPrimitive)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { + return Newtonian_ConvertL2Error(ctx, Q, in, out); +} CEED_QFUNCTION(Newtonian_L2ErrorAsConservative)(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_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; - const CeedScalar(*q_source_true)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; - const CeedScalar(*q_source_soln)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; - - // Outputs - CeedScalar(*q_target_error)[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_p_t[5] = {0}, q_p_s[5] = {0}, q_c_t[5] = {0}, q_c_s[5] = {0}; - for (CeedInt j = 0; j < 5; j++) { - q_p_s[j] = q_source_soln[j][i]; - q_p_t[j] = q_source_true[j][i]; // Automatically provided in the same variables as the solution. - } - State s_t = StateFromY(context, q_p_t, x); - StateToU(context, s_t, q_c_t); - State s_s = StateFromY(context, q_p_s, x); - StateToU(context, s_s, q_c_s); - for (CeedInt j = 0; j < 5; j++) q_target_error[j][i] = q_data[0][i] * Square(q_c_t[j] - q_c_s[j]); - } // End of Quadrature Point Loop - return 0; + return Newtonian_ConvertL2Error(ctx, Q, in, out); } // *****************************************************************************