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