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 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 */ 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 73 PetscCall(PHASTADatFileOpen(comm, path, char_array_len, dims, &fp)); 74 75 CeedScalar rij[6][stg_ctx->nprofs]; 76 CeedScalar *wall_dist = &stg_ctx->data[stg_ctx->offsets.wall_dist]; 77 CeedScalar *eps = &stg_ctx->data[stg_ctx->offsets.eps]; 78 CeedScalar *lt = &stg_ctx->data[stg_ctx->offsets.lt]; 79 CeedScalar(*ubar)[stg_ctx->nprofs] = (CeedScalar(*)[stg_ctx->nprofs]) & stg_ctx->data[stg_ctx->offsets.ubar]; 80 81 for (PetscInt i = 0; i < stg_ctx->nprofs; i++) { 82 PetscCall(PetscSynchronizedFGets(comm, fp, char_array_len, line)); 83 PetscCall(PetscStrToArray(line, ' ', &ndims, &array)); 84 PetscCheck(ndims == dims[1], comm, PETSC_ERR_FILE_UNEXPECTED, 85 "Line %" PetscInt_FMT " of %s does not have correct number of columns (%d instead of %" PetscInt_FMT ")", i, path, ndims, dims[1]); 86 87 wall_dist[i] = (CeedScalar)atof(array[0]); 88 ubar[0][i] = (CeedScalar)atof(array[1]); 89 ubar[1][i] = (CeedScalar)atof(array[2]); 90 ubar[2][i] = (CeedScalar)atof(array[3]); 91 rij[0][i] = (CeedScalar)atof(array[4]); 92 rij[1][i] = (CeedScalar)atof(array[5]); 93 rij[2][i] = (CeedScalar)atof(array[6]); 94 rij[3][i] = (CeedScalar)atof(array[7]); 95 rij[4][i] = (CeedScalar)atof(array[8]); 96 rij[5][i] = (CeedScalar)atof(array[9]); 97 lt[i] = (CeedScalar)atof(array[12]); 98 eps[i] = (CeedScalar)atof(array[13]); 99 100 PetscCheck(wall_dist[i] >= 0, comm, PETSC_ERR_FILE_UNEXPECTED, "Distance to wall in %s cannot be negative", path); 101 PetscCheck(lt[i] >= 0, comm, PETSC_ERR_FILE_UNEXPECTED, "Turbulent length scale in %s cannot be negative", path); 102 PetscCheck(eps[i] >= 0, comm, PETSC_ERR_FILE_UNEXPECTED, "Turbulent dissipation in %s cannot be negative", path); 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 108 PetscFunctionReturn(PETSC_SUCCESS); 109 } 110 111 /* 112 * @brief Read the STGRand 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 as the only two entries, separated by a single space. 115 * Assumes there are 7 columns in the file. 116 * 117 * @param[in] comm MPI_Comm for the program 118 * @param[in] path Path to the STGRand.dat file 119 * @param[in,out] stg_ctx STGShur14Context where the data will be loaded into 120 */ 121 static PetscErrorCode ReadSTGRand(const MPI_Comm comm, const char path[PETSC_MAX_PATH_LEN], STGShur14Context stg_ctx) { 122 PetscInt dims[2]; 123 int ndims; 124 FILE *fp; 125 const PetscInt char_array_len = 512; 126 char line[char_array_len]; 127 char **array; 128 129 PetscFunctionBeginUser; 130 PetscCall(PHASTADatFileOpen(comm, path, char_array_len, dims, &fp)); 131 132 CeedScalar *phi = &stg_ctx->data[stg_ctx->offsets.phi]; 133 CeedScalar(*d)[stg_ctx->nmodes] = (CeedScalar(*)[stg_ctx->nmodes]) & stg_ctx->data[stg_ctx->offsets.d]; 134 CeedScalar(*sigma)[stg_ctx->nmodes] = (CeedScalar(*)[stg_ctx->nmodes]) & stg_ctx->data[stg_ctx->offsets.sigma]; 135 136 for (PetscInt i = 0; i < stg_ctx->nmodes; i++) { 137 PetscCall(PetscSynchronizedFGets(comm, fp, char_array_len, line)); 138 PetscCall(PetscStrToArray(line, ' ', &ndims, &array)); 139 PetscCheck(ndims == dims[1], comm, PETSC_ERR_FILE_UNEXPECTED, 140 "Line %" PetscInt_FMT " of %s does not have correct number of columns (%d instead of %" PetscInt_FMT ")", i, path, ndims, dims[1]); 141 142 d[0][i] = (CeedScalar)atof(array[0]); 143 d[1][i] = (CeedScalar)atof(array[1]); 144 d[2][i] = (CeedScalar)atof(array[2]); 145 phi[i] = (CeedScalar)atof(array[3]); 146 sigma[0][i] = (CeedScalar)atof(array[4]); 147 sigma[1][i] = (CeedScalar)atof(array[5]); 148 sigma[2][i] = (CeedScalar)atof(array[6]); 149 } 150 PetscCall(PetscFClose(comm, fp)); 151 152 PetscFunctionReturn(PETSC_SUCCESS); 153 } 154 155 /* 156 * @brief Read STG data from input paths and put in STGShur14Context 157 * 158 * Reads data from input paths and puts them into a STGShur14Context object. 159 * Data stored initially in `*stg_ctx` will be copied over to the new STGShur14Context instance. 160 * 161 * @param[in] comm MPI_Comm for the program 162 * @param[in] dm DM for the program 163 * @param[in] stg_inflow_path Path to STGInflow.dat file 164 * @param[in] stg_rand_path Path to STGRand.dat file 165 * @param[in,out] stg_ctx Pointer to STGShur14Context where the data will be loaded into 166 */ 167 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], 168 STGShur14Context *stg_ctx) { 169 PetscInt nmodes, nprofs; 170 PetscFunctionBeginUser; 171 172 // Get options 173 PetscCall(PHASTADatFileGetNRows(comm, stg_rand_path, &nmodes)); 174 PetscCall(PHASTADatFileGetNRows(comm, stg_inflow_path, &nprofs)); 175 PetscCheck(nmodes < STG_NMODES_MAX, comm, PETSC_ERR_SUP, 176 "Number of wavemodes in %s (%" PetscInt_FMT ") exceeds STG_NMODES_MAX (%d). Change size of STG_NMODES_MAX and recompile", stg_rand_path, 177 nmodes, STG_NMODES_MAX); 178 179 { 180 STGShur14Context temp_ctx; 181 PetscCall(PetscCalloc1(1, &temp_ctx)); 182 *temp_ctx = **stg_ctx; 183 temp_ctx->nmodes = nmodes; 184 temp_ctx->nprofs = nprofs; 185 temp_ctx->offsets.sigma = 0; 186 temp_ctx->offsets.d = nmodes * 3; 187 temp_ctx->offsets.phi = temp_ctx->offsets.d + nmodes * 3; 188 temp_ctx->offsets.kappa = temp_ctx->offsets.phi + nmodes; 189 temp_ctx->offsets.wall_dist = temp_ctx->offsets.kappa + nmodes; 190 temp_ctx->offsets.ubar = temp_ctx->offsets.wall_dist + nprofs; 191 temp_ctx->offsets.cij = temp_ctx->offsets.ubar + nprofs * 3; 192 temp_ctx->offsets.eps = temp_ctx->offsets.cij + nprofs * 6; 193 temp_ctx->offsets.lt = temp_ctx->offsets.eps + nprofs; 194 PetscInt total_num_scalars = temp_ctx->offsets.lt + nprofs; 195 temp_ctx->total_bytes = sizeof(*temp_ctx) + total_num_scalars * sizeof(temp_ctx->data[0]); 196 PetscCall(PetscFree(*stg_ctx)); 197 PetscCall(PetscMalloc(temp_ctx->total_bytes, stg_ctx)); 198 **stg_ctx = *temp_ctx; 199 PetscCall(PetscFree(temp_ctx)); 200 } 201 202 PetscCall(ReadSTGInflow(comm, stg_inflow_path, *stg_ctx)); 203 PetscCall(ReadSTGRand(comm, stg_rand_path, *stg_ctx)); 204 205 { // -- Calculate kappa 206 CeedScalar *kappa = &(*stg_ctx)->data[(*stg_ctx)->offsets.kappa]; 207 CeedScalar *wall_dist = &(*stg_ctx)->data[(*stg_ctx)->offsets.wall_dist]; 208 CeedScalar *lt = &(*stg_ctx)->data[(*stg_ctx)->offsets.lt]; 209 CeedScalar le, le_max = 0; 210 211 CeedPragmaSIMD for (PetscInt i = 0; i < (*stg_ctx)->nprofs; i++) { 212 le = PetscMin(2 * wall_dist[i], 3 * lt[i]); 213 if (le_max < le) le_max = le; 214 } 215 CeedScalar kmin = M_PI / le_max; 216 217 CeedPragmaSIMD for (PetscInt i = 0; i < (*stg_ctx)->nmodes; i++) { kappa[i] = kmin * pow((*stg_ctx)->alpha, i); } 218 } 219 220 PetscFunctionReturn(PETSC_SUCCESS); 221 } 222 223 PetscErrorCode SetupSTG(const MPI_Comm comm, const DM dm, ProblemData *problem, User user, const bool prescribe_T, const CeedScalar theta0, 224 const CeedScalar P0) { 225 Ceed ceed = user->ceed; 226 char stg_inflow_path[PETSC_MAX_PATH_LEN] = "./STGInflow.dat"; 227 char stg_rand_path[PETSC_MAX_PATH_LEN] = "./STGRand.dat"; 228 PetscBool mean_only = PETSC_FALSE, use_stgstrong = PETSC_FALSE, use_fluctuating_IC = PETSC_FALSE, given_stg_dx = PETSC_FALSE; 229 CeedScalar u0 = 0.0, alpha = 1.01, stg_dx = 1.0e-3; 230 CeedQFunctionContext stg_context; 231 NewtonianIdealGasContext newtonian_ig_ctx; 232 PetscFunctionBeginUser; 233 234 // Get options 235 PetscOptionsBegin(comm, NULL, "STG Boundary Condition Options", NULL); 236 PetscCall(PetscOptionsString("-stg_inflow_path", "Path to STGInflow.dat", NULL, stg_inflow_path, stg_inflow_path, sizeof(stg_inflow_path), NULL)); 237 PetscCall(PetscOptionsString("-stg_rand_path", "Path to STGInflow.dat", NULL, stg_rand_path, stg_rand_path, sizeof(stg_rand_path), NULL)); 238 PetscCall(PetscOptionsReal("-stg_alpha", "Growth rate of the wavemodes", NULL, alpha, &alpha, NULL)); 239 PetscCall(PetscOptionsReal("-stg_u0", "Advective velocity for the fluctuations", NULL, u0, &u0, NULL)); 240 PetscCall(PetscOptionsBool("-stg_mean_only", "Only apply mean profile", NULL, mean_only, &mean_only, NULL)); 241 PetscCall(PetscOptionsBool("-stg_strong", "Enforce STG inflow strongly", NULL, use_stgstrong, &use_stgstrong, NULL)); 242 PetscCall(PetscOptionsBool("-stg_fluctuating_IC", "\"Extrude\" the fluctuations through the domain as an initial condition", NULL, 243 use_fluctuating_IC, &use_fluctuating_IC, NULL)); 244 PetscCall(PetscOptionsReal("-stg_dx", "Element size in streamwise direction at inflow", NULL, stg_dx, &stg_dx, &given_stg_dx)); 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 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 } 267 268 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfunction_context, CEED_MEM_HOST, &newtonian_ig_ctx)); 269 global_stg_ctx->newtonian_ctx = *newtonian_ig_ctx; 270 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfunction_context, &newtonian_ig_ctx)); 271 272 PetscCall(GetSTGContextData(comm, dm, stg_inflow_path, stg_rand_path, &global_stg_ctx)); 273 274 PetscCallCeed(ceed, CeedQFunctionContextCreate(user->ceed, &stg_context)); 275 PetscCallCeed(ceed, CeedQFunctionContextSetData(stg_context, CEED_MEM_HOST, CEED_USE_POINTER, global_stg_ctx->total_bytes, global_stg_ctx)); 276 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(stg_context, CEED_MEM_HOST, FreeContextPetsc)); 277 PetscCallCeed(ceed, CeedQFunctionContextRegisterDouble(stg_context, "solution time", offsetof(struct STGShur14Context_, time), 1, 278 "Physical time of the solution")); 279 280 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&problem->ics.qfunction_context)); 281 problem->ics.qfunction = ICsSTG; 282 problem->ics.qfunction_loc = ICsSTG_loc; 283 problem->ics.qfunction_context = stg_context; 284 285 if (use_stgstrong) { 286 // Use default boundary integral QF (BoundaryIntegral) in newtonian.h 287 problem->use_strong_bc_ceed = PETSC_TRUE; 288 problem->bc_from_ics = PETSC_FALSE; 289 } else { 290 problem->apply_inflow.qfunction = STGShur14_Inflow; 291 problem->apply_inflow.qfunction_loc = STGShur14_Inflow_loc; 292 problem->apply_inflow_jacobian.qfunction = STGShur14_Inflow_Jacobian; 293 problem->apply_inflow_jacobian.qfunction_loc = STGShur14_Inflow_Jacobian_loc; 294 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(stg_context, &problem->apply_inflow.qfunction_context)); 295 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(stg_context, &problem->apply_inflow_jacobian.qfunction_context)); 296 problem->bc_from_ics = PETSC_TRUE; 297 } 298 299 PetscFunctionReturn(PETSC_SUCCESS); 300 } 301 302 // @brief Set STG strongly enforce components using DMAddBoundary 303 PetscErrorCode SetupStrongSTG(DM dm, SimpleBC bc, ProblemData *problem, Physics phys) { 304 DMLabel label; 305 PetscFunctionBeginUser; 306 307 PetscInt comps[5], num_comps = 4; 308 switch (phys->state_var) { 309 case STATEVAR_CONSERVATIVE: 310 // {0,1,2,3} for rho, rho*u, rho*v, rho*w 311 for (int i = 0; i < 4; i++) comps[i] = i; 312 break; 313 314 case STATEVAR_PRIMITIVE: 315 // {1,2,3,4} for u, v, w, T 316 for (int i = 0; i < 4; i++) comps[i] = i + 1; 317 break; 318 } 319 320 PetscCall(DMGetLabel(dm, "Face Sets", &label)); 321 if (bc->num_inflow > 0) { 322 PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "STG", label, bc->num_inflow, bc->inflows, 0, num_comps, comps, NULL, NULL, NULL, NULL)); 323 } 324 325 PetscFunctionReturn(PETSC_SUCCESS); 326 } 327 328 PetscErrorCode SetupStrongSTG_QF(Ceed ceed, ProblemData *problem, CeedInt num_comp_x, CeedInt num_comp_q, CeedInt stg_data_size, 329 CeedInt q_data_size_sur, CeedQFunction *qf_strongbc) { 330 PetscFunctionBeginUser; 331 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, STGShur14_Inflow_StrongQF, STGShur14_Inflow_StrongQF_loc, qf_strongbc)); 332 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_strongbc, "surface qdata", q_data_size_sur, CEED_EVAL_NONE)); 333 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_strongbc, "x", num_comp_x, CEED_EVAL_NONE)); 334 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_strongbc, "scale", 1, CEED_EVAL_NONE)); 335 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_strongbc, "stg data", stg_data_size, CEED_EVAL_NONE)); 336 PetscCallCeed(ceed, CeedQFunctionAddOutput(*qf_strongbc, "q", num_comp_q, CEED_EVAL_NONE)); 337 338 PetscCallCeed(ceed, CeedQFunctionSetContext(*qf_strongbc, problem->ics.qfunction_context)); 339 PetscFunctionReturn(PETSC_SUCCESS); 340 } 341 342 PetscErrorCode SetupStrongSTG_PreProcessing(Ceed ceed, ProblemData *problem, CeedInt num_comp_x, CeedInt stg_data_size, CeedInt q_data_size_sur, 343 CeedQFunction *qf_strongbc) { 344 PetscFunctionBeginUser; 345 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, Preprocess_STGShur14, Preprocess_STGShur14_loc, qf_strongbc)); 346 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_strongbc, "surface qdata", q_data_size_sur, CEED_EVAL_NONE)); 347 PetscCallCeed(ceed, CeedQFunctionAddInput(*qf_strongbc, "x", num_comp_x, CEED_EVAL_NONE)); 348 PetscCallCeed(ceed, CeedQFunctionAddOutput(*qf_strongbc, "stg data", stg_data_size, CEED_EVAL_NONE)); 349 350 PetscCallCeed(ceed, CeedQFunctionSetContext(*qf_strongbc, problem->ics.qfunction_context)); 351 PetscFunctionReturn(PETSC_SUCCESS); 352 } 353