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