1 // Copyright (c) 2017-2022, 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 */ 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 if (isnan(Cij[0][i]) || isnan(Cij[1][i]) || isnan(Cij[2][i])) { 46 SETERRQ(comm, -1, "Cholesky decomposition failed at profile point %" PetscInt_FMT ". Either STGInflow has non-SPD matrix or contains nan.", 47 i + 1); 48 } 49 } 50 PetscFunctionReturn(0); 51 } 52 53 /* 54 * @brief Read the STGInflow file and load the contents into stg_ctx 55 * 56 * 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. 57 * Assumes there are 14 columns in the file. 58 * 59 * Function calculates the Cholesky decomposition from the Reynolds stress profile found in the file. 60 * 61 * @param[in] comm MPI_Comm for the program 62 * @param[in] path Path to the STGInflow.dat file 63 * @param[in,out] stg_ctx STGShur14Context where the data will be loaded into 64 */ 65 static PetscErrorCode ReadSTGInflow(const MPI_Comm comm, const char path[PETSC_MAX_PATH_LEN], STGShur14Context stg_ctx) { 66 PetscInt ndims, dims[2]; 67 FILE *fp; 68 const PetscInt char_array_len = 512; 69 char line[char_array_len]; 70 char **array; 71 72 PetscFunctionBeginUser; 73 74 PetscCall(PHASTADatFileOpen(comm, path, char_array_len, dims, &fp)); 75 76 CeedScalar rij[6][stg_ctx->nprofs]; 77 CeedScalar *wall_dist = &stg_ctx->data[stg_ctx->offsets.wall_dist]; 78 CeedScalar *eps = &stg_ctx->data[stg_ctx->offsets.eps]; 79 CeedScalar *lt = &stg_ctx->data[stg_ctx->offsets.lt]; 80 CeedScalar(*ubar)[stg_ctx->nprofs] = (CeedScalar(*)[stg_ctx->nprofs]) & stg_ctx->data[stg_ctx->offsets.ubar]; 81 82 for (PetscInt i = 0; i < stg_ctx->nprofs; i++) { 83 PetscCall(PetscSynchronizedFGets(comm, fp, char_array_len, line)); 84 PetscCall(PetscStrToArray(line, ' ', &ndims, &array)); 85 if (ndims < dims[1]) { 86 SETERRQ(comm, -1, "Line %" PetscInt_FMT " of %s does not contain enough columns (%" PetscInt_FMT " instead of %" PetscInt_FMT ")", i, path, 87 ndims, dims[1]); 88 } 89 90 wall_dist[i] = (CeedScalar)atof(array[0]); 91 ubar[0][i] = (CeedScalar)atof(array[1]); 92 ubar[1][i] = (CeedScalar)atof(array[2]); 93 ubar[2][i] = (CeedScalar)atof(array[3]); 94 rij[0][i] = (CeedScalar)atof(array[4]); 95 rij[1][i] = (CeedScalar)atof(array[5]); 96 rij[2][i] = (CeedScalar)atof(array[6]); 97 rij[3][i] = (CeedScalar)atof(array[7]); 98 rij[4][i] = (CeedScalar)atof(array[8]); 99 rij[5][i] = (CeedScalar)atof(array[9]); 100 lt[i] = (CeedScalar)atof(array[12]); 101 eps[i] = (CeedScalar)atof(array[13]); 102 103 if (wall_dist[i] < 0) SETERRQ(comm, -1, "Distance to wall in %s cannot be negative", path); 104 if (lt[i] < 0) SETERRQ(comm, -1, "Turbulent length scale in %s cannot be negative", path); 105 if (eps[i] < 0) SETERRQ(comm, -1, "Turbulent dissipation in %s cannot be negative", path); 106 } 107 CeedScalar(*cij)[stg_ctx->nprofs] = (CeedScalar(*)[stg_ctx->nprofs]) & stg_ctx->data[stg_ctx->offsets.cij]; 108 PetscCall(CalcCholeskyDecomp(comm, stg_ctx->nprofs, rij, cij)); 109 PetscCall(PetscFClose(comm, fp)); 110 111 PetscFunctionReturn(0); 112 } 113 114 /* 115 * @brief Read the STGRand file and load the contents into stg_ctx 116 * 117 * 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. 118 * Assumes there are 7 columns in the file. 119 * 120 * @param[in] comm MPI_Comm for the program 121 * @param[in] path Path to the STGRand.dat file 122 * @param[in,out] stg_ctx STGShur14Context where the data will be loaded into 123 */ 124 static PetscErrorCode ReadSTGRand(const MPI_Comm comm, const char path[PETSC_MAX_PATH_LEN], STGShur14Context stg_ctx) { 125 PetscInt ndims, dims[2]; 126 FILE *fp; 127 const PetscInt char_array_len = 512; 128 char line[char_array_len]; 129 char **array; 130 131 PetscFunctionBeginUser; 132 PetscCall(PHASTADatFileOpen(comm, path, char_array_len, dims, &fp)); 133 134 CeedScalar *phi = &stg_ctx->data[stg_ctx->offsets.phi]; 135 CeedScalar(*d)[stg_ctx->nmodes] = (CeedScalar(*)[stg_ctx->nmodes]) & stg_ctx->data[stg_ctx->offsets.d]; 136 CeedScalar(*sigma)[stg_ctx->nmodes] = (CeedScalar(*)[stg_ctx->nmodes]) & stg_ctx->data[stg_ctx->offsets.sigma]; 137 138 for (PetscInt i = 0; i < stg_ctx->nmodes; i++) { 139 PetscCall(PetscSynchronizedFGets(comm, fp, char_array_len, line)); 140 PetscCall(PetscStrToArray(line, ' ', &ndims, &array)); 141 if (ndims < dims[1]) { 142 SETERRQ(comm, -1, "Line %" PetscInt_FMT " of %s does not contain enough columns (%" PetscInt_FMT " instead of %" PetscInt_FMT ")", i, path, 143 ndims, dims[1]); 144 } 145 146 d[0][i] = (CeedScalar)atof(array[0]); 147 d[1][i] = (CeedScalar)atof(array[1]); 148 d[2][i] = (CeedScalar)atof(array[2]); 149 phi[i] = (CeedScalar)atof(array[3]); 150 sigma[0][i] = (CeedScalar)atof(array[4]); 151 sigma[1][i] = (CeedScalar)atof(array[5]); 152 sigma[2][i] = (CeedScalar)atof(array[6]); 153 } 154 PetscCall(PetscFClose(comm, fp)); 155 156 PetscFunctionReturn(0); 157 } 158 159 /* 160 * @brief Read STG data from input paths and put in STGShur14Context 161 * 162 * Reads data from input paths and puts them into a STGShur14Context object. 163 * Data stored initially in `*pstg_ctx` will be copied over to the new STGShur14Context instance. 164 * 165 * @param[in] comm MPI_Comm for the program 166 * @param[in] dm DM for the program 167 * @param[in] stg_inflow_path Path to STGInflow.dat file 168 * @param[in] stg_rand_path Path to STGRand.dat file 169 * @param[in,out] pstg_ctx Pointer to STGShur14Context where the data will be loaded into 170 */ 171 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], 172 STGShur14Context *pstg_ctx, const CeedScalar ynodes[]) { 173 PetscInt nmodes, nprofs; 174 STGShur14Context stg_ctx; 175 PetscFunctionBeginUser; 176 177 // Get options 178 PetscCall(PHASTADatFileGetNRows(comm, stg_rand_path, &nmodes)); 179 PetscCall(PHASTADatFileGetNRows(comm, stg_inflow_path, &nprofs)); 180 if (nmodes > STG_NMODES_MAX) { 181 SETERRQ(comm, 1, 182 "Number of wavemodes in %s (%" PetscInt_FMT ") exceeds STG_NMODES_MAX (%" PetscInt_FMT 183 "). " 184 "Change size of STG_NMODES_MAX and recompile", 185 stg_rand_path, nmodes, STG_NMODES_MAX); 186 } 187 188 { 189 STGShur14Context s; 190 PetscCall(PetscCalloc1(1, &s)); 191 *s = **pstg_ctx; 192 s->nmodes = nmodes; 193 s->nprofs = nprofs; 194 s->offsets.sigma = 0; 195 s->offsets.d = nmodes * 3; 196 s->offsets.phi = s->offsets.d + nmodes * 3; 197 s->offsets.kappa = s->offsets.phi + nmodes; 198 s->offsets.wall_dist = s->offsets.kappa + nmodes; 199 s->offsets.ubar = s->offsets.wall_dist + nprofs; 200 s->offsets.cij = s->offsets.ubar + nprofs * 3; 201 s->offsets.eps = s->offsets.cij + nprofs * 6; 202 s->offsets.lt = s->offsets.eps + nprofs; 203 s->offsets.ynodes = s->offsets.lt + nprofs; 204 PetscInt total_num_scalars = s->offsets.ynodes + s->nynodes; 205 s->total_bytes = sizeof(*stg_ctx) + total_num_scalars * sizeof(stg_ctx->data[0]); 206 PetscCall(PetscMalloc(s->total_bytes, &stg_ctx)); 207 *stg_ctx = *s; 208 PetscCall(PetscFree(s)); 209 } 210 211 PetscCall(ReadSTGInflow(comm, stg_inflow_path, stg_ctx)); 212 PetscCall(ReadSTGRand(comm, stg_rand_path, stg_ctx)); 213 214 if (stg_ctx->nynodes > 0) { 215 CeedScalar *ynodes_ctx = &stg_ctx->data[stg_ctx->offsets.ynodes]; 216 for (PetscInt i = 0; i < stg_ctx->nynodes; i++) ynodes_ctx[i] = ynodes[i]; 217 } 218 219 // -- Calculate kappa 220 { 221 CeedScalar *kappa = &stg_ctx->data[stg_ctx->offsets.kappa]; 222 CeedScalar *wall_dist = &stg_ctx->data[stg_ctx->offsets.wall_dist]; 223 CeedScalar *lt = &stg_ctx->data[stg_ctx->offsets.lt]; 224 CeedScalar le, le_max = 0; 225 226 CeedPragmaSIMD for (PetscInt i = 0; i < stg_ctx->nprofs; i++) { 227 le = PetscMin(2 * wall_dist[i], 3 * lt[i]); 228 if (le_max < le) le_max = le; 229 } 230 CeedScalar kmin = M_PI / le_max; 231 232 CeedPragmaSIMD for (PetscInt i = 0; i < stg_ctx->nmodes; i++) { kappa[i] = kmin * pow(stg_ctx->alpha, i); } 233 } // end calculate kappa 234 235 PetscCall(PetscFree(*pstg_ctx)); 236 *pstg_ctx = stg_ctx; 237 PetscFunctionReturn(0); 238 } 239 240 PetscErrorCode SetupSTG(const MPI_Comm comm, const DM dm, ProblemData *problem, User user, const bool prescribe_T, const CeedScalar theta0, 241 const CeedScalar P0, const CeedScalar ynodes[], const CeedInt nynodes) { 242 char stg_inflow_path[PETSC_MAX_PATH_LEN] = "./STGInflow.dat"; 243 char stg_rand_path[PETSC_MAX_PATH_LEN] = "./STGRand.dat"; 244 PetscBool mean_only = PETSC_FALSE, use_stgstrong = PETSC_FALSE, use_fluctuating_IC = PETSC_FALSE; 245 CeedScalar u0 = 0.0, alpha = 1.01; 246 CeedQFunctionContext stg_context; 247 NewtonianIdealGasContext newtonian_ig_ctx; 248 PetscFunctionBeginUser; 249 250 // Get options 251 PetscOptionsBegin(comm, NULL, "STG Boundary Condition Options", NULL); 252 PetscCall(PetscOptionsString("-stg_inflow_path", "Path to STGInflow.dat", NULL, stg_inflow_path, stg_inflow_path, sizeof(stg_inflow_path), NULL)); 253 PetscCall(PetscOptionsString("-stg_rand_path", "Path to STGInflow.dat", NULL, stg_rand_path, stg_rand_path, sizeof(stg_rand_path), NULL)); 254 PetscCall(PetscOptionsReal("-stg_alpha", "Growth rate of the wavemodes", NULL, alpha, &alpha, NULL)); 255 PetscCall(PetscOptionsReal("-stg_u0", "Advective velocity for the fluctuations", NULL, u0, &u0, NULL)); 256 PetscCall(PetscOptionsBool("-stg_mean_only", "Only apply mean profile", NULL, mean_only, &mean_only, NULL)); 257 PetscCall(PetscOptionsBool("-stg_strong", "Enforce STG inflow strongly", NULL, use_stgstrong, &use_stgstrong, NULL)); 258 PetscCall(PetscOptionsBool("-stg_fluctuating_IC", "\"Extrude\" the fluctuations through the domain as an initial condition", NULL, 259 use_fluctuating_IC, &use_fluctuating_IC, 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 = user->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->nynodes = nynodes; 272 273 { 274 // Calculate dx assuming constant spacing 275 PetscReal domain_min[3], domain_max[3], domain_size[3]; 276 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max)); 277 for (PetscInt i = 0; i < 3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 278 279 PetscInt nmax = 3, faces[3]; 280 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &nmax, NULL)); 281 global_stg_ctx->dx = domain_size[0] / faces[0]; 282 global_stg_ctx->dz = domain_size[2] / faces[2]; 283 } 284 285 CeedQFunctionContextGetData(problem->apply_vol_rhs.qfunction_context, CEED_MEM_HOST, &newtonian_ig_ctx); 286 global_stg_ctx->newtonian_ctx = *newtonian_ig_ctx; 287 CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfunction_context, &newtonian_ig_ctx); 288 289 PetscCall(GetSTGContextData(comm, dm, stg_inflow_path, stg_rand_path, &global_stg_ctx, ynodes)); 290 291 CeedQFunctionContextCreate(user->ceed, &stg_context); 292 CeedQFunctionContextSetData(stg_context, CEED_MEM_HOST, CEED_USE_POINTER, global_stg_ctx->total_bytes, global_stg_ctx); 293 CeedQFunctionContextSetDataDestroy(stg_context, CEED_MEM_HOST, FreeContextPetsc); 294 CeedQFunctionContextRegisterDouble(stg_context, "solution time", offsetof(struct STGShur14Context_, time), 1, "Physical time of the solution"); 295 296 CeedQFunctionContextDestroy(&problem->ics.qfunction_context); 297 problem->ics.qfunction = ICsSTG; 298 problem->ics.qfunction_loc = ICsSTG_loc; 299 problem->ics.qfunction_context = stg_context; 300 301 if (use_stgstrong) { 302 // Use default boundary integral QF (BoundaryIntegral) in newtonian.h 303 problem->use_dirichlet_ceed = PETSC_TRUE; 304 problem->bc_from_ics = PETSC_FALSE; 305 } else { 306 problem->apply_inflow.qfunction = STGShur14_Inflow; 307 problem->apply_inflow.qfunction_loc = STGShur14_Inflow_loc; 308 problem->apply_inflow_jacobian.qfunction = STGShur14_Inflow_Jacobian; 309 problem->apply_inflow_jacobian.qfunction_loc = STGShur14_Inflow_Jacobian_loc; 310 CeedQFunctionContextReferenceCopy(stg_context, &problem->apply_inflow.qfunction_context); 311 CeedQFunctionContextReferenceCopy(stg_context, &problem->apply_inflow_jacobian.qfunction_context); 312 problem->bc_from_ics = PETSC_TRUE; 313 } 314 315 PetscFunctionReturn(0); 316 } 317 318 static inline PetscScalar FindDy(const PetscScalar ynodes[], const PetscInt nynodes, const PetscScalar y) { 319 const PetscScalar half_mindy = 0.5 * (ynodes[1] - ynodes[0]); 320 // ^^assuming min(dy) is first element off the wall 321 PetscInt idx = -1; // Index 322 323 for (PetscInt i = 0; i < nynodes; i++) { 324 if (y < ynodes[i] + half_mindy) { 325 idx = i; 326 break; 327 } 328 } 329 if (idx == 0) return ynodes[1] - ynodes[0]; 330 else if (idx == nynodes - 1) return ynodes[nynodes - 2] - ynodes[nynodes - 1]; 331 else return 0.5 * (ynodes[idx + 1] - ynodes[idx - 1]); 332 } 333 334 // Function passed to DMAddBoundary 335 // NOTE: Not used in favor of QFunction-based method 336 PetscErrorCode StrongSTGbcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[], void *ctx) { 337 PetscFunctionBeginUser; 338 339 const STGShur14Context stg_ctx = (STGShur14Context)ctx; 340 PetscScalar qn[stg_ctx->nmodes], u[3], ubar[3], cij[6], eps, lt; 341 const bool mean_only = stg_ctx->mean_only; 342 const PetscScalar dx = stg_ctx->dx; 343 const PetscScalar dz = stg_ctx->dz; 344 const PetscScalar mu = stg_ctx->newtonian_ctx.mu; 345 const PetscScalar theta0 = stg_ctx->theta0; 346 const PetscScalar P0 = stg_ctx->P0; 347 const PetscScalar cv = stg_ctx->newtonian_ctx.cv; 348 const PetscScalar cp = stg_ctx->newtonian_ctx.cp; 349 const PetscScalar Rd = cp - cv; 350 351 const CeedScalar rho = P0 / (Rd * theta0); 352 InterpolateProfile(x[1], ubar, cij, &eps, <, stg_ctx); 353 if (!mean_only) { 354 const PetscInt nynodes = stg_ctx->nynodes; 355 const PetscScalar *ynodes = &stg_ctx->data[stg_ctx->offsets.ynodes]; 356 const PetscScalar h[3] = {dx, FindDy(ynodes, nynodes, x[1]), dz}; 357 CalcSpectrum(x[1], eps, lt, h, mu / rho, qn, stg_ctx); 358 STGShur14_Calc(x, time, ubar, cij, qn, u, stg_ctx); 359 } else { 360 for (CeedInt j = 0; j < 3; j++) u[j] = ubar[j]; 361 } 362 363 bcval[0] = rho; 364 bcval[1] = rho * u[0]; 365 bcval[2] = rho * u[1]; 366 bcval[3] = rho * u[2]; 367 PetscFunctionReturn(0); 368 } 369 370 PetscErrorCode SetupStrongSTG(DM dm, SimpleBC bc, ProblemData *problem, Physics phys) { 371 DMLabel label; 372 PetscFunctionBeginUser; 373 374 PetscInt comps[5], num_comps = 4; 375 switch (phys->state_var) { 376 case STATEVAR_CONSERVATIVE: 377 // {0,1,2,3} for rho, rho*u, rho*v, rho*w 378 for (int i = 0; i < 4; i++) comps[i] = i; 379 break; 380 381 case STATEVAR_PRIMITIVE: 382 // {1,2,3,4} for u, v, w, T 383 for (int i = 0; i < 4; i++) comps[i] = i + 1; 384 break; 385 } 386 387 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 388 // Set wall BCs 389 if (bc->num_inflow > 0) { 390 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "STG", label, bc->num_inflow, bc->inflows, 0, num_comps, comps, (void (*)(void))StrongSTGbcFunc, 391 NULL, global_stg_ctx, NULL)); 392 } 393 394 PetscFunctionReturn(0); 395 } 396 397 PetscErrorCode SetupStrongSTG_QF(Ceed ceed, ProblemData *problem, CeedInt num_comp_x, CeedInt num_comp_q, CeedInt stg_data_size, 398 CeedInt q_data_size_sur, CeedQFunction *pqf_strongbc) { 399 CeedQFunction qf_strongbc; 400 PetscFunctionBeginUser; 401 CeedQFunctionCreateInterior(ceed, 1, STGShur14_Inflow_StrongQF, STGShur14_Inflow_StrongQF_loc, &qf_strongbc); 402 CeedQFunctionAddInput(qf_strongbc, "surface qdata", q_data_size_sur, CEED_EVAL_NONE); 403 CeedQFunctionAddInput(qf_strongbc, "x", num_comp_x, CEED_EVAL_NONE); 404 CeedQFunctionAddInput(qf_strongbc, "scale", 1, CEED_EVAL_NONE); 405 CeedQFunctionAddInput(qf_strongbc, "stg data", stg_data_size, CEED_EVAL_NONE); 406 CeedQFunctionAddOutput(qf_strongbc, "q", num_comp_q, CEED_EVAL_NONE); 407 408 CeedQFunctionSetContext(qf_strongbc, problem->ics.qfunction_context); 409 *pqf_strongbc = qf_strongbc; 410 PetscFunctionReturn(0); 411 } 412 413 PetscErrorCode SetupStrongSTG_PreProcessing(Ceed ceed, ProblemData *problem, CeedInt num_comp_x, CeedInt stg_data_size, CeedInt q_data_size_sur, 414 CeedQFunction *pqf_strongbc) { 415 CeedQFunction qf_strongbc; 416 PetscFunctionBeginUser; 417 CeedQFunctionCreateInterior(ceed, 1, Preprocess_STGShur14, Preprocess_STGShur14_loc, &qf_strongbc); 418 CeedQFunctionAddInput(qf_strongbc, "surface qdata", q_data_size_sur, CEED_EVAL_NONE); 419 CeedQFunctionAddInput(qf_strongbc, "x", num_comp_x, CEED_EVAL_NONE); 420 CeedQFunctionAddOutput(qf_strongbc, "stg data", stg_data_size, CEED_EVAL_NONE); 421 422 CeedQFunctionSetContext(qf_strongbc, problem->ics.qfunction_context); 423 *pqf_strongbc = qf_strongbc; 424 PetscFunctionReturn(0); 425 } 426