xref: /honee/problems/stg_shur14.c (revision ea091b8e43283667a8e2f5f647273e58272b0e63)
1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3 
4 /// @file
5 /// Implementation of the Synthetic Turbulence Generation (STG) algorithm
6 /// presented in Shur et al. 2014
7 
8 #include "stg_shur14.h"
9 
10 #include <ceed.h>
11 #include <math.h>
12 #include <petscdm.h>
13 #include <stdlib.h>
14 
15 #include <navierstokes.h>
16 #include "../qfunctions/stg_shur14.h"
17 
18 StgShur14Context global_stg_ctx;
19 
20 /*
21  * @brief Perform Cholesky decomposition on array of symmetric 3x3 matrices
22  *
23  * This assumes the input matrices are in order [11,22,33,12,13,23].
24  * This format is also used for the output.
25  *
26  * @param[in]  comm   MPI_Comm
27  * @param[in]  nprofs Number of matrices in Rij
28  * @param[in]  Rij    Array of the symmetric matrices [6,nprofs]
29  * @param[out] Cij    Array of the Cholesky Decomposition matrices, [6,nprofs]
30  */
31 PetscErrorCode CalcCholeskyDecomp(MPI_Comm comm, PetscInt nprofs, const CeedScalar Rij[6][nprofs], CeedScalar Cij[6][nprofs]) {
32   PetscFunctionBeginUser;
33   for (PetscInt i = 0; i < nprofs; i++) {
34     Cij[0][i] = sqrt(Rij[0][i]);
35     Cij[3][i] = Rij[3][i] / Cij[0][i];
36     Cij[1][i] = sqrt(Rij[1][i] - Square(Cij[3][i]));
37     Cij[4][i] = Rij[4][i] / Cij[0][i];
38     Cij[5][i] = (Rij[5][i] - Cij[3][i] * Cij[4][i]) / Cij[1][i];
39     Cij[2][i] = sqrt(Rij[2][i] - Square(Cij[4][i]) - Square(Cij[5][i]));
40 
41     PetscCheck(!isnan(Cij[0][i]) && !isnan(Cij[1][i]) && !isnan(Cij[2][i]), comm, PETSC_ERR_FP,
42                "Cholesky decomposition failed at profile point %" PetscInt_FMT ". Either STGInflow has non-SPD matrix or contains nan.", i + 1);
43   }
44   PetscFunctionReturn(PETSC_SUCCESS);
45 }
46 
47 /*
48  * @brief Read the STGInflow file and load the contents into stg_ctx
49  *
50  * Assumes that the first line of the file has the number of rows and columns as the only two entries, separated by a single space.
51  * Assumes there are 14 columns in the file.
52  *
53  * Function calculates the Cholesky decomposition from the Reynolds stress profile found in the file.
54  *
55  * @param[in]     comm    MPI_Comm for the program
56  * @param[in]     path    Path to the STGInflow.dat file
57  * @param[in,out] stg_ctx STGShur14Context where the data will be loaded into
58  */
59 static PetscErrorCode ReadStgInflow(const MPI_Comm comm, const char path[PETSC_MAX_PATH_LEN], StgShur14Context stg_ctx) {
60   PetscInt       dims[2];
61   int            ndims;
62   FILE          *fp;
63   const PetscInt char_array_len = 512;
64   char           line[char_array_len];
65   char         **array;
66 
67   PetscFunctionBeginUser;
68   PetscCall(PhastaDatFileOpen(comm, path, char_array_len, dims, &fp));
69 
70   CeedScalar  rij[6][stg_ctx->nprofs];
71   CeedScalar *wall_dist              = &stg_ctx->data[stg_ctx->offsets.wall_dist];
72   CeedScalar *eps                    = &stg_ctx->data[stg_ctx->offsets.eps];
73   CeedScalar *lt                     = &stg_ctx->data[stg_ctx->offsets.lt];
74   CeedScalar(*ubar)[stg_ctx->nprofs] = (CeedScalar(*)[stg_ctx->nprofs]) & stg_ctx->data[stg_ctx->offsets.ubar];
75 
76   for (PetscInt i = 0; i < stg_ctx->nprofs; i++) {
77     PetscCall(PetscSynchronizedFGets(comm, fp, char_array_len, line));
78     PetscCall(PetscStrToArray(line, ' ', &ndims, &array));
79     PetscCheck(ndims == dims[1], comm, PETSC_ERR_FILE_UNEXPECTED,
80                "Line %" PetscInt_FMT " of %s does not have correct number of columns (%d instead of %" PetscInt_FMT ")", i, path, ndims, dims[1]);
81 
82     wall_dist[i] = (CeedScalar)atof(array[0]);
83     ubar[0][i]   = (CeedScalar)atof(array[1]);
84     ubar[1][i]   = (CeedScalar)atof(array[2]);
85     ubar[2][i]   = (CeedScalar)atof(array[3]);
86     rij[0][i]    = (CeedScalar)atof(array[4]);
87     rij[1][i]    = (CeedScalar)atof(array[5]);
88     rij[2][i]    = (CeedScalar)atof(array[6]);
89     rij[3][i]    = (CeedScalar)atof(array[7]);
90     rij[4][i]    = (CeedScalar)atof(array[8]);
91     rij[5][i]    = (CeedScalar)atof(array[9]);
92     lt[i]        = (CeedScalar)atof(array[12]);
93     eps[i]       = (CeedScalar)atof(array[13]);
94 
95     PetscCheck(wall_dist[i] >= 0, comm, PETSC_ERR_FILE_UNEXPECTED, "Distance to wall in %s cannot be negative", path);
96     PetscCheck(lt[i] >= 0, comm, PETSC_ERR_FILE_UNEXPECTED, "Turbulent length scale in %s cannot be negative", path);
97     PetscCheck(eps[i] >= 0, comm, PETSC_ERR_FILE_UNEXPECTED, "Turbulent dissipation in %s cannot be negative", path);
98     PetscCall(PetscStrToArrayDestroy(ndims, array));
99   }
100   CeedScalar(*cij)[stg_ctx->nprofs] = (CeedScalar(*)[stg_ctx->nprofs]) & stg_ctx->data[stg_ctx->offsets.cij];
101   PetscCall(CalcCholeskyDecomp(comm, stg_ctx->nprofs, rij, cij));
102   PetscCall(PetscFClose(comm, fp));
103   PetscFunctionReturn(PETSC_SUCCESS);
104 }
105 
106 /*
107  * @brief Read the STGRand file and load the contents into stg_ctx
108  *
109  * Assumes that the first line of the file has the number of rows and columns as the only two entries, separated by a single space.
110  * Assumes there are 7 columns in the file.
111  *
112  * @param[in]     comm    MPI_Comm for the program
113  * @param[in]     path    Path to the STGRand.dat file
114  * @param[in,out] stg_ctx STGShur14Context where the data will be loaded into
115  */
116 static PetscErrorCode ReadStgRand(const MPI_Comm comm, const char path[PETSC_MAX_PATH_LEN], StgShur14Context stg_ctx) {
117   PetscInt       dims[2];
118   int            ndims;
119   FILE          *fp;
120   const PetscInt char_array_len = 512;
121   char           line[char_array_len];
122   char         **array;
123 
124   PetscFunctionBeginUser;
125   PetscCall(PhastaDatFileOpen(comm, path, char_array_len, dims, &fp));
126 
127   CeedScalar *phi                     = &stg_ctx->data[stg_ctx->offsets.phi];
128   CeedScalar(*d)[stg_ctx->nmodes]     = (CeedScalar(*)[stg_ctx->nmodes]) & stg_ctx->data[stg_ctx->offsets.d];
129   CeedScalar(*sigma)[stg_ctx->nmodes] = (CeedScalar(*)[stg_ctx->nmodes]) & stg_ctx->data[stg_ctx->offsets.sigma];
130 
131   for (PetscInt i = 0; i < stg_ctx->nmodes; i++) {
132     PetscCall(PetscSynchronizedFGets(comm, fp, char_array_len, line));
133     PetscCall(PetscStrToArray(line, ' ', &ndims, &array));
134     PetscCheck(ndims == dims[1], comm, PETSC_ERR_FILE_UNEXPECTED,
135                "Line %" PetscInt_FMT " of %s does not have correct number of columns (%d instead of %" PetscInt_FMT ")", i, path, ndims, dims[1]);
136 
137     d[0][i]     = (CeedScalar)atof(array[0]);
138     d[1][i]     = (CeedScalar)atof(array[1]);
139     d[2][i]     = (CeedScalar)atof(array[2]);
140     phi[i]      = (CeedScalar)atof(array[3]);
141     sigma[0][i] = (CeedScalar)atof(array[4]);
142     sigma[1][i] = (CeedScalar)atof(array[5]);
143     sigma[2][i] = (CeedScalar)atof(array[6]);
144     PetscCall(PetscStrToArrayDestroy(ndims, array));
145   }
146   PetscCall(PetscFClose(comm, fp));
147   PetscFunctionReturn(PETSC_SUCCESS);
148 }
149 
150 /*
151  * @brief Read STG data from input paths and put in STGShur14Context
152  *
153  * Reads data from input paths and puts them into a STGShur14Context object.
154  * Data stored initially in `*stg_ctx` will be copied over to the new STGShur14Context instance.
155  *
156  * @param[in]     comm            MPI_Comm for the program
157  * @param[in]     dm              DM for the program
158  * @param[in]     stg_inflow_path Path to STGInflow.dat file
159  * @param[in]     stg_rand_path   Path to STGRand.dat file
160  * @param[in,out] stg_ctx         Pointer to STGShur14Context where the data will be loaded into
161  */
162 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],
163                                  StgShur14Context *stg_ctx) {
164   PetscInt nmodes, nprofs;
165 
166   PetscFunctionBeginUser;
167   PetscCall(PhastaDatFileGetNRows(comm, stg_rand_path, &nmodes));
168   PetscCall(PhastaDatFileGetNRows(comm, stg_inflow_path, &nprofs));
169   PetscCheck(nmodes < STG_NMODES_MAX, comm, PETSC_ERR_SUP,
170              "Number of wavemodes in %s (%" PetscInt_FMT ") exceeds STG_NMODES_MAX (%d). Change size of STG_NMODES_MAX and recompile", stg_rand_path,
171              nmodes, STG_NMODES_MAX);
172 
173   {
174     StgShur14Context temp_ctx;
175     PetscCall(PetscCalloc1(1, &temp_ctx));
176     *temp_ctx                   = **stg_ctx;
177     temp_ctx->nmodes            = nmodes;
178     temp_ctx->nprofs            = nprofs;
179     temp_ctx->offsets.sigma     = 0;
180     temp_ctx->offsets.d         = nmodes * 3;
181     temp_ctx->offsets.phi       = temp_ctx->offsets.d + nmodes * 3;
182     temp_ctx->offsets.kappa     = temp_ctx->offsets.phi + nmodes;
183     temp_ctx->offsets.wall_dist = temp_ctx->offsets.kappa + nmodes;
184     temp_ctx->offsets.ubar      = temp_ctx->offsets.wall_dist + nprofs;
185     temp_ctx->offsets.cij       = temp_ctx->offsets.ubar + nprofs * 3;
186     temp_ctx->offsets.eps       = temp_ctx->offsets.cij + nprofs * 6;
187     temp_ctx->offsets.lt        = temp_ctx->offsets.eps + nprofs;
188     PetscInt total_num_scalars  = temp_ctx->offsets.lt + nprofs;
189     temp_ctx->total_bytes       = sizeof(*temp_ctx) + total_num_scalars * sizeof(temp_ctx->data[0]);
190     PetscCall(PetscFree(*stg_ctx));
191     PetscCall(PetscMalloc(temp_ctx->total_bytes, stg_ctx));
192     **stg_ctx = *temp_ctx;
193     PetscCall(PetscFree(temp_ctx));
194   }
195 
196   PetscCall(ReadStgInflow(comm, stg_inflow_path, *stg_ctx));
197   PetscCall(ReadStgRand(comm, stg_rand_path, *stg_ctx));
198 
199   {  // -- Calculate kappa
200     CeedScalar *kappa     = &(*stg_ctx)->data[(*stg_ctx)->offsets.kappa];
201     CeedScalar *wall_dist = &(*stg_ctx)->data[(*stg_ctx)->offsets.wall_dist];
202     CeedScalar *lt        = &(*stg_ctx)->data[(*stg_ctx)->offsets.lt];
203     CeedScalar  le, le_max = 0;
204 
205     CeedPragmaSIMD for (PetscInt i = 0; i < (*stg_ctx)->nprofs; i++) {
206       le = PetscMin(2 * wall_dist[i], 3 * lt[i]);
207       if (le_max < le) le_max = le;
208     }
209     CeedScalar kmin = M_PI / le_max;
210 
211     CeedPragmaSIMD for (PetscInt i = 0; i < (*stg_ctx)->nmodes; i++) { kappa[i] = kmin * pow((*stg_ctx)->alpha, i); }
212   }
213   PetscFunctionReturn(PETSC_SUCCESS);
214 }
215 
216 static PetscErrorCode STGWeakInflowBCSetup_CreateIFunctionQF(BCDefinition bc_def, CeedQFunction *qf) {
217   HoneeBCStruct honee_bc;
218 
219   PetscFunctionBeginUser;
220   PetscCall(BCDefinitionGetContext(bc_def, &honee_bc));
221   PetscCall(HoneeBCCreateIFunctionQF(bc_def, StgShur14Inflow, StgShur14Inflow_loc, honee_bc->qfctx, qf));
222   PetscFunctionReturn(PETSC_SUCCESS);
223 }
224 
225 static PetscErrorCode STGWeakInflowBCSetup_CreateIJacobianQF(BCDefinition bc_def, CeedQFunction *qf) {
226   HoneeBCStruct honee_bc;
227 
228   PetscFunctionBeginUser;
229   PetscCall(BCDefinitionGetContext(bc_def, &honee_bc));
230   PetscCall(HoneeBCCreateIJacobianQF(bc_def, StgShur14Inflow_Jacobian, StgShur14Inflow_Jacobian_loc, honee_bc->qfctx, qf));
231   PetscFunctionReturn(PETSC_SUCCESS);
232 }
233 
234 PetscErrorCode SetupStg(const MPI_Comm comm, const DM dm, ProblemData problem, Honee honee, const bool prescribe_T, const CeedScalar theta0,
235                         const CeedScalar P0) {
236   Ceed                     ceed                                = honee->ceed;
237   char                     stg_inflow_path[PETSC_MAX_PATH_LEN] = "./STGInflow.dat";
238   char                     stg_rand_path[PETSC_MAX_PATH_LEN]   = "./STGRand.dat";
239   PetscBool                mean_only = PETSC_FALSE, use_stgstrong = PETSC_FALSE, use_fluctuating_IC = PETSC_FALSE, given_stg_dx = PETSC_FALSE;
240   CeedScalar               u0 = 0.0, alpha = 1.01, stg_dx = -1, stg_h_scale_factor = 1 / honee->app_ctx->degree;
241   CeedQFunctionContext     stg_qfctx;
242   NewtonianIdealGasContext newtonian_ig_ctx;
243 
244   PetscFunctionBeginUser;
245   PetscOptionsBegin(comm, NULL, "STG Boundary Condition Options", NULL);
246   PetscCall(PetscOptionsString("-stg_inflow_path", "Path to STGInflow.dat", NULL, stg_inflow_path, stg_inflow_path, sizeof(stg_inflow_path), NULL));
247   PetscCall(PetscOptionsString("-stg_rand_path", "Path to STGInflow.dat", NULL, stg_rand_path, stg_rand_path, sizeof(stg_rand_path), NULL));
248   PetscCall(PetscOptionsReal("-stg_alpha", "Growth rate of the wavemodes", NULL, alpha, &alpha, NULL));
249   PetscCall(PetscOptionsReal("-stg_u0", "Advective velocity for the fluctuations", NULL, u0, &u0, NULL));
250   PetscCall(PetscOptionsBool("-stg_mean_only", "Only apply mean profile", NULL, mean_only, &mean_only, NULL));
251   PetscCall(PetscOptionsBool("-stg_strong", "Enforce STG inflow strongly", NULL, use_stgstrong, &use_stgstrong, NULL));
252   PetscCall(PetscOptionsBool("-stg_fluctuating_IC", "\"Extrude\" the fluctuations through the domain as an initial condition", NULL,
253                              use_fluctuating_IC, &use_fluctuating_IC, NULL));
254   PetscCall(PetscOptionsReal("-stg_dx", "Element length in x direction at inflow", NULL, stg_dx, &stg_dx, &given_stg_dx));
255   if (given_stg_dx && use_stgstrong) PetscCall(PetscPrintf(comm, "WARNING: -stg_dx is ignored for -stg_strong\n"));
256   PetscCall(PetscOptionsReal("-stg_h_scale_factor", "Scale element size for cutoff frequency calculation", NULL, stg_h_scale_factor,
257                              &stg_h_scale_factor, NULL));
258   PetscCall(PetscOptionsDeprecated("-stg_dyScale", NULL, "libCEED 0.12.0", "Use -stg_h_scale_factor to scale all the element dimensions"));
259   PetscCall(PetscOptionsDeprecated("-stg_dz", NULL, "libCEED 0.12.0", NULL));
260   PetscOptionsEnd();
261 
262   PetscCall(PetscCalloc1(1, &global_stg_ctx));
263   global_stg_ctx->alpha              = alpha;
264   global_stg_ctx->u0                 = u0;
265   global_stg_ctx->is_implicit        = honee->phys->implicit;
266   global_stg_ctx->prescribe_T        = prescribe_T;
267   global_stg_ctx->mean_only          = mean_only;
268   global_stg_ctx->use_fluctuating_IC = use_fluctuating_IC;
269   global_stg_ctx->theta0             = theta0;
270   global_stg_ctx->P0                 = P0;
271   global_stg_ctx->h_scale_factor     = stg_h_scale_factor;
272 
273   if (!use_stgstrong) {  // Calculate dx assuming constant spacing
274     PetscReal domain_min[3], domain_max[3], domain_size[3];
275     PetscCall(DMGetBoundingBox(dm, domain_min, domain_max));
276     for (PetscInt i = 0; i < 3; i++) domain_size[i] = domain_max[i] - domain_min[i];
277 
278     PetscInt nmax = 3, faces[3];
279     PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &nmax, NULL));
280     global_stg_ctx->dx = given_stg_dx ? stg_dx : domain_size[0] / faces[0];
281     PetscCheck((global_stg_ctx->dx > 0) && PetscIsNormalReal((PetscReal)global_stg_ctx->dx), comm, PETSC_ERR_LIB,
282                "STG dx must be positive normal number, got %g", global_stg_ctx->dx);
283   }
284 
285   PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfctx, CEED_MEM_HOST, &newtonian_ig_ctx));
286   global_stg_ctx->newtonian_ctx = *newtonian_ig_ctx;
287   PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfctx, &newtonian_ig_ctx));
288 
289   PetscCall(GetStgContextData(comm, dm, stg_inflow_path, stg_rand_path, &global_stg_ctx));
290 
291   PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &stg_qfctx));
292   PetscCallCeed(ceed, CeedQFunctionContextSetData(stg_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, global_stg_ctx->total_bytes, global_stg_ctx));
293   PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(stg_qfctx, CEED_MEM_HOST, FreeContextPetsc));
294   PetscCallCeed(ceed, CeedQFunctionContextRegisterDouble(stg_qfctx, "solution time", offsetof(struct STGShur14Context_, time), 1,
295                                                          "Physical time of the solution"));
296 
297   PetscCallCeed(ceed, CeedQFunctionContextDestroy(&problem->ics.qfctx));
298   problem->ics.qf_func_ptr = ICsStg;
299   problem->ics.qf_loc      = ICsStg_loc;
300   problem->ics.qfctx       = stg_qfctx;
301 
302   if (use_stgstrong) {
303     // Use default boundary integral QF (BoundaryIntegral) in newtonian.h
304     problem->use_strong_bc_ceed = PETSC_TRUE;
305     problem->set_bc_from_ics    = PETSC_FALSE;
306   } else {
307     for (PetscCount b = 0; b < problem->num_bc_defs; b++) {
308       BCDefinition bc_def = problem->bc_defs[b];
309       const char  *name;
310 
311       PetscCall(BCDefinitionGetInfo(bc_def, &name, NULL, NULL));
312       if (!strcmp(name, "inflow")) {
313         HoneeBCStruct honee_bc;
314 
315         PetscCall(PetscNew(&honee_bc));
316         PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(stg_qfctx, &honee_bc->qfctx));
317         honee_bc->honee             = honee;
318         honee_bc->jac_data_size_sur = honee->phys->implicit ? 11 : 0;
319         PetscCall(BCDefinitionSetContext(bc_def, HoneeBCDestroy, honee_bc));
320 
321         PetscCall(BCDefinitionSetIFunction(bc_def, STGWeakInflowBCSetup_CreateIFunctionQF, HoneeBCAddIFunctionOp));
322         PetscCall(BCDefinitionSetIJacobian(bc_def, STGWeakInflowBCSetup_CreateIJacobianQF, HoneeBCAddIJacobianOp));
323       }
324     }
325     problem->set_bc_from_ics = PETSC_TRUE;
326   }
327   PetscFunctionReturn(PETSC_SUCCESS);
328 }
329 
330 // @brief Set STG strongly enforce components using DMAddBoundary
331 PetscErrorCode SetupStrongStg(DM dm, ProblemData problem, Physics phys) {
332   DMLabel  label;
333   PetscInt comps[5], num_comps = 4;
334 
335   PetscFunctionBeginUser;
336   switch (phys->state_var) {
337     case STATEVAR_CONSERVATIVE:
338       // {0,1,2,3} for rho, rho*u, rho*v, rho*w
339       for (int i = 0; i < 4; i++) comps[i] = i;
340       break;
341 
342     case STATEVAR_PRIMITIVE:
343       // {1,2,3,4} for u, v, w, T
344       for (int i = 0; i < 4; i++) comps[i] = i + 1;
345       break;
346 
347     case STATEVAR_ENTROPY:
348       // {1,2,3,4}
349       for (int i = 0; i < 4; i++) comps[i] = i + 1;
350       break;
351   }
352 
353   PetscCall(DMGetLabel(dm, "Face Sets", &label));
354   for (PetscCount b = 0; b < problem->num_bc_defs; b++) {
355     BCDefinition bc_def = problem->bc_defs[b];
356     const char  *name;
357 
358     PetscCall(BCDefinitionGetInfo(bc_def, &name, NULL, NULL));
359     if (!strcmp(name, "inflow")) PetscCall(BCDefinitionSetEssential(bc_def, num_comps, comps));
360   }
361   PetscFunctionReturn(PETSC_SUCCESS);
362 }
363 
364 PetscErrorCode SetupStrongStg_QF(Ceed ceed, ProblemData problem, CeedInt num_comp_x, CeedInt num_comp_q, CeedInt stg_data_size, CeedInt dXdx_size,
365                                  CeedQFunction *qf_strongbc) {
366   PetscFunctionBeginUser;
367   PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, StgShur14InflowStrongQF, StgShur14InflowStrongQF_loc, qf_strongbc));
368   PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_strongbc, "dXdx", dXdx_size, CEED_EVAL_NONE));
369   PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_strongbc, "x", num_comp_x, CEED_EVAL_NONE));
370   PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_strongbc, "scale", 1, CEED_EVAL_NONE));
371   PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_strongbc, "stg data", stg_data_size, CEED_EVAL_NONE));
372   PetscCallCeed(ceed, CeedQFunctionAddOutput(*qf_strongbc, "q", num_comp_q, CEED_EVAL_NONE));
373 
374   PetscCallCeed(ceed, CeedQFunctionSetContext(*qf_strongbc, problem->ics.qfctx));
375   PetscFunctionReturn(PETSC_SUCCESS);
376 }
377 
378 PetscErrorCode SetupStrongStg_PreProcessing(Ceed ceed, ProblemData problem, CeedInt num_comp_x, CeedInt stg_data_size, CeedInt dXdx_size,
379                                             CeedQFunction *qf_strongbc) {
380   PetscFunctionBeginUser;
381   PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, StgShur14Preprocess, StgShur14Preprocess_loc, qf_strongbc));
382   PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_strongbc, "dXdx", dXdx_size, CEED_EVAL_NONE));
383   PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_strongbc, "x", num_comp_x, CEED_EVAL_NONE));
384   PetscCallCeed(ceed, CeedQFunctionAddOutput(*qf_strongbc, "stg data", stg_data_size, CEED_EVAL_NONE));
385 
386   PetscCallCeed(ceed, CeedQFunctionSetContext(*qf_strongbc, problem->ics.qfctx));
387   PetscFunctionReturn(PETSC_SUCCESS);
388 }
389