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