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_ newtonian_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 // ***************************************************************************** 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 // ***************************************************************************** 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->newtonian_ctx.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->newtonian_ctx); 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->newtonian_ctx, Y); 91 } 92 93 // ***************************************************************************** 94 // This QFunction sets a Blasius boundary layer for the initial condition 95 // ***************************************************************************** 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 NewtonianIdealGasContext gas = &context->newtonian_ctx; 102 const CeedScalar mu = context->newtonian_ctx.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, gas->state_var); 118 for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j]; 119 } 120 return 0; 121 } 122 123 // ***************************************************************************** 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 CeedScalar(*jac_data_sur) = context->newtonian_ctx.is_implicit ? out[1] : NULL; 131 132 const bool is_implicit = context->implicit; 133 const NewtonianIdealGasContext gas = &context->newtonian_ctx; 134 State S_infty = context->S_infty; 135 const CeedScalar rho_0 = S_infty.U.density; 136 const CeedScalar U_infty = Norm3(S_infty.Y.velocity); 137 const CeedScalar x0 = U_infty * rho_0 / (gas->mu * 25 / Square(context->delta0)); 138 const CeedScalar zeros[11] = {0.}; 139 140 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 141 CeedScalar wdetJb, normal[3]; 142 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, normal); 143 wdetJb *= is_implicit ? -1. : 1.; 144 145 // Calculate inflow values 146 const CeedScalar x[3] = {X[0][i], X[1][i], 0.}; 147 CeedScalar t12; 148 State s = BlasiusSolution(context, x, x0, context->x_inflow, rho_0, &t12); 149 CeedScalar qi[5]; 150 for (CeedInt j = 0; j < 5; j++) qi[j] = q[j][i]; 151 State s_int = StateFromU(gas, qi); 152 153 // enabling user to choose between weak T and weak rho inflow 154 if (context->weakT) { // density from the current solution 155 s.U.density = s_int.U.density; 156 s.Y = StatePrimitiveFromConservative(gas, s.U); 157 } else { // Total energy from current solution 158 s.U.E_total = s_int.U.E_total; 159 s.Y = StatePrimitiveFromConservative(gas, s.U); 160 } 161 162 StateConservative Flux_inviscid[3]; 163 FluxInviscid(&context->newtonian_ctx, s, Flux_inviscid); 164 165 const CeedScalar stress[3][3] = { 166 {0, t12, 0}, 167 {t12, 0, 0}, 168 {0, 0, 0} 169 }; 170 const CeedScalar Fe[3] = {0}; // TODO: viscous energy flux needs grad temperature 171 CeedScalar Flux[5]; 172 FluxTotal_Boundary(Flux_inviscid, stress, Fe, normal, Flux); 173 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j]; 174 if (is_implicit) StoredValuesPack(Q, i, 0, 11, zeros, jac_data_sur); 175 } 176 return 0; 177 } 178 179 // ***************************************************************************** 180 CEED_QFUNCTION(Blasius_Inflow_Jacobian)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 181 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 182 const CeedScalar(*q_data_sur) = in[2]; 183 const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 184 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 185 186 const BlasiusContext context = (BlasiusContext)ctx; 187 const NewtonianIdealGasContext gas = &context->newtonian_ctx; 188 const bool is_implicit = context->implicit; 189 const CeedScalar Rd = GasConstant(gas); 190 const CeedScalar gamma = HeatCapacityRatio(gas); 191 const State S_infty = context->S_infty; 192 const CeedScalar rho_0 = S_infty.U.density; 193 const CeedScalar U_infty = Norm3(S_infty.Y.velocity); 194 const CeedScalar x0 = U_infty * rho_0 / (gas->mu * 25 / Square(context->delta0)); 195 196 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 197 CeedScalar wdetJb, normal[3]; 198 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, normal); 199 wdetJb *= is_implicit ? -1. : 1.; 200 201 // Calculate inflow values 202 const CeedScalar x[3] = {X[0][i], X[1][i], X[2][i]}; 203 CeedScalar t12; 204 State s = BlasiusSolution(context, x, x0, 0, rho_0, &t12); 205 206 // enabling user to choose between weak T and weak rho inflow 207 CeedScalar drho, dE, dP; 208 if (context->weakT) { 209 // rho should be from the current solution 210 drho = dq[0][i]; 211 CeedScalar dE_internal = drho * gas->cv * S_infty.Y.temperature; 212 CeedScalar dE_kinetic = .5 * drho * Dot3(s.Y.velocity, s.Y.velocity); 213 dE = dE_internal + dE_kinetic; 214 dP = drho * Rd * S_infty.Y.temperature; // interior rho with exterior T 215 } else { 216 // rho specified, E_internal from solution 217 drho = 0; 218 dE = dq[4][i]; 219 dP = dE * (gamma - 1.); 220 } 221 222 const CeedScalar u_normal = Dot3(normal, s.Y.velocity); 223 224 v[0][i] = -wdetJb * drho * u_normal; 225 for (int j = 0; j < 3; j++) { 226 v[j + 1][i] = -wdetJb * (drho * u_normal * s.Y.velocity[j] + normal[j] * dP); 227 } 228 v[4][i] = -wdetJb * u_normal * (dE + dP); 229 } 230 return 0; 231 } 232