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