Merge pull request CEED#1278 from CEED/jrwrigh/blasius-simplex

fluids: Fix STG and Blasius for non-box meshes

examples/fluids/README.md

@@ -160,6 +160,10 @@ The following options are common among all problem types:
   - Viewer for the force on each no-slip wall, e.g., `ascii:force.csv:ascii_csv` to write a CSV file.
   -
 
+* - `-mesh_transform`
+  - Transform the mesh, usually for an initial box mesh.
+  - `none`
+
 * - `-snes_view`
   - View PETSc `SNES` nonlinear solver configuration
   -
@@ -963,6 +967,11 @@ The Blasius problem has the following command-line options in addition to the Ne
   - `1.01E5`
   - `Pa`
 
+* - `-platemesh_modify_mesh`
+  - Whether to modify the mesh using the given options below.
+  - `false`
+  -
+
 * - `-platemesh_refine_height`
   - Height at which `-platemesh_Ndelta` number of elements should refined into
   - `5.9E-4`
@@ -983,16 +992,16 @@ The Blasius problem has the following command-line options in addition to the Ne
   - `5`
   - `degrees`
 
-* - `-stg_use`
-  - Whether to use stg for the inflow conditions
-  - `false`
-  -
-
 * - `-platemesh_y_node_locs_path`
   - Path to file with y node locations. If empty, will use mesh warping instead.
   - `""`
   -
 
+* - `-stg_use`
+  - Whether to use STG for the inflow conditions
+  - `false`
+  -
+
 * - `-n_chebyshev`
   - Number of Chebyshev terms
   - `20`

examples/fluids/blasius.yaml

@@ -13,12 +13,14 @@ checkpoint_interval: 10
 ## Linear Settings:
 degree: 1
 dm_plex_box_faces: 40,60,1
+mesh_transform: platemesh
 platemesh_nDelta: 45
 
 # # Quadratic Settings:
 # degree: 2
 # dm_plex_box_faces: 20,30,1
 # platemesh:
+#   modify_mesh: true
 #   nDelta: 22
 #   growth: 1.1664 # 1.08^2
 

examples/fluids/index.md

@@ -812,10 +812,10 @@ flow state. $P_\mathrm{ref}$ is defined via the `-reference_pressure` flag.
 
 ### Meshing
 
-The flat plate boundary layer example has custom meshing features to better
-resolve the flow. One of those is tilting the top of the domain, allowing for
-it to be a outflow boundary condition. The angle of this tilt is controlled by
-`-platemesh_top_angle`
+The flat plate boundary layer example has custom meshing features to better resolve the flow when using a generated box mesh.
+These meshing features modify the nodal layout of the default, equispaced box mesh and are enabled via `-mesh_transform platemesh`.
+One of those is tilting the top of the domain, allowing for it to be a outflow boundary condition.
+The angle of this tilt is controlled by `-platemesh_top_angle`.
 
 The primary meshing feature is the ability to grade the mesh, providing better
 resolution near the wall. There are two methods to do this; algorithmically, or

examples/fluids/navierstokes.h

@@ -87,13 +87,20 @@ typedef enum {
 } TestType;
 static const char *const TestTypes[] = {"none", "solver", "turb_spanstats", "diff_filter", "TestType", "TESTTYPE_", NULL};
 
-// Test mode type
+// Subgrid-Stress mode type
 typedef enum {
   SGS_MODEL_NONE        = 0,
   SGS_MODEL_DATA_DRIVEN = 1,
 } SGSModelType;
 static const char *const SGSModelTypes[] = {"none", "data_driven", "SGSModelType", "SGS_MODEL_", NULL};
 
+// Mesh transformation type
+typedef enum {
+  MESH_TRANSFORM_NONE      = 0,
+  MESH_TRANSFORM_PLATEMESH = 1,
+} MeshTransformType;
+static const char *const MeshTransformTypes[] = {"none", "platemesh", "MeshTransformType", "MESH_TRANSFORM_", NULL};
+
 static const char *const DifferentialFilterDampingFunctions[] = {
     "none", "van_driest", "mms", "DifferentialFilterDampingFunction", "DIFF_FILTER_DAMP_", NULL};
 
@@ -160,7 +167,8 @@ struct AppCtx_private {
   // Subgrid Stress Model
   SGSModelType sgs_model_type;
   // Differential Filtering
-  PetscBool diff_filter_monitor;
+  PetscBool         diff_filter_monitor;
+  MeshTransformType mesh_transform_type;
 };
 
 // libCEED data struct

examples/fluids/problems/blasius.c

@@ -271,16 +271,18 @@ PetscErrorCode NS_BLASIUS(ProblemData *problem, DM dm, void *ctx, SimpleBC bc) {
   PetscCall(PetscOptionsInt("-n_chebyshev", "Number of Chebyshev terms", NULL, N, &N, NULL));
   PetscCheck(3 <= N && N <= BLASIUS_MAX_N_CHEBYSHEV, comm, PETSC_ERR_ARG_OUTOFRANGE, "-n_chebyshev %" PetscInt_FMT " must be in range [3, %d]", N,
              BLASIUS_MAX_N_CHEBYSHEV);
-  PetscCall(PetscOptionsBoundedInt("-platemesh_Ndelta", "Velocity at boundary layer edge", NULL, mesh_Ndelta, &mesh_Ndelta, NULL, 1));
-  PetscCall(PetscOptionsScalar("-platemesh_refine_height", "Height of boundary layer mesh refinement", NULL, mesh_refine_height, &mesh_refine_height,
-                               NULL));
-  PetscCall(PetscOptionsScalar("-platemesh_growth", "Geometric growth rate of boundary layer mesh", NULL, mesh_growth, &mesh_growth, NULL));
-  PetscCall(
-      PetscOptionsScalar("-platemesh_top_angle", "Geometric top_angle rate of boundary layer mesh", NULL, mesh_top_angle, &mesh_top_angle, NULL));
-  PetscCall(PetscOptionsString("-platemesh_y_node_locs_path",
-                               "Path to file with y node locations. "
-                               "If empty, will use the algorithmic mesh warping.",
-                               NULL, mesh_ynodes_path, mesh_ynodes_path, sizeof(mesh_ynodes_path), NULL));
+  if (user->app_ctx->mesh_transform_type == MESH_TRANSFORM_PLATEMESH) {
+    PetscCall(PetscOptionsBoundedInt("-platemesh_Ndelta", "Velocity at boundary layer edge", NULL, mesh_Ndelta, &mesh_Ndelta, NULL, 1));
+    PetscCall(PetscOptionsScalar("-platemesh_refine_height", "Height of boundary layer mesh refinement", NULL, mesh_refine_height,
+                                 &mesh_refine_height, NULL));
+    PetscCall(PetscOptionsScalar("-platemesh_growth", "Geometric growth rate of boundary layer mesh", NULL, mesh_growth, &mesh_growth, NULL));
+    PetscCall(
+        PetscOptionsScalar("-platemesh_top_angle", "Geometric top_angle rate of boundary layer mesh", NULL, mesh_top_angle, &mesh_top_angle, NULL));
+    PetscCall(PetscOptionsString("-platemesh_y_node_locs_path",
+                                 "Path to file with y node locations. "
+                                 "If empty, will use the algorithmic mesh warping.",
+                                 NULL, mesh_ynodes_path, mesh_ynodes_path, sizeof(mesh_ynodes_path), NULL));
+  }
   PetscCall(PetscOptionsBool("-stg_use", "Use STG inflow boundary condition", NULL, use_stg, &use_stg, NULL));
   PetscOptionsEnd();
 
@@ -295,12 +297,13 @@ PetscErrorCode NS_BLASIUS(ProblemData *problem, DM dm, void *ctx, SimpleBC bc) {
   U_inf *= meter / second;
   delta0 *= meter;
 
-  PetscReal *mesh_ynodes  = NULL;
-  PetscInt   mesh_nynodes = 0;
-  if (strcmp(mesh_ynodes_path, "")) {
-    PetscCall(GetYNodeLocs(comm, mesh_ynodes_path, &mesh_ynodes, &mesh_nynodes));
+  if (user->app_ctx->mesh_transform_type == MESH_TRANSFORM_PLATEMESH) {
+    PetscReal *mesh_ynodes  = NULL;
+    PetscInt   mesh_nynodes = 0;
+    if (strcmp(mesh_ynodes_path, "")) PetscCall(GetYNodeLocs(comm, mesh_ynodes_path, &mesh_ynodes, &mesh_nynodes));
+    PetscCall(ModifyMesh(comm, dm, problem->dim, mesh_growth, mesh_Ndelta, mesh_refine_height, mesh_top_angle, &mesh_ynodes, &mesh_nynodes));
+    PetscCall(PetscFree(mesh_ynodes));
   }
-  PetscCall(ModifyMesh(comm, dm, problem->dim, mesh_growth, mesh_Ndelta, mesh_refine_height, mesh_top_angle, &mesh_ynodes, &mesh_nynodes));
 
   // Some properties depend on parameters from NewtonianIdealGas
   PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfunction_context, CEED_MEM_HOST, &newtonian_ig_ctx));
@@ -335,7 +338,7 @@ PetscErrorCode NS_BLASIUS(ProblemData *problem, DM dm, void *ctx, SimpleBC bc) {
   PetscCallCeed(ceed, CeedQFunctionContextDestroy(&problem->ics.qfunction_context));
   problem->ics.qfunction_context = blasius_context;
   if (use_stg) {
-    PetscCall(SetupSTG(comm, dm, problem, user, weakT, T_inf, P0, mesh_ynodes, mesh_nynodes));
+    PetscCall(SetupSTG(comm, dm, problem, user, weakT, T_inf, P0));
   } else if (diff_filter_mms) {
     PetscCall(DifferentialFilter_MMS_ICSetup(problem));
   } else {
@@ -346,6 +349,5 @@ PetscErrorCode NS_BLASIUS(ProblemData *problem, DM dm, void *ctx, SimpleBC bc) {
     PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(blasius_context, &problem->apply_inflow.qfunction_context));
     PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(blasius_context, &problem->apply_inflow_jacobian.qfunction_context));
   }
-  PetscCall(PetscFree(mesh_ynodes));
   PetscFunctionReturn(PETSC_SUCCESS);
 }

examples/fluids/problems/stg_shur14.c

@@ -165,7 +165,7 @@ static PetscErrorCode ReadSTGRand(const MPI_Comm comm, const char path[PETSC_MAX
  * @param[in,out] stg_ctx         Pointer to STGShur14Context where the data will be loaded into
  */
 PetscErrorCode GetSTGContextData(const MPI_Comm comm, const DM dm, char stg_inflow_path[PETSC_MAX_PATH_LEN], char stg_rand_path[PETSC_MAX_PATH_LEN],
-                                 STGShur14Context *stg_ctx, const CeedScalar ynodes[]) {
+                                 STGShur14Context *stg_ctx) {
   PetscInt nmodes, nprofs;
   PetscFunctionBeginUser;
 
@@ -191,8 +191,7 @@ PetscErrorCode GetSTGContextData(const MPI_Comm comm, const DM dm, char stg_infl
     temp_ctx->offsets.cij       = temp_ctx->offsets.ubar + nprofs * 3;
     temp_ctx->offsets.eps       = temp_ctx->offsets.cij + nprofs * 6;
     temp_ctx->offsets.lt        = temp_ctx->offsets.eps + nprofs;
-    temp_ctx->offsets.ynodes    = temp_ctx->offsets.lt + nprofs;
-    PetscInt total_num_scalars  = temp_ctx->offsets.ynodes + temp_ctx->nynodes;
+    PetscInt total_num_scalars  = temp_ctx->offsets.lt + nprofs;
     temp_ctx->total_bytes       = sizeof(*temp_ctx) + total_num_scalars * sizeof(temp_ctx->data[0]);
     PetscCall(PetscFree(*stg_ctx));
     PetscCall(PetscMalloc(temp_ctx->total_bytes, stg_ctx));
@@ -203,11 +202,6 @@ PetscErrorCode GetSTGContextData(const MPI_Comm comm, const DM dm, char stg_infl
   PetscCall(ReadSTGInflow(comm, stg_inflow_path, *stg_ctx));
   PetscCall(ReadSTGRand(comm, stg_rand_path, *stg_ctx));
 
-  if ((*stg_ctx)->nynodes > 0) {
-    CeedScalar *ynodes_ctx = &(*stg_ctx)->data[(*stg_ctx)->offsets.ynodes];
-    for (PetscInt i = 0; i < (*stg_ctx)->nynodes; i++) ynodes_ctx[i] = ynodes[i];
-  }
-
   {  // -- Calculate kappa
     CeedScalar *kappa     = &(*stg_ctx)->data[(*stg_ctx)->offsets.kappa];
     CeedScalar *wall_dist = &(*stg_ctx)->data[(*stg_ctx)->offsets.wall_dist];
@@ -227,12 +221,12 @@ PetscErrorCode GetSTGContextData(const MPI_Comm comm, const DM dm, char stg_infl
 }
 
 PetscErrorCode SetupSTG(const MPI_Comm comm, const DM dm, ProblemData *problem, User user, const bool prescribe_T, const CeedScalar theta0,
-                        const CeedScalar P0, const CeedScalar ynodes[], const CeedInt nynodes) {
+                        const CeedScalar P0) {
   Ceed                     ceed                                = user->ceed;
   char                     stg_inflow_path[PETSC_MAX_PATH_LEN] = "./STGInflow.dat";
   char                     stg_rand_path[PETSC_MAX_PATH_LEN]   = "./STGRand.dat";
-  PetscBool                mean_only = PETSC_FALSE, use_stgstrong = PETSC_FALSE, use_fluctuating_IC = PETSC_FALSE;
-  CeedScalar               u0 = 0.0, alpha = 1.01;
+  PetscBool                mean_only = PETSC_FALSE, use_stgstrong = PETSC_FALSE, use_fluctuating_IC = PETSC_FALSE, given_stg_dx = PETSC_FALSE;
+  CeedScalar               u0 = 0.0, alpha = 1.01, stg_dx = 1.0e-3;
   CeedQFunctionContext     stg_context;
   NewtonianIdealGasContext newtonian_ig_ctx;
   PetscFunctionBeginUser;
@@ -247,6 +241,7 @@ PetscErrorCode SetupSTG(const MPI_Comm comm, const DM dm, ProblemData *problem,
   PetscCall(PetscOptionsBool("-stg_strong", "Enforce STG inflow strongly", NULL, use_stgstrong, &use_stgstrong, NULL));
   PetscCall(PetscOptionsBool("-stg_fluctuating_IC", "\"Extrude\" the fluctuations through the domain as an initial condition", NULL,
                              use_fluctuating_IC, &use_fluctuating_IC, NULL));
+  PetscCall(PetscOptionsReal("-stg_dx", "Element size in streamwise direction at inflow", NULL, stg_dx, &stg_dx, &given_stg_dx));
   PetscOptionsEnd();
 
   PetscCall(PetscCalloc1(1, &global_stg_ctx));
@@ -258,7 +253,6 @@ PetscErrorCode SetupSTG(const MPI_Comm comm, const DM dm, ProblemData *problem,
   global_stg_ctx->use_fluctuating_IC = use_fluctuating_IC;
   global_stg_ctx->theta0             = theta0;
   global_stg_ctx->P0                 = P0;
-  global_stg_ctx->nynodes            = nynodes;
 
   {
     // Calculate dx assuming constant spacing
@@ -268,15 +262,14 @@ PetscErrorCode SetupSTG(const MPI_Comm comm, const DM dm, ProblemData *problem,
 
     PetscInt nmax = 3, faces[3];
     PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &nmax, NULL));
-    global_stg_ctx->dx = domain_size[0] / faces[0];
-    global_stg_ctx->dz = domain_size[2] / faces[2];
+    global_stg_ctx->dx = given_stg_dx ? stg_dx : domain_size[0] / faces[0];
   }
 
   PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfunction_context, CEED_MEM_HOST, &newtonian_ig_ctx));
   global_stg_ctx->newtonian_ctx = *newtonian_ig_ctx;
   PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfunction_context, &newtonian_ig_ctx));
 
-  PetscCall(GetSTGContextData(comm, dm, stg_inflow_path, stg_rand_path, &global_stg_ctx, ynodes));
+  PetscCall(GetSTGContextData(comm, dm, stg_inflow_path, stg_rand_path, &global_stg_ctx));
 
   PetscCallCeed(ceed, CeedQFunctionContextCreate(user->ceed, &stg_context));
   PetscCallCeed(ceed, CeedQFunctionContextSetData(stg_context, CEED_MEM_HOST, CEED_USE_POINTER, global_stg_ctx->total_bytes, global_stg_ctx));
@@ -306,58 +299,7 @@ PetscErrorCode SetupSTG(const MPI_Comm comm, const DM dm, ProblemData *problem,
   PetscFunctionReturn(PETSC_SUCCESS);
 }
 
-static inline PetscScalar FindDy(const PetscScalar ynodes[], const PetscInt nynodes, const PetscScalar y) {
-  const PetscScalar half_mindy = 0.5 * (ynodes[1] - ynodes[0]);
-  // ^^assuming min(dy) is first element off the wall
-  PetscInt idx = -1;  // Index
-
-  for (PetscInt i = 0; i < nynodes; i++) {
-    if (y < ynodes[i] + half_mindy) {
-      idx = i;
-      break;
-    }
-  }
-  if (idx == 0) return ynodes[1] - ynodes[0];
-  else if (idx == nynodes - 1) return ynodes[nynodes - 2] - ynodes[nynodes - 1];
-  else return 0.5 * (ynodes[idx + 1] - ynodes[idx - 1]);
-}
-
-// Function passed to DMAddBoundary
-// NOTE: Not used in favor of QFunction-based method
-PetscErrorCode StrongSTGbcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[], void *ctx) {
-  PetscFunctionBeginUser;
-
-  const STGShur14Context stg_ctx = (STGShur14Context)ctx;
-  PetscScalar            qn[stg_ctx->nmodes], u[3], ubar[3], cij[6], eps, lt;
-  const bool             mean_only = stg_ctx->mean_only;
-  const PetscScalar      dx        = stg_ctx->dx;
-  const PetscScalar      dz        = stg_ctx->dz;
-  const PetscScalar      mu        = stg_ctx->newtonian_ctx.mu;
-  const PetscScalar      theta0    = stg_ctx->theta0;
-  const PetscScalar      P0        = stg_ctx->P0;
-  const PetscScalar      cv        = stg_ctx->newtonian_ctx.cv;
-  const PetscScalar      cp        = stg_ctx->newtonian_ctx.cp;
-  const PetscScalar      Rd        = cp - cv;
-
-  const CeedScalar rho = P0 / (Rd * theta0);
-  InterpolateProfile(x[1], ubar, cij, &eps, &lt, stg_ctx);
-  if (!mean_only) {
-    const PetscInt     nynodes = stg_ctx->nynodes;
-    const PetscScalar *ynodes  = &stg_ctx->data[stg_ctx->offsets.ynodes];
-    const PetscScalar  h[3]    = {dx, FindDy(ynodes, nynodes, x[1]), dz};
-    CalcSpectrum(x[1], eps, lt, h, mu / rho, qn, stg_ctx);
-    STGShur14_Calc(x, time, ubar, cij, qn, u, stg_ctx);
-  } else {
-    for (CeedInt j = 0; j < 3; j++) u[j] = ubar[j];
-  }
-
-  bcval[0] = rho;
-  bcval[1] = rho * u[0];
-  bcval[2] = rho * u[1];
-  bcval[3] = rho * u[2];
-  PetscFunctionReturn(PETSC_SUCCESS);
-}
-
+// @brief Set STG strongly enforce components using DMAddBoundary
 PetscErrorCode SetupStrongSTG(DM dm, SimpleBC bc, ProblemData *problem, Physics phys) {
   DMLabel label;
   PetscFunctionBeginUser;
@@ -376,10 +318,8 @@ PetscErrorCode SetupStrongSTG(DM dm, SimpleBC bc, ProblemData *problem, Physics
   }
 
   PetscCall(DMGetLabel(dm, "Face Sets", &label));
-  // Set wall BCs
   if (bc->num_inflow > 0) {
-    PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "STG", label, bc->num_inflow, bc->inflows, 0, num_comps, comps, (void (*)(void))StrongSTGbcFunc,
-                            NULL, global_stg_ctx, NULL));
+    PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "STG", label, bc->num_inflow, bc->inflows, 0, num_comps, comps, NULL, NULL, NULL, NULL));
   }
 
   PetscFunctionReturn(PETSC_SUCCESS);

examples/fluids/problems/stg_shur14.h

@@ -12,7 +12,7 @@
 #include "../qfunctions/stg_shur14_type.h"
 
 extern PetscErrorCode SetupSTG(const MPI_Comm comm, const DM dm, ProblemData *problem, User user, const bool prescribe_T, const CeedScalar theta0,
-                               const CeedScalar P0, const CeedScalar ynodes[], const CeedInt nynodes);
+                               const CeedScalar P0);
 
 extern PetscErrorCode SetupStrongSTG(DM dm, SimpleBC bc, ProblemData *problem, Physics phys);
 

examples/fluids/qfunctions/differential_filter.h

@@ -201,7 +201,8 @@ CEED_QFUNCTION(DifferentialFilter_MMS_RHS)(void *ctx, CeedInt Q, const CeedScala
 // It can be run via:
 // ./navierstokes -options_file tests-output/blasius_test.yaml -diff_filter_enable -diff_filter_view cgns:filtered_solution.cgns -ts_max_steps 0
 // -diff_filter_mms -diff_filter_kernel_scaling 1 -diff_filter_wall_damping_function mms -dm_plex_box_upper 0.5,0.5,0.5 -dm_plex_box_faces 75,75,1
-// -platemesh_y_node_locs_path tests-output/diff_filter_mms_y_spacing.dat -platemesh_top_angle 0 -diff_filter_grid_based_width
+// -platemesh_modify_mesh true -platemesh_y_node_locs_path tests-output/diff_filter_mms_y_spacing.dat -platemesh_top_angle 0
+// -diff_filter_grid_based_width
 CEED_QFUNCTION(DifferentialFilter_MMS_IC)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
   const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
   CeedScalar(*q0)[CEED_Q_VLA]      = (CeedScalar(*)[CEED_Q_VLA])out[0];

examples/fluids/qfunctions/stg_shur14.h

@@ -90,13 +90,13 @@ CEED_QFUNCTION_HELPER void InterpolateProfile(const CeedScalar wall_dist, CeedSc
  *
  * Calculates q_n at a given distance to the wall
  *
- * @param[in]  kappa  nth wavenumber
- * @param[in]  dkappa Difference between wavenumbers
- * @param[in]  keta   Dissipation wavenumber
- * @param[in]  kcut   Mesh-induced cutoff wavenumber
- * @param[in]  ke     Energy-containing wavenumber
- * @param[in]  Ektot  Total turbulent kinetic energy of spectrum
- * @returns    qn     Spectrum coefficient
+ * @param[in]  kappa     nth wavenumber
+ * @param[in]  dkappa    Difference between wavenumbers
+ * @param[in]  keta      Dissipation wavenumber
+ * @param[in]  kcut      Mesh-induced cutoff wavenumber
+ * @param[in]  ke        Energy-containing wavenumber
+ * @param[in]  Ektot_inv Inverse of total turbulent kinetic energy of spectrum
+ * @returns    qn        Spectrum coefficient
  */
 CEED_QFUNCTION_HELPER CeedScalar Calc_qn(const CeedScalar kappa, const CeedScalar dkappa, const CeedScalar keta, const CeedScalar kcut,
                                          const CeedScalar ke, const CeedScalar Ektot_inv) {
@@ -523,7 +523,7 @@ CEED_QFUNCTION(STGShur14_Inflow_StrongQF)(void *ctx, CeedInt Q, const CeedScalar
   const CeedScalar(*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
   const CeedScalar(*coords)[CEED_Q_VLA]     = (const CeedScalar(*)[CEED_Q_VLA])in[1];
   const CeedScalar(*scale)                  = (const CeedScalar(*))in[2];
-  const CeedScalar(*stg_data)[CEED_Q_VLA]   = (const CeedScalar(*)[CEED_Q_VLA])in[3];
+  const CeedScalar(*inv_Ektotal)            = (const CeedScalar(*))in[3];
 
   CeedScalar(*bcval)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
 
@@ -551,7 +551,7 @@ CEED_QFUNCTION(STGShur14_Inflow_StrongQF)(void *ctx, CeedInt Q, const CeedScalar
     InterpolateProfile(coords[1][i], ubar, cij, &eps, &lt, stg_ctx);
     if (!mean_only) {
       if (1) {
-        STGShur14_Calc_PrecompEktot(x, time, ubar, cij, stg_data[0][i], h, x[1], eps, lt, nu, u, stg_ctx);
+        STGShur14_Calc_PrecompEktot(x, time, ubar, cij, inv_Ektotal[i], h, x[1], eps, lt, nu, u, stg_ctx);
       } else {  // Original way
         CeedScalar qn[STG_NMODES_MAX];
         CalcSpectrum(coords[1][i], eps, lt, h, nu, qn, stg_ctx);