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