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 */
CalcCholeskyDecomp(MPI_Comm comm,PetscInt nprofs,const CeedScalar Rij[6][nprofs],CeedScalar Cij[6][nprofs])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 */
ReadStgInflow(const MPI_Comm comm,const char path[PETSC_MAX_PATH_LEN],StgShur14Context stg_ctx)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 */
ReadStgRand(const MPI_Comm comm,const char path[PETSC_MAX_PATH_LEN],StgShur14Context stg_ctx)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 */
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)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 = 0, nprofs;
165
166 PetscFunctionBeginUser;
167 PetscCall(PhastaDatFileGetNRows(comm, stg_inflow_path, &nprofs));
168 const PetscBool need_rand = (!(*stg_ctx)->mean_only) || (*stg_ctx)->use_fluctuating_IC;
169 if (need_rand) {
170 PetscCall(PhastaDatFileGetNRows(comm, stg_rand_path, &nmodes));
171 PetscCheck(nmodes < STG_NMODES_MAX, comm, PETSC_ERR_SUP,
172 "Number of wavemodes in %s (%" PetscInt_FMT ") exceeds STG_NMODES_MAX (%d). Change size of STG_NMODES_MAX and recompile",
173 stg_rand_path, nmodes, STG_NMODES_MAX);
174 }
175
176 {
177 StgShur14Context temp_ctx;
178 PetscCall(PetscNew(&temp_ctx));
179 *temp_ctx = **stg_ctx;
180 temp_ctx->nmodes = nmodes;
181 temp_ctx->nprofs = nprofs;
182 // nmode = 0 if random numbers are not read from file, therefore offsets will be correctly handled
183 temp_ctx->offsets.sigma = 0;
184 temp_ctx->offsets.d = nmodes * 3;
185 temp_ctx->offsets.phi = temp_ctx->offsets.d + nmodes * 3;
186 temp_ctx->offsets.kappa = temp_ctx->offsets.phi + nmodes;
187 temp_ctx->offsets.wall_dist = temp_ctx->offsets.kappa + nmodes;
188 temp_ctx->offsets.ubar = temp_ctx->offsets.wall_dist + nprofs;
189 temp_ctx->offsets.cij = temp_ctx->offsets.ubar + nprofs * 3;
190 temp_ctx->offsets.eps = temp_ctx->offsets.cij + nprofs * 6;
191 temp_ctx->offsets.lt = temp_ctx->offsets.eps + nprofs;
192 PetscInt total_num_scalars = temp_ctx->offsets.lt + nprofs;
193 temp_ctx->total_bytes = sizeof(*temp_ctx) + total_num_scalars * sizeof(temp_ctx->data[0]);
194 PetscCall(PetscFree(*stg_ctx));
195 PetscCall(PetscMalloc(temp_ctx->total_bytes, stg_ctx));
196 **stg_ctx = *temp_ctx;
197 PetscCall(PetscFree(temp_ctx));
198 }
199
200 PetscCall(ReadStgInflow(comm, stg_inflow_path, *stg_ctx));
201 if (need_rand) {
202 PetscCall(ReadStgRand(comm, stg_rand_path, *stg_ctx));
203 CeedScalar *kappa = &(*stg_ctx)->data[(*stg_ctx)->offsets.kappa];
204 CeedScalar *wall_dist = &(*stg_ctx)->data[(*stg_ctx)->offsets.wall_dist];
205 CeedScalar *lt = &(*stg_ctx)->data[(*stg_ctx)->offsets.lt];
206 CeedScalar le, le_max = 0;
207
208 CeedPragmaSIMD for (PetscInt i = 0; i < (*stg_ctx)->nprofs; i++) {
209 le = PetscMin(2 * wall_dist[i], 3 * lt[i]);
210 if (le_max < le) le_max = le;
211 }
212 CeedScalar kmin = M_PI / le_max;
213
214 CeedPragmaSIMD for (PetscInt i = 0; i < (*stg_ctx)->nmodes; i++) kappa[i] = kmin * pow((*stg_ctx)->alpha, i);
215 }
216 PetscFunctionReturn(PETSC_SUCCESS);
217 }
218
STGWeakInflowBCSetup_CreateIFunctionQF(BCDefinition bc_def,CeedQFunction * qf)219 static PetscErrorCode STGWeakInflowBCSetup_CreateIFunctionQF(BCDefinition bc_def, CeedQFunction *qf) {
220 HoneeBCStruct honee_bc;
221
222 PetscFunctionBeginUser;
223 PetscCall(BCDefinitionGetContext(bc_def, &honee_bc));
224 PetscCall(HoneeBCCreateIFunctionQF(bc_def, StgShur14Inflow, StgShur14Inflow_loc, honee_bc->qfctx, qf));
225 PetscFunctionReturn(PETSC_SUCCESS);
226 }
227
STGWeakInflowBCSetup_CreateIJacobianQF(BCDefinition bc_def,CeedQFunction * qf)228 static PetscErrorCode STGWeakInflowBCSetup_CreateIJacobianQF(BCDefinition bc_def, CeedQFunction *qf) {
229 HoneeBCStruct honee_bc;
230
231 PetscFunctionBeginUser;
232 PetscCall(BCDefinitionGetContext(bc_def, &honee_bc));
233 PetscCall(HoneeBCCreateIJacobianQF(bc_def, StgShur14Inflow_Jacobian, StgShur14Inflow_Jacobian_loc, honee_bc->qfctx, qf));
234 PetscFunctionReturn(PETSC_SUCCESS);
235 }
236
SetupStg(const MPI_Comm comm,const DM dm,ProblemData problem,Honee honee,const bool prescribe_T,const CeedScalar theta0,const CeedScalar P0)237 PetscErrorCode SetupStg(const MPI_Comm comm, const DM dm, ProblemData problem, Honee honee, const bool prescribe_T, const CeedScalar theta0,
238 const CeedScalar P0) {
239 Ceed ceed = honee->ceed;
240 char stg_inflow_path[PETSC_MAX_PATH_LEN] = "./STGInflow.dat";
241 char stg_rand_path[PETSC_MAX_PATH_LEN] = "./STGRand.dat";
242 PetscBool mean_only = PETSC_FALSE, use_stgstrong = PETSC_FALSE, use_fluctuating_IC = PETSC_FALSE, given_stg_dx = PETSC_FALSE;
243 CeedScalar u0 = 0.0, alpha = 1.01, stg_dx = -1, stg_h_scale_factor = 1 / honee->app_ctx->degree;
244 CeedQFunctionContext stg_qfctx;
245 NewtonianIdealGasContext newtonian_ig_ctx;
246
247 PetscFunctionBeginUser;
248 PetscOptionsBegin(comm, NULL, "STG Boundary Condition Options", NULL);
249 PetscCall(PetscOptionsString("-stg_inflow_path", "Path to STGInflow.dat", NULL, stg_inflow_path, stg_inflow_path, sizeof(stg_inflow_path), NULL));
250 PetscCall(PetscOptionsString("-stg_rand_path", "Path to STGRand.dat", NULL, stg_rand_path, stg_rand_path, sizeof(stg_rand_path), NULL));
251 PetscCall(PetscOptionsReal("-stg_alpha", "Growth rate of the wavemodes", NULL, alpha, &alpha, NULL));
252 PetscCall(PetscOptionsReal("-stg_u0", "Advective velocity for the fluctuations", NULL, u0, &u0, NULL));
253 PetscCall(PetscOptionsBool("-stg_mean_only", "Only apply mean profile", NULL, mean_only, &mean_only, NULL));
254 PetscCall(PetscOptionsBool("-stg_strong", "Enforce STG inflow strongly", NULL, use_stgstrong, &use_stgstrong, NULL));
255 PetscCall(PetscOptionsBool("-stg_fluctuating_IC", "\"Extrude\" the fluctuations through the domain as an initial condition", NULL,
256 use_fluctuating_IC, &use_fluctuating_IC, NULL));
257 PetscCall(PetscOptionsReal("-stg_dx", "Element length in x direction at inflow", NULL, stg_dx, &stg_dx, &given_stg_dx));
258 if (given_stg_dx && use_stgstrong) PetscCall(PetscPrintf(comm, "WARNING: -stg_dx is ignored for -stg_strong\n"));
259 PetscCall(PetscOptionsReal("-stg_h_scale_factor", "Scale element size for cutoff frequency calculation", NULL, stg_h_scale_factor,
260 &stg_h_scale_factor, NULL));
261 PetscCall(PetscOptionsDeprecated("-stg_dyScale", NULL, "libCEED 0.12.0", "Use -stg_h_scale_factor to scale all the element dimensions"));
262 PetscCall(PetscOptionsDeprecated("-stg_dz", NULL, "libCEED 0.12.0", NULL));
263 PetscOptionsEnd();
264
265 PetscCall(PetscNew(&global_stg_ctx));
266 *global_stg_ctx = (struct STGShur14Context_){
267 .alpha = alpha,
268 .u0 = u0,
269 .is_implicit = honee->phys->implicit,
270 .prescribe_T = prescribe_T,
271 .mean_only = mean_only,
272 .use_fluctuating_IC = use_fluctuating_IC,
273 .theta0 = theta0,
274 .P0 = P0,
275 .h_scale_factor = stg_h_scale_factor,
276 };
277
278 if (!use_stgstrong) { // Calculate dx assuming constant spacing
279 PetscReal domain_min[3], domain_max[3], domain_size[3];
280 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max));
281 for (PetscInt i = 0; i < 3; i++) domain_size[i] = domain_max[i] - domain_min[i];
282
283 PetscInt nmax = 3, faces[3];
284 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &nmax, NULL));
285 global_stg_ctx->dx = given_stg_dx ? stg_dx : domain_size[0] / faces[0];
286 PetscCheck((global_stg_ctx->dx > 0) && PetscIsNormalReal((PetscReal)global_stg_ctx->dx), comm, PETSC_ERR_LIB,
287 "STG dx must be positive normal number, got %g", global_stg_ctx->dx);
288 }
289
290 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfctx, CEED_MEM_HOST, &newtonian_ig_ctx));
291 global_stg_ctx->newt_ctx = *newtonian_ig_ctx;
292 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfctx, &newtonian_ig_ctx));
293
294 PetscCall(GetStgContextData(comm, dm, stg_inflow_path, stg_rand_path, &global_stg_ctx));
295
296 PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &stg_qfctx));
297 PetscCallCeed(ceed, CeedQFunctionContextSetData(stg_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, global_stg_ctx->total_bytes, global_stg_ctx));
298 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(stg_qfctx, CEED_MEM_HOST, FreeContextPetsc));
299 PetscCallCeed(ceed, CeedQFunctionContextRegisterDouble(stg_qfctx, "solution time", offsetof(struct STGShur14Context_, time), 1,
300 "Physical time of the solution"));
301
302 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&problem->ics.qfctx));
303 problem->ics = (HoneeQFSpec){.qf_func_ptr = ICsStg, .qf_loc = ICsStg_loc, .qfctx = stg_qfctx};
304
305 if (use_stgstrong) {
306 // Use default boundary integral QF (BoundaryIntegral) in newtonian.h
307 problem->use_strong_bc_ceed = PETSC_TRUE;
308 problem->set_bc_from_ics = PETSC_FALSE;
309 } else {
310 for (PetscCount b = 0; b < problem->num_bc_defs; b++) {
311 BCDefinition bc_def = problem->bc_defs[b];
312 const char *name;
313
314 PetscCall(BCDefinitionGetInfo(bc_def, &name, NULL, NULL));
315 if (!strcmp(name, "inflow")) {
316 HoneeBCStruct honee_bc;
317
318 PetscCall(PetscNew(&honee_bc));
319 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(stg_qfctx, &honee_bc->qfctx));
320 honee_bc->honee = honee;
321 honee_bc->num_comps_jac_data = honee->phys->implicit ? 11 : 0;
322 PetscCall(BCDefinitionSetContext(bc_def, (PetscCtxDestroyFn *)HoneeBCDestroy, honee_bc));
323
324 PetscCall(BCDefinitionSetIFunction(bc_def, STGWeakInflowBCSetup_CreateIFunctionQF, HoneeBCAddIFunctionOp));
325 PetscCall(BCDefinitionSetIJacobian(bc_def, STGWeakInflowBCSetup_CreateIJacobianQF, HoneeBCAddIJacobianOp));
326 }
327 }
328 problem->set_bc_from_ics = PETSC_TRUE;
329 }
330 PetscFunctionReturn(PETSC_SUCCESS);
331 }
332
333 // @brief Set STG strongly enforce components using DMAddBoundary
SetupStrongStg(DM dm,ProblemData problem,Physics phys)334 PetscErrorCode SetupStrongStg(DM dm, ProblemData problem, Physics phys) {
335 DMLabel label;
336 PetscInt comps[5], num_comps = 4;
337
338 PetscFunctionBeginUser;
339 switch (phys->state_var) {
340 case STATEVAR_CONSERVATIVE:
341 // {0,1,2,3} for rho, rho*u, rho*v, rho*w
342 for (int i = 0; i < 4; i++) comps[i] = i;
343 break;
344
345 case STATEVAR_PRIMITIVE:
346 // {1,2,3,4} for u, v, w, T
347 for (int i = 0; i < 4; i++) comps[i] = i + 1;
348 break;
349
350 case STATEVAR_ENTROPY:
351 // {1,2,3,4}
352 for (int i = 0; i < 4; i++) comps[i] = i + 1;
353 break;
354 }
355
356 PetscCall(DMGetLabel(dm, "Face Sets", &label));
357 for (PetscCount b = 0; b < problem->num_bc_defs; b++) {
358 BCDefinition bc_def = problem->bc_defs[b];
359 const char *name;
360
361 PetscCall(BCDefinitionGetInfo(bc_def, &name, NULL, NULL));
362 if (!strcmp(name, "inflow")) PetscCall(BCDefinitionSetEssential(bc_def, num_comps, comps));
363 }
364 PetscFunctionReturn(PETSC_SUCCESS);
365 }
366
SetupStrongStg_QF(Ceed ceed,ProblemData problem,CeedInt num_comp_x,CeedInt num_comp_q,CeedInt stg_data_size,CeedInt dXdx_size,CeedQFunction * qf_strongbc)367 PetscErrorCode SetupStrongStg_QF(Ceed ceed, ProblemData problem, CeedInt num_comp_x, CeedInt num_comp_q, CeedInt stg_data_size, CeedInt dXdx_size,
368 CeedQFunction *qf_strongbc) {
369 PetscFunctionBeginUser;
370 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, StgShur14InflowStrongQF, StgShur14InflowStrongQF_loc, qf_strongbc));
371 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_strongbc, "dXdx", dXdx_size, CEED_EVAL_NONE));
372 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_strongbc, "x", num_comp_x, CEED_EVAL_NONE));
373 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_strongbc, "scale", 1, CEED_EVAL_NONE));
374 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_strongbc, "stg data", stg_data_size, CEED_EVAL_NONE));
375 PetscCallCeed(ceed, CeedQFunctionAddOutput(*qf_strongbc, "q", num_comp_q, CEED_EVAL_NONE));
376
377 PetscCallCeed(ceed, CeedQFunctionSetContext(*qf_strongbc, problem->ics.qfctx));
378 PetscFunctionReturn(PETSC_SUCCESS);
379 }
380
SetupStrongStg_PreProcessing(Ceed ceed,ProblemData problem,CeedInt num_comp_x,CeedInt stg_data_size,CeedInt dXdx_size,CeedQFunction * qf_strongbc)381 PetscErrorCode SetupStrongStg_PreProcessing(Ceed ceed, ProblemData problem, CeedInt num_comp_x, CeedInt stg_data_size, CeedInt dXdx_size,
382 CeedQFunction *qf_strongbc) {
383 PetscFunctionBeginUser;
384 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, StgShur14Preprocess, StgShur14Preprocess_loc, qf_strongbc));
385 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_strongbc, "dXdx", dXdx_size, CEED_EVAL_NONE));
386 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_strongbc, "x", num_comp_x, CEED_EVAL_NONE));
387 PetscCallCeed(ceed, CeedQFunctionAddOutput(*qf_strongbc, "stg data", stg_data_size, CEED_EVAL_NONE));
388
389 PetscCallCeed(ceed, CeedQFunctionSetContext(*qf_strongbc, problem->ics.qfctx));
390 PetscFunctionReturn(PETSC_SUCCESS);
391 }
392