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 #ifndef M_PI 20 #define M_PI 3.14159265358979323846 21 #endif 22 23 /* 24 * @brief Perform Cholesky decomposition on array of symmetric 3x3 matrices 25 * 26 * This assumes the input matrices are in order [11,22,33,12,13,23]. This 27 * format is also used for the output. 28 * 29 * @param[in] comm MPI_Comm 30 * @param[in] nprofs Number of matrices in Rij 31 * @param[in] Rij Array of the symmetric matrices [6,nprofs] 32 * @param[out] Cij Array of the Cholesky Decomposition matrices, [6,nprofs] 33 */ 34 PetscErrorCode CalcCholeskyDecomp(MPI_Comm comm, PetscInt nprofs, 35 const CeedScalar Rij[6][nprofs], CeedScalar Cij[6][nprofs]) { 36 37 PetscFunctionBeginUser; 38 for (PetscInt i=0; i<nprofs; i++) { 39 Cij[0][i] = sqrt(Rij[0][i]); 40 Cij[3][i] = Rij[3][i] / Cij[0][i]; 41 Cij[1][i] = sqrt(Rij[1][i] - pow(Cij[3][i], 2) ); 42 Cij[4][i] = Rij[4][i] / Cij[0][i]; 43 Cij[5][i] = (Rij[5][i] - Cij[3][i]*Cij[4][i]) / Cij[1][i]; 44 Cij[2][i] = sqrt(Rij[2][i] - pow(Cij[4][i], 2) - pow(Cij[5][i], 2)); 45 46 if (isnan(Cij[0][i]) || isnan(Cij[1][i]) || isnan(Cij[2][i])) 47 SETERRQ(comm, -1, "Cholesky decomposition failed at profile point %d. " 48 "Either STGInflow has non-SPD matrix or contains nan.", i+1); 49 } 50 PetscFunctionReturn(0); 51 } 52 53 54 /* 55 * @brief Open a PHASTA *.dat file, grabbing dimensions and file pointer 56 * 57 * This function opens the file specified by `path` using `PetscFOpen` and 58 * passes the file pointer in `fp`. It is not closed in this function, thus 59 * `fp` must be closed sometime after this function has been called (using 60 * `PetscFClose` for example). 61 * 62 * Assumes that the first line of the file has the number of rows and columns 63 * as the only two entries, separated by a single space 64 * 65 * @param[in] comm MPI_Comm for the program 66 * @param[in] path Path to the file 67 * @param[in] char_array_len Length of the character array that should contain each line 68 * @param[out] dims Dimensions of the file, taken from the first line of the file 69 * @param[out] fp File pointer to the opened file 70 */ 71 static PetscErrorCode OpenPHASTADatFile(const MPI_Comm comm, 72 const char path[PETSC_MAX_PATH_LEN], const PetscInt char_array_len, 73 PetscInt dims[2], FILE **fp) { 74 PetscErrorCode ierr; 75 PetscInt ndims; 76 char line[char_array_len]; 77 char **array; 78 79 PetscFunctionBeginUser; 80 ierr = PetscFOpen(comm, path, "r", fp); CHKERRQ(ierr); 81 ierr = PetscSynchronizedFGets(comm, *fp, char_array_len, line); CHKERRQ(ierr); 82 ierr = PetscStrToArray(line, ' ', &ndims, &array); CHKERRQ(ierr); 83 if (ndims != 2) SETERRQ(comm, -1, 84 "Found %d dimensions instead of 2 on the first line of %s", 85 ndims, path); 86 87 for (PetscInt i=0; i<ndims; i++) dims[i] = atoi(array[i]); 88 ierr = PetscStrToArrayDestroy(ndims, array); CHKERRQ(ierr); 89 PetscFunctionReturn(0); 90 } 91 92 /* 93 * @brief Get the number of rows for the PHASTA file at path 94 * 95 * Assumes that the first line of the file has the number of rows and columns 96 * as the only two entries, separated by a single space 97 * 98 * @param[in] comm MPI_Comm for the program 99 * @param[in] path Path to the file 100 * @param[out] nrows Number of rows 101 */ 102 static PetscErrorCode GetNRows(const MPI_Comm comm, 103 const char path[PETSC_MAX_PATH_LEN], PetscInt *nrows) { 104 PetscErrorCode ierr; 105 const PetscInt char_array_len = 512; 106 PetscInt dims[2]; 107 FILE *fp; 108 109 PetscFunctionBeginUser; 110 ierr = OpenPHASTADatFile(comm, path, char_array_len, dims, &fp); CHKERRQ(ierr); 111 *nrows = dims[0]; 112 ierr = PetscFClose(comm, fp); CHKERRQ(ierr); 113 PetscFunctionReturn(0); 114 } 115 116 /* 117 * @brief Read the STGInflow file and load the contents into stg_ctx 118 * 119 * Assumes that the first line of the file has the number of rows and columns 120 * as the only two entries, separated by a single space. 121 * Assumes there are 14 columns in the file 122 * 123 * Function calculates the Cholesky decomposition from the Reynolds stress 124 * profile found in the file 125 * 126 * @param[in] comm MPI_Comm for the program 127 * @param[in] path Path to the STGInflow.dat file 128 * @param[inout] stg_ctx STGShur14Context where the data will be loaded into 129 */ 130 static PetscErrorCode ReadSTGInflow(const MPI_Comm comm, 131 const char path[PETSC_MAX_PATH_LEN], STGShur14Context stg_ctx) { 132 PetscErrorCode ierr; 133 PetscInt ndims, dims[2]; 134 FILE *fp; 135 const PetscInt char_array_len=512; 136 char line[char_array_len]; 137 char **array; 138 139 PetscFunctionBeginUser; 140 141 ierr = OpenPHASTADatFile(comm, path, char_array_len, dims, &fp); CHKERRQ(ierr); 142 143 CeedScalar rij[6][stg_ctx->nprofs]; 144 CeedScalar *prof_dw = &stg_ctx->data[stg_ctx->offsets.prof_dw]; 145 CeedScalar *eps = &stg_ctx->data[stg_ctx->offsets.eps]; 146 CeedScalar *lt = &stg_ctx->data[stg_ctx->offsets.lt]; 147 CeedScalar (*ubar)[stg_ctx->nprofs] = (CeedScalar (*)[stg_ctx->nprofs]) 148 &stg_ctx->data[stg_ctx->offsets.ubar]; 149 150 for (PetscInt i=0; i<stg_ctx->nprofs; i++) { 151 ierr = PetscSynchronizedFGets(comm, fp, char_array_len, line); CHKERRQ(ierr); 152 ierr = PetscStrToArray(line, ' ', &ndims, &array); CHKERRQ(ierr); 153 if (ndims < dims[1]) SETERRQ(comm, -1, 154 "Line %d of %s does not contain enough columns (%d instead of %d)", i, 155 path, ndims, dims[1]); 156 157 prof_dw[i] = (CeedScalar) atof(array[0]); 158 ubar[0][i] = (CeedScalar) atof(array[1]); 159 ubar[1][i] = (CeedScalar) atof(array[2]); 160 ubar[2][i] = (CeedScalar) atof(array[3]); 161 rij[0][i] = (CeedScalar) atof(array[4]); 162 rij[1][i] = (CeedScalar) atof(array[5]); 163 rij[2][i] = (CeedScalar) atof(array[6]); 164 rij[3][i] = (CeedScalar) atof(array[7]); 165 rij[4][i] = (CeedScalar) atof(array[8]); 166 rij[5][i] = (CeedScalar) atof(array[9]); 167 lt[i] = (CeedScalar) atof(array[12]); 168 eps[i] = (CeedScalar) atof(array[13]); 169 170 if (prof_dw[i] < 0) SETERRQ(comm, -1, 171 "Distance to wall in %s cannot be negative", path); 172 if (lt[i] < 0) SETERRQ(comm, -1, 173 "Turbulent length scale in %s cannot be negative", path); 174 if (eps[i] < 0) SETERRQ(comm, -1, 175 "Turbulent dissipation in %s cannot be negative", path); 176 177 } 178 CeedScalar (*cij)[stg_ctx->nprofs] = (CeedScalar (*)[stg_ctx->nprofs]) 179 &stg_ctx->data[stg_ctx->offsets.cij]; 180 ierr = CalcCholeskyDecomp(comm, stg_ctx->nprofs, rij, cij); CHKERRQ(ierr); 181 ierr = PetscFClose(comm, fp); CHKERRQ(ierr); 182 PetscFunctionReturn(0); 183 } 184 185 /* 186 * @brief Read the STGRand file and load the contents into stg_ctx 187 * 188 * Assumes that the first line of the file has the number of rows and columns 189 * as the only two entries, separated by a single space. 190 * Assumes there are 7 columns in the file 191 * 192 * @param[in] comm MPI_Comm for the program 193 * @param[in] path Path to the STGRand.dat file 194 * @param[inout] stg_ctx STGShur14Context where the data will be loaded into 195 */ 196 static PetscErrorCode ReadSTGRand(const MPI_Comm comm, 197 const char path[PETSC_MAX_PATH_LEN], 198 STGShur14Context stg_ctx) { 199 200 PetscErrorCode ierr; 201 PetscInt ndims, dims[2]; 202 FILE *fp; 203 const PetscInt char_array_len = 512; 204 char line[char_array_len]; 205 char **array; 206 207 PetscFunctionBeginUser; 208 ierr = OpenPHASTADatFile(comm, path, char_array_len, dims, &fp); CHKERRQ(ierr); 209 210 CeedScalar *phi = &stg_ctx->data[stg_ctx->offsets.phi]; 211 CeedScalar (*d)[stg_ctx->nmodes] = (CeedScalar (*)[stg_ctx->nmodes]) 212 &stg_ctx->data[stg_ctx->offsets.d]; 213 CeedScalar (*sigma)[stg_ctx->nmodes] = (CeedScalar (*)[stg_ctx->nmodes]) 214 &stg_ctx->data[stg_ctx->offsets.sigma]; 215 216 for (PetscInt i=0; i<stg_ctx->nmodes; i++) { 217 ierr = PetscSynchronizedFGets(comm, fp, char_array_len, line); CHKERRQ(ierr); 218 ierr = PetscStrToArray(line, ' ', &ndims, &array); CHKERRQ(ierr); 219 if (ndims < dims[1]) SETERRQ(comm, -1, 220 "Line %d of %s does not contain enough columns (%d instead of %d)", i, 221 path, ndims, dims[1]); 222 223 d[0][i] = (CeedScalar) atof(array[0]); 224 d[1][i] = (CeedScalar) atof(array[1]); 225 d[2][i] = (CeedScalar) atof(array[2]); 226 phi[i] = (CeedScalar) atof(array[3]); 227 sigma[0][i] = (CeedScalar) atof(array[4]); 228 sigma[1][i] = (CeedScalar) atof(array[5]); 229 sigma[2][i] = (CeedScalar) atof(array[6]); 230 } 231 ierr = PetscFClose(comm, fp); CHKERRQ(ierr); 232 PetscFunctionReturn(0); 233 } 234 235 /* 236 * @brief Read STG data from input paths and put in STGShur14Context 237 * 238 * Reads data from input paths and puts them into a STGShur14Context object. 239 * Data stored initially in `*pstg_ctx` will be copied over to the new 240 * STGShur14Context instance. 241 * 242 * @param[in] comm MPI_Comm for the program 243 * @param[in] dm DM for the program 244 * @param[in] stg_inflow_path Path to STGInflow.dat file 245 * @param[in] stg_rand_path Path to STGRand.dat file 246 * @param[inout] pstg_ctx Pointer to STGShur14Context where the data will be loaded into 247 */ 248 PetscErrorCode GetSTGContextData(const MPI_Comm comm, const DM dm, 249 char stg_inflow_path[PETSC_MAX_PATH_LEN], 250 char stg_rand_path[PETSC_MAX_PATH_LEN], 251 STGShur14Context *pstg_ctx) { 252 PetscErrorCode ierr; 253 PetscInt nmodes, nprofs; 254 STGShur14Context stg_ctx; 255 PetscFunctionBeginUser; 256 257 // Get options 258 ierr = GetNRows(comm, stg_rand_path, &nmodes); CHKERRQ(ierr); 259 ierr = GetNRows(comm, stg_inflow_path, &nprofs); CHKERRQ(ierr); 260 if (nmodes > STG_NMODES_MAX) 261 SETERRQ(comm, 1, "Number of wavemodes in %s (%d) exceeds STG_NMODES_MAX (%d). " 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.prof_dw = s->offsets.kappa + nmodes; 276 s->offsets.ubar = s->offsets.prof_dw + 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 PetscInt total_num_scalars = s->offsets.lt + nprofs; 281 s->total_bytes = sizeof(*stg_ctx) + total_num_scalars*sizeof(stg_ctx->data[0]); 282 ierr = PetscMalloc(s->total_bytes, &stg_ctx); CHKERRQ(ierr); 283 *stg_ctx = *s; 284 ierr = PetscFree(s); CHKERRQ(ierr); 285 } 286 287 ierr = ReadSTGInflow(comm, stg_inflow_path, stg_ctx); CHKERRQ(ierr); 288 ierr = ReadSTGRand(comm, stg_rand_path, stg_ctx); CHKERRQ(ierr); 289 290 // -- Calculate kappa 291 { 292 CeedScalar *kappa = &stg_ctx->data[stg_ctx->offsets.kappa]; 293 CeedScalar *prof_dw = &stg_ctx->data[stg_ctx->offsets.prof_dw]; 294 CeedScalar *lt = &stg_ctx->data[stg_ctx->offsets.lt]; 295 CeedScalar le, le_max=0; 296 297 CeedPragmaSIMD 298 for (PetscInt i=0; i<stg_ctx->nprofs; i++) { 299 le = PetscMin(2*prof_dw[i], 3*lt[i]); 300 if (le_max < le) le_max = le; 301 } 302 CeedScalar kmin = M_PI/le_max; 303 304 CeedPragmaSIMD 305 for (PetscInt i=0; i<stg_ctx->nmodes; i++) { 306 kappa[i] = kmin*pow(stg_ctx->alpha, i); 307 } 308 } //end calculate kappa 309 310 *pstg_ctx = stg_ctx; 311 PetscFunctionReturn(0); 312 } 313 314 PetscErrorCode SetupSTG(const MPI_Comm comm, const DM dm, ProblemData *problem, 315 User user, const bool prescribe_T, 316 const CeedScalar theta0, const CeedScalar P0) { 317 PetscErrorCode ierr; 318 char stg_inflow_path[PETSC_MAX_PATH_LEN] = "./STGInflow.dat"; 319 char stg_rand_path[PETSC_MAX_PATH_LEN] = "./STGRand.dat"; 320 PetscBool mean_only = PETSC_FALSE; 321 CeedScalar u0=0.0, alpha=1.01; 322 STGShur14Context stg_ctx; 323 CeedQFunctionContext stg_context; 324 NewtonianIdealGasContext newtonian_ig_ctx; 325 PetscFunctionBeginUser; 326 327 // Get options 328 PetscOptionsBegin(comm, NULL, "STG Boundary Condition Options", NULL); 329 ierr = PetscOptionsString("-stg_inflow_path", "Path to STGInflow.dat", NULL, 330 stg_inflow_path, stg_inflow_path, 331 sizeof(stg_inflow_path), NULL); CHKERRQ(ierr); 332 ierr = PetscOptionsString("-stg_rand_path", "Path to STGInflow.dat", NULL, 333 stg_rand_path,stg_rand_path, 334 sizeof(stg_rand_path), NULL); CHKERRQ(ierr); 335 ierr = PetscOptionsReal("-stg_alpha", "Growth rate of the wavemodes", NULL, 336 alpha, &alpha, NULL); CHKERRQ(ierr); 337 ierr = PetscOptionsReal("-stg_u0", "Advective velocity for the fluctuations", 338 NULL, u0, &u0, NULL); CHKERRQ(ierr); 339 ierr = PetscOptionsBool("-stg_mean_only", "Only apply mean profile", 340 NULL, mean_only, &mean_only, NULL); CHKERRQ(ierr); 341 PetscOptionsEnd(); 342 343 ierr = PetscCalloc1(1, &stg_ctx); CHKERRQ(ierr); 344 stg_ctx->alpha = alpha; 345 stg_ctx->u0 = u0; 346 stg_ctx->is_implicit = user->phys->implicit; 347 stg_ctx->prescribe_T = prescribe_T; 348 stg_ctx->mean_only = mean_only; 349 stg_ctx->theta0 = theta0; 350 stg_ctx->P0 = P0; 351 352 { 353 // Calculate dx assuming constant spacing 354 PetscReal domain_min[3], domain_max[3], domain_size[3]; 355 ierr = DMGetBoundingBox(dm, domain_min, domain_max); CHKERRQ(ierr); 356 for (PetscInt i=0; i<3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 357 358 PetscInt nmax = 3, faces[3]; 359 ierr = PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &nmax, 360 NULL); CHKERRQ(ierr); 361 stg_ctx->dx = domain_size[0]/faces[0]; 362 } 363 364 CeedQFunctionContextGetData(problem->apply_vol_rhs.qfunction_context, 365 CEED_MEM_HOST, &newtonian_ig_ctx); 366 stg_ctx->newtonian_ctx = *newtonian_ig_ctx; 367 CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfunction_context, 368 &newtonian_ig_ctx); 369 370 ierr = GetSTGContextData(comm, dm, stg_inflow_path, stg_rand_path, &stg_ctx); 371 CHKERRQ(ierr); 372 373 CeedQFunctionContextDestroy(&problem->apply_inflow.qfunction_context); 374 CeedQFunctionContextCreate(user->ceed, &stg_context); 375 CeedQFunctionContextSetData(stg_context, CEED_MEM_HOST, 376 CEED_USE_POINTER, stg_ctx->total_bytes, stg_ctx); 377 CeedQFunctionContextSetDataDestroy(stg_context, CEED_MEM_HOST, 378 FreeContextPetsc); 379 CeedQFunctionContextRegisterDouble(stg_context, "solution time", 380 offsetof(struct STGShur14Context_, time), 1, 381 "Phyiscal time of the solution"); 382 383 problem->apply_inflow.qfunction = STGShur14_Inflow; 384 problem->apply_inflow.qfunction_loc = STGShur14_Inflow_loc; 385 problem->apply_inflow.qfunction_context = stg_context; 386 387 PetscFunctionReturn(0); 388 } 389