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