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
CompressibleBlasiusResidual(SNES snes,Vec X,Vec R,void * ctx)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
ComputeChebyshevCoefficients(BlasiusContext blasius)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.");
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
BlasiusInflowBCSetup_CreateIFunctionQF(BCDefinition bc_def,CeedQFunction * qf)114 static PetscErrorCode BlasiusInflowBCSetup_CreateIFunctionQF(BCDefinition bc_def, CeedQFunction *qf) {
115 HoneeBCStruct honee_bc;
116
117 PetscFunctionBeginUser;
118 PetscCall(BCDefinitionGetContext(bc_def, &honee_bc));
119 PetscCall(HoneeBCCreateIFunctionQF(bc_def, Blasius_Inflow, Blasius_Inflow_loc, honee_bc->qfctx, qf));
120 PetscFunctionReturn(PETSC_SUCCESS);
121 }
122
BlasiusInflowBCSetup_CreateIJacobianQF(BCDefinition bc_def,CeedQFunction * qf)123 static PetscErrorCode BlasiusInflowBCSetup_CreateIJacobianQF(BCDefinition bc_def, CeedQFunction *qf) {
124 HoneeBCStruct honee_bc;
125
126 PetscFunctionBeginUser;
127 PetscCall(BCDefinitionGetContext(bc_def, &honee_bc));
128 PetscCall(HoneeBCCreateIJacobianQF(bc_def, Blasius_Inflow_Jacobian, Blasius_Inflow_Jacobian_loc, honee_bc->qfctx, qf));
129 PetscFunctionReturn(PETSC_SUCCESS);
130 }
131
NS_BLASIUS(ProblemData problem,DM dm,void * ctx)132 PetscErrorCode NS_BLASIUS(ProblemData problem, DM dm, void *ctx) {
133 Honee honee = *(Honee *)ctx;
134 MPI_Comm comm = honee->comm;
135 Ceed ceed = honee->ceed;
136 PetscBool use_stg = PETSC_FALSE;
137 BlasiusContext blasius_ctx;
138 NewtonianIdealGasContext newtonian_ig_ctx;
139 CeedQFunctionContext blasius_qfctx;
140
141 PetscFunctionBeginUser;
142 PetscCall(NS_NEWTONIAN_IG(problem, dm, ctx));
143
144 // ------------------------------------------------------
145 // SET UP Blasius
146 // ------------------------------------------------------
147 problem->ics = (HoneeQFSpec){.qf_func_ptr = ICsBlasius, .qf_loc = ICsBlasius_loc, .qfctx = problem->ics.qfctx};
148
149 CeedScalar U_inf = 40; // m/s
150 CeedScalar T_inf = 288.; // K
151 CeedScalar T_wall = 288.; // K
152 CeedScalar delta0 = 4.2e-3; // m
153 CeedScalar P_inf = 1.01e5; // Pa
154 PetscInt N = 20; // Number of Chebyshev terms
155 PetscBool weakT = PETSC_FALSE; // weak density or temperature
156 PetscBool P0_set;
157
158 PetscOptionsBegin(comm, NULL, "Options for BLASIUS problem", NULL);
159 PetscCall(PetscOptionsBool("-weakT", "Change from rho weak to T weak at inflow", NULL, weakT, &weakT, NULL));
160 PetscCall(PetscOptionsScalar("-velocity_infinity", "Velocity at boundary layer edge", NULL, U_inf, &U_inf, NULL));
161 PetscCall(PetscOptionsScalar("-temperature_infinity", "Temperature at boundary layer edge", NULL, T_inf, &T_inf, NULL));
162 PetscCall(PetscOptionsHasName(NULL, NULL, "-P0", &P0_set)); // For maintaining behavior of -P0 flag (which is deprecated)
163 PetscCall(PetscOptionsDeprecated("-P0", "-pressure_infinity", "libCEED 0.12.0",
164 "Use -pressure_infinity to set pressure at boundary layer edge and -idl_pressure to set the IDL reference "
165 "pressure"));
166 PetscCall(PetscOptionsScalar("-pressure_infinity", "Pressure at boundary layer edge", NULL, P_inf, &P_inf, NULL));
167 PetscCall(PetscOptionsScalar("-temperature_wall", "Temperature at wall", NULL, T_wall, &T_wall, NULL));
168 PetscCall(PetscOptionsScalar("-delta0", "Boundary layer height at inflow", NULL, delta0, &delta0, NULL));
169 PetscCall(PetscOptionsInt("-n_chebyshev", "Number of Chebyshev terms", NULL, N, &N, NULL));
170 PetscCheck(3 <= N && N <= BLASIUS_MAX_N_CHEBYSHEV, comm, PETSC_ERR_ARG_OUTOFRANGE, "-n_chebyshev %" PetscInt_FMT " must be in range [3, %d]", N,
171 BLASIUS_MAX_N_CHEBYSHEV);
172 PetscCall(PetscOptionsBool("-stg_use", "Use STG inflow boundary condition", NULL, use_stg, &use_stg, NULL));
173 PetscOptionsEnd();
174
175 Units units = honee->units;
176
177 T_inf *= units->Kelvin;
178 T_wall *= units->Kelvin;
179 P_inf *= units->Pascal;
180 U_inf *= units->meter / units->second;
181 delta0 *= units->meter;
182
183 // Some properties depend on parameters from NewtonianIdealGas
184 PetscCallCeed(ceed, CeedQFunctionContextGetData(problem->apply_vol_rhs.qfctx, CEED_MEM_HOST, &newtonian_ig_ctx));
185
186 StatePrimitive Y_inf = {
187 .pressure = P_inf, .velocity = {U_inf, 0, 0},
188 .temperature = T_inf
189 };
190 State S_infty = StateFromPrimitive(newtonian_ig_ctx->gas, Y_inf);
191
192 PetscCall(PetscNew(&blasius_ctx));
193 blasius_ctx->weakT = weakT;
194 blasius_ctx->T_wall = T_wall;
195 blasius_ctx->delta0 = delta0;
196 blasius_ctx->S_infty = S_infty;
197 blasius_ctx->n_cheb = N;
198 blasius_ctx->implicit = honee->phys->implicit;
199 if (P0_set) newtonian_ig_ctx->idl_pressure = P_inf; // For maintaining behavior of -P0 flag (which is deprecated)
200 blasius_ctx->newt_ctx = *newtonian_ig_ctx;
201
202 {
203 PetscReal domain_min[3], domain_max[3];
204 PetscCall(DMGetBoundingBox(dm, domain_min, domain_max));
205 blasius_ctx->x_inflow = domain_min[0];
206 blasius_ctx->eta_max = 5 * domain_max[1] / blasius_ctx->delta0;
207 }
208 PetscBool diff_filter_mms = PETSC_FALSE;
209 PetscCall(PetscOptionsGetBool(NULL, NULL, "-diff_filter_mms", &diff_filter_mms, NULL));
210 if (!use_stg && !diff_filter_mms) PetscCall(ComputeChebyshevCoefficients(blasius_ctx));
211
212 PetscCallCeed(ceed, CeedQFunctionContextRestoreData(problem->apply_vol_rhs.qfctx, &newtonian_ig_ctx));
213
214 PetscCallCeed(ceed, CeedQFunctionContextCreate(honee->ceed, &blasius_qfctx));
215 PetscCallCeed(ceed, CeedQFunctionContextSetData(blasius_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*blasius_ctx), blasius_ctx));
216 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(blasius_qfctx, CEED_MEM_HOST, FreeContextPetsc));
217
218 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&problem->ics.qfctx));
219 problem->ics.qfctx = blasius_qfctx;
220 if (use_stg) {
221 PetscCall(SetupStg(comm, dm, problem, honee, weakT, S_infty.Y.temperature, S_infty.Y.pressure));
222 } else if (diff_filter_mms) {
223 PetscCall(DifferentialFilterMmsICSetup(honee));
224 } else {
225 PetscCheck((honee->phys->state_var == STATEVAR_CONSERVATIVE) || (honee->app_ctx->test_type == TESTTYPE_DIFF_FILTER), honee->comm,
226 PETSC_ERR_ARG_INCOMP, "Can only use conservative variables with Blasius and weak inflow");
227 for (PetscCount b = 0; b < problem->num_bc_defs; b++) {
228 BCDefinition bc_def = problem->bc_defs[b];
229 const char *name;
230
231 PetscCall(BCDefinitionGetInfo(bc_def, &name, NULL, NULL));
232 if (!strcmp(name, "inflow")) {
233 HoneeBCStruct honee_bc;
234
235 PetscCall(PetscNew(&honee_bc));
236 PetscCallCeed(ceed, CeedQFunctionContextReferenceCopy(blasius_qfctx, &honee_bc->qfctx));
237 honee_bc->honee = honee;
238 honee_bc->num_comps_jac_data = 0;
239 PetscCall(BCDefinitionSetContext(bc_def, (PetscCtxDestroyFn *)HoneeBCDestroy, honee_bc));
240
241 PetscCall(BCDefinitionSetIFunction(bc_def, BlasiusInflowBCSetup_CreateIFunctionQF, HoneeBCAddIFunctionOp));
242 PetscCall(BCDefinitionSetIJacobian(bc_def, BlasiusInflowBCSetup_CreateIJacobianQF, HoneeBCAddIJacobianOp));
243 }
244 }
245 }
246 PetscFunctionReturn(PETSC_SUCCESS);
247 }
248