1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3
4 /// @file
5 /// Operator for HONEE
6 #include <ceed/types.h>
7
8 #include "newtonian_state.h"
9 #include "newtonian_types.h"
10 #include "utils.h"
11
12 #define BLASIUS_MAX_N_CHEBYSHEV 50
13
14 typedef struct BlasiusContext_ *BlasiusContext;
15 struct BlasiusContext_ {
16 bool implicit; // !< Using implicit timesteping or not
17 bool weakT; // !< flag to set Temperature weakly at inflow
18 CeedScalar delta0; // !< Boundary layer height at inflow
19 State S_infty;
20 CeedScalar T_wall; // !< Temperature at the wall
21 CeedScalar x_inflow; // !< Location of inflow in x
22 CeedScalar n_cheb; // !< Number of Chebyshev terms
23 CeedScalar *X; // !< Chebyshev polynomial coordinate vector (CPU only)
24 CeedScalar eta_max; // !< Maximum eta in the domain
25 CeedScalar Tf_cheb[BLASIUS_MAX_N_CHEBYSHEV]; // !< Chebyshev coefficient for f
26 CeedScalar Th_cheb[BLASIUS_MAX_N_CHEBYSHEV - 1]; // !< Chebyshev coefficient for h
27 struct NewtonianIdealGasContext_ newt_ctx;
28 };
29
30 // *****************************************************************************
31 // This helper function evaluates Chebyshev polynomials with a set of coefficients with all their derivatives represented as a recurrence table.
32 // *****************************************************************************
ChebyshevEval(int N,const double * Tf,double x,double eta_max,double * f)33 CEED_QFUNCTION_HELPER void ChebyshevEval(int N, const double *Tf, double x, double eta_max, double *f) {
34 double dX_deta = 2 / eta_max;
35 // Chebyshev polynomials T_0, T_1, T_2 of the first kind in (-1,1)
36 double table[4][3] = {
37 {1, x, 2 * x * x - 1},
38 {0, 1, 4 * x },
39 {0, 0, 4 },
40 {0, 0, 0 }
41 };
42 for (int i = 0; i < 4; i++) {
43 // i-th derivative of f
44 f[i] = table[i][0] * Tf[0] + table[i][1] * Tf[1] + table[i][2] * Tf[2];
45 }
46 for (int i = 3; i < N; i++) {
47 // T_n(x) = 2xT_{n-1}(x) - T_{n-2}(x)
48 table[0][i % 3] = 2 * x * table[0][(i - 1) % 3] - table[0][(i - 2) % 3];
49 // Differentiate Chebyshev polynomials with the recurrence relation
50 for (int j = 1; j < 4; j++) {
51 // T'_{n}(x)/n = 2T_{n-1}(x) + T'_{n-2}(x)/n-2
52 table[j][i % 3] = i * (2 * table[j - 1][(i - 1) % 3] + table[j][(i - 2) % 3] / (i - 2));
53 }
54 for (int j = 0; j < 4; j++) {
55 f[j] += table[j][i % 3] * Tf[i];
56 }
57 }
58 for (int i = 1; i < 4; i++) {
59 // Transform derivatives from Chebyshev [-1, 1] to [0, eta_max].
60 for (int j = 0; j < i; j++) f[i] *= dX_deta;
61 }
62 }
63
64 // *****************************************************************************
65 // This helper function computes the Blasius boundary layer solution.
66 // *****************************************************************************
BlasiusSolution(const BlasiusContext blasius,const CeedScalar x[3],const CeedScalar x0,const CeedScalar x_inflow,const CeedScalar rho_infty,CeedScalar * t12)67 State CEED_QFUNCTION_HELPER(BlasiusSolution)(const BlasiusContext blasius, const CeedScalar x[3], const CeedScalar x0, const CeedScalar x_inflow,
68 const CeedScalar rho_infty, CeedScalar *t12) {
69 CeedInt N = blasius->n_cheb;
70 CeedScalar mu = blasius->newt_ctx.gas.mu;
71 State S_infty = blasius->S_infty;
72 CeedScalar nu = mu / rho_infty;
73 CeedScalar U_infty = Norm3(S_infty.Y.velocity);
74 CeedScalar eta = x[1] * sqrt(U_infty / (nu * (x0 + x[0] - x_inflow)));
75 CeedScalar X = 2 * (eta / blasius->eta_max) - 1.;
76 CeedScalar Rd = GasConstant(blasius->newt_ctx.gas);
77
78 CeedScalar f[4], h[4];
79 ChebyshevEval(N, blasius->Tf_cheb, X, blasius->eta_max, f);
80 ChebyshevEval(N - 1, blasius->Th_cheb, X, blasius->eta_max, h);
81
82 *t12 = mu * U_infty * f[2] * sqrt(U_infty / (nu * (x0 + x[0] - x_inflow)));
83
84 CeedScalar Y[5];
85 Y[1] = U_infty * f[1];
86 Y[2] = 0.5 * sqrt(nu * U_infty / (x0 + x[0] - x_inflow)) * (eta * f[1] - f[0]);
87 Y[3] = 0.;
88 Y[4] = S_infty.Y.temperature * h[0];
89 Y[0] = rho_infty / h[0] * Rd * Y[4];
90 return StateFromY(blasius->newt_ctx.gas, Y);
91 }
92
93 // *****************************************************************************
94 // This QFunction sets a Blasius boundary layer for the initial condition
95 // *****************************************************************************
ICsBlasius(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)96 CEED_QFUNCTION(ICsBlasius)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
97 const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
98 CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
99
100 const BlasiusContext context = (BlasiusContext)ctx;
101 const NewtonianIGProperties gas = context->newt_ctx.gas;
102 const CeedScalar mu = context->newt_ctx.gas.mu;
103 const CeedScalar delta0 = context->delta0;
104 const CeedScalar x_inflow = context->x_inflow;
105 CeedScalar t12;
106
107 const State S_infty = context->S_infty;
108 const CeedScalar U_infty = Norm3(S_infty.Y.velocity);
109
110 const CeedScalar x0 = U_infty * S_infty.U.density / (mu * 25 / Square(delta0));
111
112 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
113 const CeedScalar x[3] = {X[0][i], X[1][i], X[2][i]};
114 State s = BlasiusSolution(context, x, x0, x_inflow, S_infty.U.density, &t12);
115 CeedScalar q[5];
116
117 StateToQ(gas, s, q, context->newt_ctx.state_var);
118 for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j];
119 }
120 return 0;
121 }
122
123 // *****************************************************************************
Blasius_Inflow(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)124 CEED_QFUNCTION(Blasius_Inflow)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
125 const BlasiusContext context = (BlasiusContext)ctx;
126 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
127 const CeedScalar(*q_data_sur) = in[2];
128 const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3];
129 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
130
131 const bool is_implicit = context->implicit;
132 const NewtonianIGProperties gas = context->newt_ctx.gas;
133 State S_infty = context->S_infty;
134 const CeedScalar rho_0 = S_infty.U.density;
135 const CeedScalar U_infty = Norm3(S_infty.Y.velocity);
136 const CeedScalar x0 = U_infty * rho_0 / (gas.mu * 25 / Square(context->delta0));
137
138 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
139 CeedScalar wdetJb, normal[3];
140 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, normal);
141 wdetJb *= is_implicit ? -1. : 1.;
142
143 // Calculate inflow values
144 const CeedScalar x[3] = {X[0][i], X[1][i], 0.};
145 CeedScalar t12;
146 State s = BlasiusSolution(context, x, x0, context->x_inflow, rho_0, &t12);
147 CeedScalar qi[5];
148 for (CeedInt j = 0; j < 5; j++) qi[j] = q[j][i];
149 State s_int = StateFromU(gas, qi);
150
151 // enabling user to choose between weak T and weak rho inflow
152 if (context->weakT) { // density from the current solution
153 s.U.density = s_int.U.density;
154 s.Y = StatePrimitiveFromConservative(gas, s.U);
155 } else { // Total energy from current solution
156 s.U.E_total = s_int.U.E_total;
157 s.Y = StatePrimitiveFromConservative(gas, s.U);
158 }
159
160 StateConservative Flux_inviscid[3];
161 FluxInviscid(gas, s, Flux_inviscid);
162
163 const CeedScalar stress[3][3] = {
164 {0, t12, 0},
165 {t12, 0, 0},
166 {0, 0, 0}
167 };
168 const CeedScalar Fe[3] = {0}; // TODO: viscous energy flux needs grad temperature
169 CeedScalar Flux[5];
170 FluxTotal_Boundary(Flux_inviscid, stress, Fe, normal, Flux);
171 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j];
172 }
173 return 0;
174 }
175
176 // *****************************************************************************
Blasius_Inflow_Jacobian(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)177 CEED_QFUNCTION(Blasius_Inflow_Jacobian)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
178 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
179 const CeedScalar(*q_data_sur) = in[2];
180 const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3];
181 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
182
183 const BlasiusContext context = (BlasiusContext)ctx;
184 const NewtonianIGProperties gas = context->newt_ctx.gas;
185 const bool is_implicit = context->implicit;
186 const CeedScalar Rd = GasConstant(gas);
187 const CeedScalar gamma = HeatCapacityRatio(gas);
188 const State S_infty = context->S_infty;
189 const CeedScalar rho_0 = S_infty.U.density;
190 const CeedScalar U_infty = Norm3(S_infty.Y.velocity);
191 const CeedScalar x0 = U_infty * rho_0 / (gas.mu * 25 / Square(context->delta0));
192
193 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
194 CeedScalar wdetJb, normal[3];
195 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, normal);
196 wdetJb *= is_implicit ? -1. : 1.;
197
198 // Calculate inflow values
199 const CeedScalar x[3] = {X[0][i], X[1][i], X[2][i]};
200 CeedScalar t12;
201 State s = BlasiusSolution(context, x, x0, 0, rho_0, &t12);
202
203 // enabling user to choose between weak T and weak rho inflow
204 CeedScalar drho, dE, dP;
205 if (context->weakT) {
206 // rho should be from the current solution
207 drho = dq[0][i];
208 CeedScalar dE_internal = drho * gas.cv * S_infty.Y.temperature;
209 CeedScalar dE_kinetic = .5 * drho * Dot3(s.Y.velocity, s.Y.velocity);
210 dE = dE_internal + dE_kinetic;
211 dP = drho * Rd * S_infty.Y.temperature; // interior rho with exterior T
212 } else {
213 // rho specified, E_internal from solution
214 drho = 0;
215 dE = dq[4][i];
216 dP = dE * (gamma - 1.);
217 }
218
219 const CeedScalar u_normal = Dot3(normal, s.Y.velocity);
220
221 v[0][i] = -wdetJb * drho * u_normal;
222 for (int j = 0; j < 3; j++) {
223 v[j + 1][i] = -wdetJb * (drho * u_normal * s.Y.velocity[j] + normal[j] * dP);
224 }
225 v[4][i] = -wdetJb * u_normal * (dE + dP);
226 }
227 return 0;
228 }
229