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 <stdlib.h> 13 #include <math.h> 14 #include <petsc.h> 15 #include "../navierstokes.h" 16 #include "stg_shur14.h" 17 #include "../qfunctions/stg_shur14.h" 18 19 STGShur14Context global_stg_ctx; 20 21 /* 22 * @brief Perform Cholesky decomposition on array of symmetric 3x3 matrices 23 * 24 * This assumes the input matrices are in order [11,22,33,12,13,23]. This 25 * format is also used for the output. 26 * 27 * @param[in] comm MPI_Comm 28 * @param[in] nprofs Number of matrices in Rij 29 * @param[in] Rij Array of the symmetric matrices [6,nprofs] 30 * @param[out] Cij Array of the Cholesky Decomposition matrices, [6,nprofs] 31 */ 32 PetscErrorCode CalcCholeskyDecomp(MPI_Comm comm, PetscInt nprofs, 33 const CeedScalar Rij[6][nprofs], CeedScalar Cij[6][nprofs]) { 34 PetscFunctionBeginUser; 35 for (PetscInt i=0; i<nprofs; i++) { 36 Cij[0][i] = sqrt(Rij[0][i]); 37 Cij[3][i] = Rij[3][i] / Cij[0][i]; 38 Cij[1][i] = sqrt(Rij[1][i] - pow(Cij[3][i], 2) ); 39 Cij[4][i] = Rij[4][i] / Cij[0][i]; 40 Cij[5][i] = (Rij[5][i] - Cij[3][i]*Cij[4][i]) / Cij[1][i]; 41 Cij[2][i] = sqrt(Rij[2][i] - pow(Cij[4][i], 2) - pow(Cij[5][i], 2)); 42 43 if (isnan(Cij[0][i]) || isnan(Cij[1][i]) || isnan(Cij[2][i])) 44 SETERRQ(comm, -1, "Cholesky decomposition failed at profile point %" 45 PetscInt_FMT ". Either STGInflow has non-SPD matrix or contains nan.", i+1); 46 } 47 PetscFunctionReturn(0); 48 } 49 50 51 /* 52 * @brief Open a PHASTA *.dat file, grabbing dimensions and file pointer 53 * 54 * This function opens the file specified by `path` using `PetscFOpen` and 55 * passes the file pointer in `fp`. It is not closed in this function, thus 56 * `fp` must be closed sometime after this function has been called (using 57 * `PetscFClose` for example). 58 * 59 * Assumes that the first line of the file has the number of rows and columns 60 * as the only two entries, separated by a single space 61 * 62 * @param[in] comm MPI_Comm for the program 63 * @param[in] path Path to the file 64 * @param[in] char_array_len Length of the character array that should contain each line 65 * @param[out] dims Dimensions of the file, taken from the first line of the file 66 * @param[out] fp File pointer to the opened file 67 */ 68 static PetscErrorCode OpenPHASTADatFile(const MPI_Comm comm, 69 const char path[PETSC_MAX_PATH_LEN], const PetscInt char_array_len, 70 PetscInt dims[2], FILE **fp) { 71 PetscInt ndims; 72 char line[char_array_len]; 73 char **array; 74 75 PetscFunctionBeginUser; 76 PetscCall(PetscFOpen(comm, path, "r", fp)); 77 PetscCall(PetscSynchronizedFGets(comm, *fp, char_array_len, line)); 78 PetscCall(PetscStrToArray(line, ' ', &ndims, &array)); 79 if (ndims != 2) SETERRQ(comm, -1, 80 "Found %" PetscInt_FMT" dimensions instead of 2 on the first line of %s", 81 ndims, path); 82 83 for (PetscInt i=0; i<ndims; i++) dims[i] = atoi(array[i]); 84 PetscCall(PetscStrToArrayDestroy(ndims, array)); 85 PetscFunctionReturn(0); 86 } 87 88 /* 89 * @brief Get the number of rows for the PHASTA file at path 90 * 91 * Assumes that the first line of the file has the number of rows and columns 92 * as the only two entries, separated by a single space 93 * 94 * @param[in] comm MPI_Comm for the program 95 * @param[in] path Path to the file 96 * @param[out] nrows Number of rows 97 */ 98 static PetscErrorCode GetNRows(const MPI_Comm comm, 99 const char path[PETSC_MAX_PATH_LEN], PetscInt *nrows) { 100 const PetscInt char_array_len = 512; 101 PetscInt dims[2]; 102 FILE *fp; 103 104 PetscFunctionBeginUser; 105 PetscCall(OpenPHASTADatFile(comm, path, char_array_len, dims, &fp)); 106 *nrows = dims[0]; 107 PetscCall(PetscFClose(comm, fp)); 108 PetscFunctionReturn(0); 109 } 110 111 /* 112 * @brief Read the STGInflow file and load the contents into stg_ctx 113 * 114 * Assumes that the first line of the file has the number of rows and columns 115 * as the only two entries, separated by a single space. 116 * Assumes there are 14 columns in the file 117 * 118 * Function calculates the Cholesky decomposition from the Reynolds stress 119 * profile found in the file 120 * 121 * @param[in] comm MPI_Comm for the program 122 * @param[in] path Path to the STGInflow.dat file 123 * @param[inout] stg_ctx STGShur14Context where the data will be loaded into 124 */ 125 static PetscErrorCode ReadSTGInflow(const MPI_Comm comm, 126 const char path[PETSC_MAX_PATH_LEN], STGShur14Context stg_ctx) { 127 PetscInt ndims, dims[2]; 128 FILE *fp; 129 const PetscInt char_array_len=512; 130 char line[char_array_len]; 131 char **array; 132 133 PetscFunctionBeginUser; 134 135 PetscCall(OpenPHASTADatFile(comm, path, char_array_len, dims, &fp)); 136 137 CeedScalar rij[6][stg_ctx->nprofs]; 138 CeedScalar *wall_dist = &stg_ctx->data[stg_ctx->offsets.wall_dist]; 139 CeedScalar *eps = &stg_ctx->data[stg_ctx->offsets.eps]; 140 CeedScalar *lt = &stg_ctx->data[stg_ctx->offsets.lt]; 141 CeedScalar (*ubar)[stg_ctx->nprofs] = (CeedScalar (*)[stg_ctx->nprofs]) 142 &stg_ctx->data[stg_ctx->offsets.ubar]; 143 144 for (PetscInt i=0; i<stg_ctx->nprofs; i++) { 145 PetscCall(PetscSynchronizedFGets(comm, fp, char_array_len, line)); 146 PetscCall(PetscStrToArray(line, ' ', &ndims, &array)); 147 if (ndims < dims[1]) SETERRQ(comm, -1, 148 "Line %" PetscInt_FMT" of %s does not contain enough columns (%" 149 PetscInt_FMT" instead of %" PetscInt_FMT")", i, 150 path, ndims, dims[1]); 151 152 wall_dist[i] = (CeedScalar) atof(array[0]); 153 ubar[0][i] = (CeedScalar) atof(array[1]); 154 ubar[1][i] = (CeedScalar) atof(array[2]); 155 ubar[2][i] = (CeedScalar) atof(array[3]); 156 rij[0][i] = (CeedScalar) atof(array[4]); 157 rij[1][i] = (CeedScalar) atof(array[5]); 158 rij[2][i] = (CeedScalar) atof(array[6]); 159 rij[3][i] = (CeedScalar) atof(array[7]); 160 rij[4][i] = (CeedScalar) atof(array[8]); 161 rij[5][i] = (CeedScalar) atof(array[9]); 162 lt[i] = (CeedScalar) atof(array[12]); 163 eps[i] = (CeedScalar) atof(array[13]); 164 165 if (wall_dist[i] < 0) SETERRQ(comm, -1, 166 "Distance to wall in %s cannot be negative", path); 167 if (lt[i] < 0) SETERRQ(comm, -1, 168 "Turbulent length scale in %s cannot be negative", path); 169 if (eps[i] < 0) SETERRQ(comm, -1, 170 "Turbulent dissipation in %s cannot be negative", path); 171 172 } 173 CeedScalar (*cij)[stg_ctx->nprofs] = (CeedScalar (*)[stg_ctx->nprofs]) 174 &stg_ctx->data[stg_ctx->offsets.cij]; 175 PetscCall(CalcCholeskyDecomp(comm, stg_ctx->nprofs, rij, cij)); 176 PetscCall(PetscFClose(comm, fp)); 177 PetscFunctionReturn(0); 178 } 179 180 /* 181 * @brief Read the STGRand file and load the contents into stg_ctx 182 * 183 * Assumes that the first line of the file has the number of rows and columns 184 * as the only two entries, separated by a single space. 185 * Assumes there are 7 columns in the file 186 * 187 * @param[in] comm MPI_Comm for the program 188 * @param[in] path Path to the STGRand.dat file 189 * @param[inout] stg_ctx STGShur14Context where the data will be loaded into 190 */ 191 static PetscErrorCode ReadSTGRand(const MPI_Comm comm, 192 const char path[PETSC_MAX_PATH_LEN], 193 STGShur14Context stg_ctx) { 194 PetscInt ndims, dims[2]; 195 FILE *fp; 196 const PetscInt char_array_len = 512; 197 char line[char_array_len]; 198 char **array; 199 200 PetscFunctionBeginUser; 201 PetscCall(OpenPHASTADatFile(comm, path, char_array_len, dims, &fp)); 202 203 CeedScalar *phi = &stg_ctx->data[stg_ctx->offsets.phi]; 204 CeedScalar (*d)[stg_ctx->nmodes] = (CeedScalar (*)[stg_ctx->nmodes]) 205 &stg_ctx->data[stg_ctx->offsets.d]; 206 CeedScalar (*sigma)[stg_ctx->nmodes] = (CeedScalar (*)[stg_ctx->nmodes]) 207 &stg_ctx->data[stg_ctx->offsets.sigma]; 208 209 for (PetscInt i=0; i<stg_ctx->nmodes; i++) { 210 PetscCall(PetscSynchronizedFGets(comm, fp, char_array_len, line)); 211 PetscCall(PetscStrToArray(line, ' ', &ndims, &array)); 212 if (ndims < dims[1]) SETERRQ(comm, -1, 213 "Line %" PetscInt_FMT" of %s does not contain enough columns (%" 214 PetscInt_FMT" instead of %" PetscInt_FMT")", i, 215 path, ndims, dims[1]); 216 217 d[0][i] = (CeedScalar) atof(array[0]); 218 d[1][i] = (CeedScalar) atof(array[1]); 219 d[2][i] = (CeedScalar) atof(array[2]); 220 phi[i] = (CeedScalar) atof(array[3]); 221 sigma[0][i] = (CeedScalar) atof(array[4]); 222 sigma[1][i] = (CeedScalar) atof(array[5]); 223 sigma[2][i] = (CeedScalar) atof(array[6]); 224 } 225 PetscCall(PetscFClose(comm, fp)); 226 PetscFunctionReturn(0); 227 } 228 229 /* 230 * @brief Read STG data from input paths and put in STGShur14Context 231 * 232 * Reads data from input paths and puts them into a STGShur14Context object. 233 * Data stored initially in `*pstg_ctx` will be copied over to the new 234 * STGShur14Context instance. 235 * 236 * @param[in] comm MPI_Comm for the program 237 * @param[in] dm DM for the program 238 * @param[in] stg_inflow_path Path to STGInflow.dat file 239 * @param[in] stg_rand_path Path to STGRand.dat file 240 * @param[inout] pstg_ctx Pointer to STGShur14Context where the data will be loaded into 241 */ 242 PetscErrorCode GetSTGContextData(const MPI_Comm comm, const DM dm, 243 char stg_inflow_path[PETSC_MAX_PATH_LEN], 244 char stg_rand_path[PETSC_MAX_PATH_LEN], 245 STGShur14Context *pstg_ctx, 246 const CeedScalar ynodes[]) { 247 PetscInt nmodes, nprofs; 248 STGShur14Context stg_ctx; 249 PetscFunctionBeginUser; 250 251 // Get options 252 PetscCall(GetNRows(comm, stg_rand_path, &nmodes)); 253 PetscCall(GetNRows(comm, stg_inflow_path, &nprofs)); 254 if (nmodes > STG_NMODES_MAX) 255 SETERRQ(comm, 1, "Number of wavemodes in %s (%" 256 PetscInt_FMT") exceeds STG_NMODES_MAX (%" PetscInt_FMT"). " 257 "Change size of STG_NMODES_MAX and recompile", stg_rand_path, nmodes, 258 STG_NMODES_MAX); 259 260 { 261 STGShur14Context s; 262 PetscCall(PetscCalloc1(1, &s)); 263 *s = **pstg_ctx; 264 s->nmodes = nmodes; 265 s->nprofs = nprofs; 266 s->offsets.sigma = 0; 267 s->offsets.d = nmodes*3; 268 s->offsets.phi = s->offsets.d + nmodes*3; 269 s->offsets.kappa = s->offsets.phi + nmodes; 270 s->offsets.wall_dist = s->offsets.kappa + nmodes; 271 s->offsets.ubar = s->offsets.wall_dist + nprofs; 272 s->offsets.cij = s->offsets.ubar + nprofs*3; 273 s->offsets.eps = s->offsets.cij + nprofs*6; 274 s->offsets.lt = s->offsets.eps + nprofs; 275 s->offsets.ynodes = s->offsets.lt + nprofs; 276 PetscInt total_num_scalars = s->offsets.ynodes + s->nynodes; 277 s->total_bytes = sizeof(*stg_ctx) + total_num_scalars*sizeof(stg_ctx->data[0]); 278 PetscCall(PetscMalloc(s->total_bytes, &stg_ctx)); 279 *stg_ctx = *s; 280 PetscCall(PetscFree(s)); 281 } 282 283 PetscCall(ReadSTGInflow(comm, stg_inflow_path, stg_ctx)); 284 PetscCall(ReadSTGRand(comm, stg_rand_path, stg_ctx)); 285 286 if (stg_ctx->nynodes > 0) { 287 CeedScalar *ynodes_ctx = &stg_ctx->data[stg_ctx->offsets.ynodes]; 288 for (PetscInt i=0; i<stg_ctx->nynodes; i++) ynodes_ctx[i] = ynodes[i]; 289 } 290 291 // -- Calculate kappa 292 { 293 CeedScalar *kappa = &stg_ctx->data[stg_ctx->offsets.kappa]; 294 CeedScalar *wall_dist = &stg_ctx->data[stg_ctx->offsets.wall_dist]; 295 CeedScalar *lt = &stg_ctx->data[stg_ctx->offsets.lt]; 296 CeedScalar le, le_max=0; 297 298 CeedPragmaSIMD 299 for (PetscInt i=0; i<stg_ctx->nprofs; i++) { 300 le = PetscMin(2*wall_dist[i], 3*lt[i]); 301 if (le_max < le) le_max = le; 302 } 303 CeedScalar kmin = M_PI/le_max; 304 305 CeedPragmaSIMD 306 for (PetscInt i=0; i<stg_ctx->nmodes; i++) { 307 kappa[i] = kmin*pow(stg_ctx->alpha, i); 308 } 309 } //end calculate kappa 310 311 PetscCall(PetscFree(*pstg_ctx)); 312 *pstg_ctx = stg_ctx; 313 PetscFunctionReturn(0); 314 } 315 316 PetscErrorCode SetupSTG(const MPI_Comm comm, const DM dm, ProblemData *problem, 317 User user, const bool prescribe_T, 318 const CeedScalar theta0, const CeedScalar P0, 319 const CeedScalar ynodes[], const CeedInt nynodes) { 320 char stg_inflow_path[PETSC_MAX_PATH_LEN] = "./STGInflow.dat"; 321 char stg_rand_path[PETSC_MAX_PATH_LEN] = "./STGRand.dat"; 322 PetscBool mean_only = PETSC_FALSE, 323 use_stgstrong = PETSC_FALSE, 324 use_fluctuating_IC = PETSC_FALSE; 325 CeedScalar u0 = 0.0, 326 alpha = 1.01; 327 CeedQFunctionContext stg_context; 328 NewtonianIdealGasContext newtonian_ig_ctx; 329 PetscFunctionBeginUser; 330 331 // Get options 332 PetscOptionsBegin(comm, NULL, "STG Boundary Condition Options", NULL); 333 PetscCall(PetscOptionsString("-stg_inflow_path", "Path to STGInflow.dat", NULL, 334 stg_inflow_path, stg_inflow_path, 335 sizeof(stg_inflow_path), NULL)); 336 PetscCall(PetscOptionsString("-stg_rand_path", "Path to STGInflow.dat", NULL, 337 stg_rand_path,stg_rand_path, 338 sizeof(stg_rand_path), NULL)); 339 PetscCall(PetscOptionsReal("-stg_alpha", "Growth rate of the wavemodes", NULL, 340 alpha, &alpha, NULL)); 341 PetscCall(PetscOptionsReal("-stg_u0", "Advective velocity for the fluctuations", 342 NULL, u0, &u0, NULL)); 343 PetscCall(PetscOptionsBool("-stg_mean_only", "Only apply mean profile", 344 NULL, mean_only, &mean_only, NULL)); 345 PetscCall(PetscOptionsBool("-stg_strong", "Enforce STG inflow strongly", 346 NULL, use_stgstrong, &use_stgstrong, NULL)); 347 PetscCall(PetscOptionsBool("-stg_fluctuating_IC", 348 "\"Extrude\" the fluctuations through the domain as an initial condition", 349 NULL, use_fluctuating_IC, &use_fluctuating_IC, NULL)); 350 PetscOptionsEnd(); 351 352 PetscCall(PetscCalloc1(1, &global_stg_ctx)); 353 global_stg_ctx->alpha = alpha; 354 global_stg_ctx->u0 = u0; 355 global_stg_ctx->is_implicit = user->phys->implicit; 356 global_stg_ctx->prescribe_T = prescribe_T; 357 global_stg_ctx->mean_only = mean_only; 358 global_stg_ctx->use_fluctuating_IC = use_fluctuating_IC; 359 global_stg_ctx->theta0 = theta0; 360 global_stg_ctx->P0 = P0; 361 global_stg_ctx->nynodes = nynodes; 362 363 { 364 // Calculate dx assuming constant spacing 365 PetscReal domain_min[3], domain_max[3], domain_size[3]; 366 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max)); 367 for (PetscInt i=0; i<3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 368 369 PetscInt nmax = 3, faces[3]; 370 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, 371 &nmax, NULL)); 372 global_stg_ctx->dx = domain_size[0]/faces[0]; 373 global_stg_ctx->dz = domain_size[2]/faces[2]; 374 } 375 376 CeedQFunctionContextGetData(problem->apply_vol_rhs.qfunction_context, 377 CEED_MEM_HOST, &newtonian_ig_ctx); 378 global_stg_ctx->newtonian_ctx = *newtonian_ig_ctx; 379 CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfunction_context, 380 &newtonian_ig_ctx); 381 382 PetscCall(GetSTGContextData(comm, dm, stg_inflow_path, stg_rand_path, 383 &global_stg_ctx, ynodes)); 384 385 CeedQFunctionContextCreate(user->ceed, &stg_context); 386 CeedQFunctionContextSetData(stg_context, CEED_MEM_HOST, 387 CEED_USE_POINTER, global_stg_ctx->total_bytes, global_stg_ctx); 388 CeedQFunctionContextSetDataDestroy(stg_context, CEED_MEM_HOST, 389 FreeContextPetsc); 390 CeedQFunctionContextRegisterDouble(stg_context, "solution time", 391 offsetof(struct STGShur14Context_, time), 1, 392 "Physical time of the solution"); 393 394 CeedQFunctionContextDestroy(&problem->ics.qfunction_context); 395 problem->ics.qfunction = ICsSTG; 396 problem->ics.qfunction_loc = ICsSTG_loc; 397 problem->ics.qfunction_context = stg_context; 398 399 if (use_stgstrong) { 400 // Use default boundary integral QF (BoundaryIntegral) in newtonian.h 401 problem->use_dirichlet_ceed = PETSC_TRUE; 402 problem->bc_from_ics = PETSC_FALSE; 403 } else { 404 problem->apply_inflow.qfunction = STGShur14_Inflow; 405 problem->apply_inflow.qfunction_loc = STGShur14_Inflow_loc; 406 problem->apply_inflow_jacobian.qfunction = STGShur14_Inflow_Jacobian; 407 problem->apply_inflow_jacobian.qfunction_loc = STGShur14_Inflow_Jacobian_loc; 408 CeedQFunctionContextReferenceCopy(stg_context, 409 &problem->apply_inflow.qfunction_context); 410 CeedQFunctionContextReferenceCopy(stg_context, 411 &problem->apply_inflow_jacobian.qfunction_context); 412 problem->bc_from_ics = PETSC_TRUE; 413 } 414 415 PetscFunctionReturn(0); 416 } 417 418 static inline PetscScalar FindDy(const PetscScalar ynodes[], 419 const PetscInt nynodes, const PetscScalar y) { 420 const PetscScalar half_mindy = 0.5 * (ynodes[1] - ynodes[0]); 421 // ^^assuming min(dy) is first element off the wall 422 PetscInt idx = -1; // Index 423 424 for (PetscInt i=0; i<nynodes; i++) { 425 if (y < ynodes[i] + half_mindy) { 426 idx = i; break; 427 } 428 } 429 if (idx == 0) return ynodes[1] - ynodes[0]; 430 else if (idx == nynodes-1) return ynodes[nynodes-2] - ynodes[nynodes-1]; 431 else return 0.5 * (ynodes[idx+1] - ynodes[idx-1]); 432 } 433 434 // Function passed to DMAddBoundary 435 // NOTE: Not used in favor of QFunction-based method 436 PetscErrorCode StrongSTGbcFunc(PetscInt dim, PetscReal time, 437 const PetscReal x[], PetscInt Nc, PetscScalar bcval[], void *ctx) { 438 PetscFunctionBeginUser; 439 const STGShur14Context stg_ctx = (STGShur14Context) ctx; 440 PetscScalar qn[stg_ctx->nmodes], u[3], ubar[3], cij[6], eps, lt; 441 const bool mean_only = stg_ctx->mean_only; 442 const PetscScalar dx = stg_ctx->dx; 443 const PetscScalar dz = stg_ctx->dz; 444 const PetscScalar mu = stg_ctx->newtonian_ctx.mu; 445 const PetscScalar theta0 = stg_ctx->theta0; 446 const PetscScalar P0 = stg_ctx->P0; 447 const PetscScalar cv = stg_ctx->newtonian_ctx.cv; 448 const PetscScalar cp = stg_ctx->newtonian_ctx.cp; 449 const PetscScalar Rd = cp - cv; 450 451 const CeedScalar rho = P0 / (Rd * theta0); 452 InterpolateProfile(x[1], ubar, cij, &eps, <, stg_ctx); 453 if (!mean_only) { 454 const PetscInt nynodes = stg_ctx->nynodes; 455 const PetscScalar *ynodes = &stg_ctx->data[stg_ctx->offsets.ynodes]; 456 const PetscScalar h[3] = {dx, FindDy(ynodes, nynodes, x[1]), dz}; 457 CalcSpectrum(x[1], eps, lt, h, mu/rho, qn, stg_ctx); 458 STGShur14_Calc(x, time, ubar, cij, qn, u, stg_ctx); 459 } else { 460 for (CeedInt j=0; j<3; j++) u[j] = ubar[j]; 461 } 462 463 bcval[0] = rho; 464 bcval[1] = rho * u[0]; 465 bcval[2] = rho * u[1]; 466 bcval[3] = rho * u[2]; 467 PetscFunctionReturn(0); 468 } 469 470 PetscErrorCode SetupStrongSTG(DM dm, SimpleBC bc, ProblemData *problem, 471 Physics phys) { 472 DMLabel label; 473 PetscFunctionBeginUser; 474 475 PetscInt comps[5], num_comps=4; 476 switch (phys->state_var) { 477 case STATEVAR_CONSERVATIVE: 478 // {0,1,2,3} for rho, rho*u, rho*v, rho*w 479 for(int i=0; i<4; i++) comps[i] = i; 480 break; 481 482 case STATEVAR_PRIMITIVE: 483 // {1,2,3,4} for u, v, w, T 484 for(int i=0; i<4; i++) comps[i] = i+1; 485 break; 486 } 487 488 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 489 // Set wall BCs 490 if (bc->num_inflow > 0) { 491 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "STG", label, 492 bc->num_inflow, bc->inflows, 0, num_comps, 493 comps, (void(*)(void))StrongSTGbcFunc, 494 NULL, global_stg_ctx, NULL)); 495 } 496 497 PetscFunctionReturn(0); 498 } 499 500 PetscErrorCode SetupStrongSTG_QF(Ceed ceed, ProblemData *problem, 501 CeedInt num_comp_x, CeedInt num_comp_q, CeedInt stg_data_size, 502 CeedInt q_data_size_sur, CeedQFunction *pqf_strongbc) { 503 504 CeedQFunction qf_strongbc; 505 PetscFunctionBeginUser; 506 CeedQFunctionCreateInterior(ceed, 1, STGShur14_Inflow_StrongQF, 507 STGShur14_Inflow_StrongQF_loc, &qf_strongbc); 508 CeedQFunctionAddInput(qf_strongbc, "surface qdata", q_data_size_sur, 509 CEED_EVAL_NONE); 510 CeedQFunctionAddInput(qf_strongbc, "x", num_comp_x, CEED_EVAL_NONE); 511 CeedQFunctionAddInput(qf_strongbc, "scale", 1, CEED_EVAL_NONE); 512 CeedQFunctionAddInput(qf_strongbc, "stg data", stg_data_size, CEED_EVAL_NONE); 513 CeedQFunctionAddOutput(qf_strongbc, "q", num_comp_q, CEED_EVAL_NONE); 514 515 CeedQFunctionSetContext(qf_strongbc, problem->ics.qfunction_context); 516 *pqf_strongbc = qf_strongbc; 517 PetscFunctionReturn(0); 518 } 519 520 PetscErrorCode SetupStrongSTG_PreProcessing(Ceed ceed, ProblemData *problem, 521 CeedInt num_comp_x, CeedInt stg_data_size, CeedInt q_data_size_sur, 522 CeedQFunction *pqf_strongbc) { 523 524 CeedQFunction qf_strongbc; 525 PetscFunctionBeginUser; 526 CeedQFunctionCreateInterior(ceed, 1, Preprocess_STGShur14, 527 Preprocess_STGShur14_loc, &qf_strongbc); 528 CeedQFunctionAddInput(qf_strongbc, "surface qdata", q_data_size_sur, 529 CEED_EVAL_NONE); 530 CeedQFunctionAddInput(qf_strongbc, "x", num_comp_x, CEED_EVAL_NONE); 531 CeedQFunctionAddOutput(qf_strongbc, "stg data", stg_data_size, CEED_EVAL_NONE); 532 533 CeedQFunctionSetContext(qf_strongbc, problem->ics.qfunction_context); 534 *pqf_strongbc = qf_strongbc; 535 PetscFunctionReturn(0); 536 } 537