1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors. 2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause 3 4 /// @file 5 /// Utility functions for setting up Blasius Boundary Layer 6 7 #include "../qfunctions/blasius.h" 8 9 #include <ceed.h> 10 #include <petscdm.h> 11 #include <petscdt.h> 12 13 #include <differential_filter.h> 14 #include <navierstokes.h> 15 #include "stg_shur14.h" 16 17 PetscErrorCode CompressibleBlasiusResidual(SNES snes, Vec X, Vec R, void *ctx) { 18 const BlasiusContext blasius = (BlasiusContext)ctx; 19 const PetscScalar *Tf, *Th; // Chebyshev coefficients 20 PetscScalar *r, f[4], h[4]; 21 PetscInt N = blasius->n_cheb; 22 State S_infty = blasius->S_infty; 23 CeedScalar U_infty = Norm3(S_infty.Y.velocity); 24 NewtonianIGProperties gas = blasius->newt_ctx.gas; 25 26 PetscFunctionBeginUser; 27 PetscScalar Ma = Mach(gas, S_infty.Y.temperature, U_infty), Pr = Prandtl(gas), gamma = HeatCapacityRatio(gas); 28 29 PetscCall(VecGetArrayRead(X, &Tf)); 30 Th = Tf + N; 31 PetscCall(VecGetArray(R, &r)); 32 33 // Left boundary conditions f = f' = 0 34 ChebyshevEval(N, Tf, -1., blasius->eta_max, f); 35 r[0] = f[0]; 36 r[1] = f[1]; 37 38 // f - right end boundary condition 39 ChebyshevEval(N, Tf, 1., blasius->eta_max, f); 40 r[2] = f[1] - 1.; 41 42 for (int i = 0; i < N - 3; i++) { 43 ChebyshevEval(N, Tf, blasius->X[i], blasius->eta_max, f); 44 ChebyshevEval(N - 1, Th, blasius->X[i], blasius->eta_max, h); 45 // mu and rho generally depend on h. 46 // We naively assume constant mu. 47 // For an ideal gas at constant pressure, density is inversely proportional to enthalpy. 48 // The *_tilde values are *relative* to their freestream values, and we proved first derivatives here. 49 const PetscScalar mu_tilde[2] = {1, 0}; 50 const PetscScalar rho_tilde[2] = {1 / h[0], -h[1] / PetscSqr(h[0])}; 51 const PetscScalar mu_rho_tilde[2] = { 52 mu_tilde[0] * rho_tilde[0], 53 mu_tilde[1] * rho_tilde[0] + mu_tilde[0] * rho_tilde[1], 54 }; 55 r[3 + i] = 2 * (mu_rho_tilde[0] * f[3] + mu_rho_tilde[1] * f[2]) + f[2] * f[0]; 56 r[N + 2 + i] = (mu_rho_tilde[0] * h[2] + mu_rho_tilde[1] * h[1]) + Pr * f[0] * h[1] + Pr * (gamma - 1) * mu_rho_tilde[0] * PetscSqr(Ma * f[2]); 57 } 58 59 // h - left end boundary condition 60 ChebyshevEval(N - 1, Th, -1., blasius->eta_max, h); 61 r[N] = h[0] - blasius->T_wall / S_infty.Y.temperature; 62 63 // h - right end boundary condition 64 ChebyshevEval(N - 1, Th, 1., blasius->eta_max, h); 65 r[N + 1] = h[0] - 1.; 66 67 // Restore vectors 68 PetscCall(VecRestoreArrayRead(X, &Tf)); 69 PetscCall(VecRestoreArray(R, &r)); 70 PetscFunctionReturn(PETSC_SUCCESS); 71 } 72 73 PetscErrorCode ComputeChebyshevCoefficients(BlasiusContext blasius) { 74 SNES snes; 75 Vec sol, res; 76 PetscReal *w; 77 PetscInt N = blasius->n_cheb; 78 SNESConvergedReason reason; 79 const PetscScalar *cheb_coefs; 80 81 PetscFunctionBeginUser; 82 // Allocate memory 83 PetscCall(PetscMalloc2(N - 3, &blasius->X, N - 3, &w)); 84 PetscCall(PetscDTGaussQuadrature(N - 3, -1., 1., blasius->X, w)); 85 86 // Snes solve 87 PetscCall(SNESCreate(PETSC_COMM_SELF, &snes)); 88 PetscCall(VecCreate(PETSC_COMM_SELF, &sol)); 89 PetscCall(VecSetSizes(sol, PETSC_DECIDE, 2 * N - 1)); 90 PetscCall(VecSetFromOptions(sol)); 91 // Constant relative enthalpy 1 as initial guess 92 PetscCall(VecSetValue(sol, N, 1., INSERT_VALUES)); 93 PetscCall(VecDuplicate(sol, &res)); 94 PetscCall(SNESSetFunction(snes, res, CompressibleBlasiusResidual, blasius)); 95 PetscCall(SNESSetOptionsPrefix(snes, "chebyshev_")); 96 PetscCall(SNESSetFromOptions(snes)); 97 PetscCall(SNESSolve(snes, NULL, sol)); 98 PetscCall(SNESGetConvergedReason(snes, &reason)); 99 PetscCheck(reason >= 0, PETSC_COMM_WORLD, PETSC_ERR_CONV_FAILED, "The Chebyshev solve failed.\n"); 100 101 // Assign Chebyshev coefficients 102 PetscCall(VecGetArrayRead(sol, &cheb_coefs)); 103 for (int i = 0; i < N; i++) blasius->Tf_cheb[i] = cheb_coefs[i]; 104 for (int i = 0; i < N - 1; i++) blasius->Th_cheb[i] = cheb_coefs[i + N]; 105 106 // Destroy objects 107 PetscCall(PetscFree2(blasius->X, w)); 108 PetscCall(VecDestroy(&sol)); 109 PetscCall(VecDestroy(&res)); 110 PetscCall(SNESDestroy(&snes)); 111 PetscFunctionReturn(PETSC_SUCCESS); 112 } 113 114 static PetscErrorCode GetYNodeLocs(const MPI_Comm comm, const char path[PETSC_MAX_PATH_LEN], PetscReal **pynodes, PetscInt *nynodes) { 115 int ndims, dims[2] = {0.}; 116 FILE *fp; 117 const PetscInt char_array_len = 512; 118 char line[char_array_len]; 119 PetscReal *node_locs; 120 121 PetscFunctionBeginUser; 122 PetscCall(PetscFOpen(comm, path, "r", &fp)); 123 PetscCall(PetscSynchronizedFGets(comm, fp, char_array_len, line)); 124 { 125 char **array; 126 127 PetscCall(PetscStrToArray(line, ' ', &ndims, &array)); 128 for (PetscInt i = 0; i < ndims; i++) dims[i] = atoi(array[i]); 129 PetscCall(PetscStrToArrayDestroy(ndims, array)); 130 } 131 if (ndims < 2) dims[1] = 1; // Assume 1 column of data is not otherwise specified 132 *nynodes = dims[0]; 133 PetscCall(PetscMalloc1(*nynodes, &node_locs)); 134 135 for (PetscInt i = 0; i < dims[0]; i++) { 136 char **array; 137 138 PetscCall(PetscSynchronizedFGets(comm, fp, char_array_len, line)); 139 PetscCall(PetscStrToArray(line, ' ', &ndims, &array)); 140 PetscCheck(ndims == dims[1], comm, PETSC_ERR_FILE_UNEXPECTED, 141 "Line %" PetscInt_FMT " of %s does not contain correct number of columns (%d instead of %d)", i, path, ndims, dims[1]); 142 143 node_locs[i] = (PetscReal)atof(array[0]); 144 PetscCall(PetscStrToArrayDestroy(ndims, array)); 145 } 146 PetscCall(PetscFClose(comm, fp)); 147 *pynodes = node_locs; 148 PetscFunctionReturn(PETSC_SUCCESS); 149 } 150 151 /* \brief Modify the domain and mesh for blasius 152 * 153 * Modifies mesh such that `N` elements are within `refine_height` with a geometric growth ratio of `growth`. Excess elements are then distributed 154 * linearly in logspace to the top surface. 155 * 156 * The top surface is also angled downwards, so that it may be used as an outflow. 157 * It's angle is controlled by `top_angle` (in units of degrees). 158 * 159 * If `node_locs` is not NULL, then the nodes will be placed at `node_locs` locations. 160 * If it is NULL, then the modified coordinate values will be set in the array, along with `num_node_locs`. 161 */ 162 static PetscErrorCode ModifyMesh(MPI_Comm comm, DM dm, PetscReal growth, PetscInt N, PetscReal refine_height, PetscReal top_angle, 163 PetscReal *node_locs[], PetscInt *num_node_locs) { 164 PetscInt narr, ncoords, dim; 165 PetscReal domain_min[3], domain_max[3], domain_size[3]; 166 PetscScalar *arr_coords; 167 Vec vec_coords; 168 169 PetscFunctionBeginUser; 170 PetscCall(DMGetDimension(dm, &dim)); 171 PetscReal angle_coeff = tan(top_angle * (M_PI / 180)); 172 // Get domain boundary information 173 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max)); 174 for (PetscInt i = 0; i < 3; i++) domain_size[i] = domain_max[i] - domain_min[i]; 175 176 // Get coords array from DM 177 PetscCall(DMGetCoordinatesLocal(dm, &vec_coords)); 178 PetscCall(VecGetLocalSize(vec_coords, &narr)); 179 PetscCall(VecGetArray(vec_coords, &arr_coords)); 180 181 PetscScalar(*coords)[dim] = (PetscScalar(*)[dim])arr_coords; 182 ncoords = narr / dim; 183 184 // Get mesh information 185 PetscInt nmax = 3, faces[3]; 186 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &nmax, NULL)); 187 // Get element size of the box mesh, for indexing each node 188 const PetscReal dybox = domain_size[1] / faces[1]; 189 190 if (!*node_locs) { 191 // Calculate the first element height 192 PetscReal dy1 = refine_height * (growth - 1) / (pow(growth, N) - 1); 193 194 // Calculate log of sizing outside BL 195 PetscReal logdy = (log(domain_max[1]) - log(refine_height)) / (faces[1] - N); 196 197 *num_node_locs = faces[1] + 1; 198 PetscReal *temp_node_locs; 199 PetscCall(PetscMalloc1(*num_node_locs, &temp_node_locs)); 200 201 for (PetscInt i = 0; i < ncoords; i++) { 202 PetscInt y_box_index = round(coords[i][1] / dybox); 203 if (y_box_index <= N) { 204 coords[i][1] = 205 (1 - (coords[i][0] - domain_min[0]) * angle_coeff / domain_max[1]) * dy1 * (pow(growth, coords[i][1] / dybox) - 1) / (growth - 1); 206 } else { 207 PetscInt j = y_box_index - N; 208 coords[i][1] = (1 - (coords[i][0] - domain_min[0]) * angle_coeff / domain_max[1]) * exp(log(refine_height) + logdy * j); 209 } 210 if (coords[i][0] == domain_min[0] && coords[i][2] == domain_min[2]) temp_node_locs[y_box_index] = coords[i][1]; 211 } 212 213 *node_locs = temp_node_locs; 214 } else { 215 PetscCheck(*num_node_locs >= faces[1] + 1, comm, PETSC_ERR_FILE_UNEXPECTED, 216 "The y_node_locs_path has too few locations; There are %" PetscInt_FMT " + 1 nodes, but only %" PetscInt_FMT " locations given", 217 faces[1] + 1, *num_node_locs); 218 if (*num_node_locs > faces[1] + 1) { 219 PetscCall(PetscPrintf(comm, 220 "WARNING: y_node_locs_path has more locations (%" PetscInt_FMT ") " 221 "than the mesh has nodes (%" PetscInt_FMT "). This maybe unintended.\n", 222 *num_node_locs, faces[1] + 1)); 223 } 224 PetscScalar max_y = (*node_locs)[faces[1]]; 225 226 for (PetscInt i = 0; i < ncoords; i++) { 227 // Determine which y-node we're at 228 PetscInt y_box_index = round(coords[i][1] / dybox); 229 coords[i][1] = (1 - (coords[i][0] - domain_min[0]) * angle_coeff / max_y) * (*node_locs)[y_box_index]; 230 } 231 } 232 233 PetscCall(VecRestoreArray(vec_coords, &arr_coords)); 234 PetscCall(DMSetCoordinatesLocal(dm, vec_coords)); 235 PetscFunctionReturn(PETSC_SUCCESS); 236 } 237 238 static PetscErrorCode BlasiusInflowBCSetup_CreateIFunctionQF(BCDefinition bc_def, CeedQFunction *qf) { 239 HoneeBCStruct honee_bc; 240 241 PetscFunctionBeginUser; 242 PetscCall(BCDefinitionGetContext(bc_def, &honee_bc)); 243 PetscCall(HoneeBCCreateIFunctionQF(bc_def, Blasius_Inflow, Blasius_Inflow_loc, honee_bc->qfctx, qf)); 244 PetscFunctionReturn(PETSC_SUCCESS); 245 } 246 247 static PetscErrorCode BlasiusInflowBCSetup_CreateIJacobianQF(BCDefinition bc_def, CeedQFunction *qf) { 248 HoneeBCStruct honee_bc; 249 250 PetscFunctionBeginUser; 251 PetscCall(BCDefinitionGetContext(bc_def, &honee_bc)); 252 PetscCall(HoneeBCCreateIJacobianQF(bc_def, Blasius_Inflow_Jacobian, Blasius_Inflow_Jacobian_loc, honee_bc->qfctx, qf)); 253 PetscFunctionReturn(PETSC_SUCCESS); 254 } 255 256 PetscErrorCode NS_BLASIUS(ProblemData problem, DM dm, void *ctx) { 257 Honee honee = *(Honee *)ctx; 258 MPI_Comm comm = honee->comm; 259 Ceed ceed = honee->ceed; 260 PetscBool use_stg = PETSC_FALSE; 261 BlasiusContext blasius_ctx; 262 NewtonianIdealGasContext newtonian_ig_ctx; 263 CeedQFunctionContext blasius_qfctx; 264 265 PetscFunctionBeginUser; 266 PetscCall(NS_NEWTONIAN_IG(problem, dm, ctx)); 267 PetscCall(PetscNew(&blasius_ctx)); 268 269 // ------------------------------------------------------ 270 // SET UP Blasius 271 // ------------------------------------------------------ 272 problem->ics = (HoneeQFSpec){.qf_func_ptr = ICsBlasius, .qf_loc = ICsBlasius_loc, .qfctx = problem->ics.qfctx}; 273 274 CeedScalar U_inf = 40; // m/s 275 CeedScalar T_inf = 288.; // K 276 CeedScalar T_wall = 288.; // K 277 CeedScalar delta0 = 4.2e-3; // m 278 CeedScalar P_inf = 1.01e5; // Pa 279 PetscInt N = 20; // Number of Chebyshev terms 280 PetscBool weakT = PETSC_FALSE; // weak density or temperature 281 PetscReal mesh_refine_height = 5.9e-4; // m 282 PetscReal mesh_growth = 1.08; // [-] 283 PetscInt mesh_Ndelta = 45; // [-] 284 PetscReal mesh_top_angle = 5; // degrees 285 char mesh_ynodes_path[PETSC_MAX_PATH_LEN] = ""; 286 PetscBool P0_set; 287 288 PetscOptionsBegin(comm, NULL, "Options for BLASIUS problem", NULL); 289 PetscCall(PetscOptionsBool("-weakT", "Change from rho weak to T weak at inflow", NULL, weakT, &weakT, NULL)); 290 PetscCall(PetscOptionsScalar("-velocity_infinity", "Velocity at boundary layer edge", NULL, U_inf, &U_inf, NULL)); 291 PetscCall(PetscOptionsScalar("-temperature_infinity", "Temperature at boundary layer edge", NULL, T_inf, &T_inf, NULL)); 292 PetscCall(PetscOptionsHasName(NULL, NULL, "-P0", &P0_set)); // For maintaining behavior of -P0 flag (which is deprecated) 293 PetscCall( 294 PetscOptionsDeprecated("-P0", "-pressure_infinity", "libCEED 0.12.0", 295 "Use -pressure_infinity to set pressure at boundary layer edge and -idl_pressure to set the IDL reference pressure")); 296 PetscCall(PetscOptionsScalar("-pressure_infinity", "Pressure at boundary layer edge", NULL, P_inf, &P_inf, NULL)); 297 PetscCall(PetscOptionsScalar("-temperature_wall", "Temperature at wall", NULL, T_wall, &T_wall, NULL)); 298 PetscCall(PetscOptionsScalar("-delta0", "Boundary layer height at inflow", NULL, delta0, &delta0, NULL)); 299 PetscCall(PetscOptionsInt("-n_chebyshev", "Number of Chebyshev terms", NULL, N, &N, NULL)); 300 PetscCheck(3 <= N && N <= BLASIUS_MAX_N_CHEBYSHEV, comm, PETSC_ERR_ARG_OUTOFRANGE, "-n_chebyshev %" PetscInt_FMT " must be in range [3, %d]", N, 301 BLASIUS_MAX_N_CHEBYSHEV); 302 if (honee->app_ctx->mesh_transform_type == MESH_TRANSFORM_PLATEMESH) { 303 PetscCall(PetscOptionsBoundedInt("-platemesh_Ndelta", "Velocity at boundary layer edge", NULL, mesh_Ndelta, &mesh_Ndelta, NULL, 1)); 304 PetscCall(PetscOptionsScalar("-platemesh_refine_height", "Height of boundary layer mesh refinement", NULL, mesh_refine_height, 305 &mesh_refine_height, NULL)); 306 PetscCall(PetscOptionsScalar("-platemesh_growth", "Geometric growth rate of boundary layer mesh", NULL, mesh_growth, &mesh_growth, NULL)); 307 PetscCall( 308 PetscOptionsScalar("-platemesh_top_angle", "Geometric top_angle rate of boundary layer mesh", NULL, mesh_top_angle, &mesh_top_angle, NULL)); 309 PetscCall(PetscOptionsString("-platemesh_y_node_locs_path", 310 "Path to file with y node locations. " 311 "If empty, will use the algorithmic mesh warping.", 312 NULL, mesh_ynodes_path, mesh_ynodes_path, sizeof(mesh_ynodes_path), NULL)); 313 } 314 PetscCall(PetscOptionsBool("-stg_use", "Use STG inflow boundary condition", NULL, use_stg, &use_stg, NULL)); 315 PetscOptionsEnd(); 316 317 Units units = honee->units; 318 319 T_inf *= units->Kelvin; 320 T_wall *= units->Kelvin; 321 P_inf *= units->Pascal; 322 U_inf *= units->meter / units->second; 323 delta0 *= units->meter; 324 325 if (honee->app_ctx->mesh_transform_type == MESH_TRANSFORM_PLATEMESH) { 326 PetscReal *mesh_ynodes = NULL; 327 PetscInt mesh_nynodes = 0; 328 if (strcmp(mesh_ynodes_path, "")) PetscCall(GetYNodeLocs(comm, mesh_ynodes_path, &mesh_ynodes, &mesh_nynodes)); 329 PetscCall(ModifyMesh(comm, dm, mesh_growth, mesh_Ndelta, mesh_refine_height, mesh_top_angle, &mesh_ynodes, &mesh_nynodes)); 330 PetscCall(PetscFree(mesh_ynodes)); 331 } 332 333 // Some properties depend on parameters from NewtonianIdealGas 334 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfctx, CEED_MEM_HOST, &newtonian_ig_ctx)); 335 336 StatePrimitive Y_inf = { 337 .pressure = P_inf, .velocity = {U_inf, 0, 0}, 338 .temperature = T_inf 339 }; 340 State S_infty = StateFromPrimitive(newtonian_ig_ctx->gas, Y_inf); 341 342 blasius_ctx->weakT = weakT; 343 blasius_ctx->T_wall = T_wall; 344 blasius_ctx->delta0 = delta0; 345 blasius_ctx->S_infty = S_infty; 346 blasius_ctx->n_cheb = N; 347 blasius_ctx->implicit = honee->phys->implicit; 348 if (P0_set) newtonian_ig_ctx->idl_pressure = P_inf; // For maintaining behavior of -P0 flag (which is deprecated) 349 blasius_ctx->newt_ctx = *newtonian_ig_ctx; 350 351 { 352 PetscReal domain_min[3], domain_max[3]; 353 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max)); 354 blasius_ctx->x_inflow = domain_min[0]; 355 blasius_ctx->eta_max = 5 * domain_max[1] / blasius_ctx->delta0; 356 } 357 PetscBool diff_filter_mms = PETSC_FALSE; 358 PetscCall(PetscOptionsGetBool(NULL, NULL, "-diff_filter_mms", &diff_filter_mms, NULL)); 359 if (!use_stg && !diff_filter_mms) PetscCall(ComputeChebyshevCoefficients(blasius_ctx)); 360 361 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfctx, &newtonian_ig_ctx)); 362 363 PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &blasius_qfctx)); 364 PetscCallCeed(ceed, CeedQFunctionContextSetData(blasius_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*blasius_ctx), blasius_ctx)); 365 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(blasius_qfctx, CEED_MEM_HOST, FreeContextPetsc)); 366 367 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&problem->ics.qfctx)); 368 problem->ics.qfctx = blasius_qfctx; 369 if (use_stg) { 370 PetscCall(SetupStg(comm, dm, problem, honee, weakT, S_infty.Y.temperature, S_infty.Y.pressure)); 371 } else if (diff_filter_mms) { 372 PetscCall(DifferentialFilterMmsICSetup(honee)); 373 } else { 374 PetscCheck((honee->phys->state_var == STATEVAR_CONSERVATIVE) || (honee->app_ctx->test_type == TESTTYPE_DIFF_FILTER), honee->comm, 375 PETSC_ERR_ARG_INCOMP, "Can only use conservative variables with Blasius and weak inflow"); 376 for (PetscCount b = 0; b < problem->num_bc_defs; b++) { 377 BCDefinition bc_def = problem->bc_defs[b]; 378 const char *name; 379 380 PetscCall(BCDefinitionGetInfo(bc_def, &name, NULL, NULL)); 381 if (!strcmp(name, "inflow")) { 382 HoneeBCStruct honee_bc; 383 384 PetscCall(PetscNew(&honee_bc)); 385 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(blasius_qfctx, &honee_bc->qfctx)); 386 honee_bc->honee = honee; 387 honee_bc->num_comps_jac_data = 0; 388 PetscCall(BCDefinitionSetContext(bc_def, HoneeBCDestroy, honee_bc)); 389 390 PetscCall(BCDefinitionSetIFunction(bc_def, BlasiusInflowBCSetup_CreateIFunctionQF, HoneeBCAddIFunctionOp)); 391 PetscCall(BCDefinitionSetIJacobian(bc_def, BlasiusInflowBCSetup_CreateIJacobianQF, HoneeBCAddIJacobianOp)); 392 } 393 } 394 } 395 PetscFunctionReturn(PETSC_SUCCESS); 396 } 397