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