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 12 #ifndef blasius_h 13 #define blasius_h 14 15 #include <ceed.h> 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 42 // coefficients with all their derivatives represented as a recurrence table. 43 // ***************************************************************************** 44 CEED_QFUNCTION_HELPER void ChebyshevEval(int N, const double *Tf, double x, 45 double eta_max, double *f) { 46 double dX_deta = 2 / eta_max; 47 double table[4][3] = { 48 // Chebyshev polynomials T_0, T_1, T_2 of the first kind in (-1,1) 49 {1, x, 2*x *x - 1}, {0, 1, 4*x}, {0, 0, 4}, {0, 0, 0} 50 }; 51 for (int i=0; i<4; i++) { 52 // i-th derivative of f 53 f[i] = table[i][0] * Tf[0] + table[i][1] * Tf[1] + table[i][2] * Tf[2]; 54 } 55 for (int i=3; i<N; i++) { 56 // T_n(x) = 2xT_{n-1}(x) - T_{n-2}(x) 57 table[0][i%3] = 2 * x * table[0][(i-1) % 3] - table[0][(i-2)%3]; 58 // Differentiate Chebyshev polynomials with the recurrence relation 59 for (int j=1; j<4; j++) { 60 // T'_{n}(x)/n = 2T_{n-1}(x) + T'_{n-2}(x)/n-2 61 table[j][i%3] = i * (2 * table[j-1][(i-1) % 3] + table[j][(i-2)%3] / (i-2)); 62 } 63 for (int j=0; j<4; j++) { 64 f[j] += table[j][i%3] * Tf[i]; 65 } 66 } 67 for (int i=1; i<4; i++) { 68 // Transform derivatives from Chebyshev [-1, 1] to [0, eta_max]. 69 for (int j=0; j<i; j++) f[i] *= dX_deta; 70 } 71 } 72 73 // ***************************************************************************** 74 // This helper function computes the Blasius boundary layer solution. 75 // ***************************************************************************** 76 State CEED_QFUNCTION_HELPER(BlasiusSolution)(const BlasiusContext blasius, 77 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, 106 const CeedScalar *const *in, CeedScalar *const *out) { 107 // Inputs 108 const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 109 110 // Outputs 111 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 112 113 const BlasiusContext context = (BlasiusContext)ctx; 114 const CeedScalar cv = context->newtonian_ctx.cv; 115 const CeedScalar mu = context->newtonian_ctx.mu; 116 const CeedScalar T_inf = context->T_inf; 117 const CeedScalar P0 = context->P0; 118 const CeedScalar delta0 = context->delta0; 119 const CeedScalar U_inf = context->U_inf; 120 const CeedScalar x_inflow = context->x_inflow; 121 const CeedScalar gamma = HeatCapacityRatio(&context->newtonian_ctx); 122 const CeedScalar e_internal = cv * T_inf; 123 const CeedScalar rho = P0 / ((gamma - 1) * e_internal); 124 const CeedScalar x0 = U_inf*rho / (mu*25/(delta0*delta0)); 125 CeedScalar t12; 126 127 // Quadrature Point Loop 128 CeedPragmaSIMD 129 for (CeedInt i=0; i<Q; i++) { 130 const CeedScalar x[3] = {X[0][i], X[1][i], 0.}; 131 State s = BlasiusSolution(context, x, x0, x_inflow, rho, &t12); 132 CeedScalar q[5] = {0}; 133 UnpackState_U(s.U, q); 134 for (CeedInt j=0; j<5; j++) q0[j][i] = q[j]; 135 136 } // End of Quadrature Point Loop 137 return 0; 138 } 139 140 // ***************************************************************************** 141 CEED_QFUNCTION(Blasius_Inflow)(void *ctx, CeedInt Q, 142 const CeedScalar *const *in, 143 CeedScalar *const *out) { 144 // *INDENT-OFF* 145 // Inputs 146 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 147 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 148 (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 149 150 // Outputs 151 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 152 // *INDENT-ON* 153 const BlasiusContext context = (BlasiusContext)ctx; 154 const bool implicit = context->implicit; 155 NewtonianIdealGasContext gas = &context->newtonian_ctx; 156 const CeedScalar mu = context->newtonian_ctx.mu; 157 const CeedScalar Rd = GasConstant(&context->newtonian_ctx); 158 const CeedScalar T_inf = context->T_inf; 159 const CeedScalar P0 = context->P0; 160 const CeedScalar delta0 = context->delta0; 161 const CeedScalar U_inf = context->U_inf; 162 const CeedScalar x_inflow = context->x_inflow; 163 const bool weakT = context->weakT; 164 const CeedScalar rho_0 = P0 / (Rd * T_inf); 165 const CeedScalar x0 = U_inf*rho_0 / (mu*25/ Square(delta0)); 166 167 CeedPragmaSIMD 168 // Quadrature Point Loop 169 for (CeedInt i=0; i<Q; i++) { 170 // Setup 171 // -- Interp-to-Interp q_data 172 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 173 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 174 // We can effect this by swapping the sign on this weight 175 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 176 177 // Calculate inflow values 178 const CeedScalar x[3] = {X[0][i], X[1][i], 0.}; 179 CeedScalar t12; 180 State s = BlasiusSolution(context, x, x0, x_inflow, rho_0, &t12); 181 CeedScalar qi[5]; 182 for (CeedInt j=0; j<5; j++) qi[j] = q[j][i]; 183 State s_int = StateFromU(gas, qi, x); 184 185 // enabling user to choose between weak T and weak rho inflow 186 if (weakT) { // density from the current solution 187 s.U.density = s_int.U.density; 188 s.Y = StatePrimitiveFromConservative(gas, s.U, x); 189 } else { // Total energy from current solution 190 s.U.E_total = s_int.U.E_total; 191 s.Y = StatePrimitiveFromConservative(gas, s.U, x); 192 } 193 194 // ---- Normal vect 195 const CeedScalar norm[3] = {q_data_sur[1][i], 196 q_data_sur[2][i], 197 q_data_sur[3][i] 198 }; 199 200 StateConservative Flux_inviscid[3]; 201 FluxInviscid(&context->newtonian_ctx, s, Flux_inviscid); 202 203 const CeedScalar stress[3][3] = {{0, t12, 0}, {t12, 0, 0}, {0, 0, 0}}; 204 const CeedScalar Fe[3] = {0}; // TODO: viscous energy flux needs grad temperature 205 CeedScalar Flux[5]; 206 FluxTotal_Boundary(Flux_inviscid, stress, Fe, norm, Flux); 207 for (CeedInt j=0; j<5; j++) 208 v[j][i] = -wdetJb * Flux[j]; 209 } // End Quadrature Point Loop 210 return 0; 211 } 212 213 // ***************************************************************************** 214 CEED_QFUNCTION(Blasius_Inflow_Jacobian)(void *ctx, CeedInt Q, 215 const CeedScalar *const *in, 216 CeedScalar *const *out) { 217 // *INDENT-OFF* 218 // Inputs 219 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 220 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 221 (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 222 223 // Outputs 224 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 225 // *INDENT-ON* 226 const BlasiusContext context = (BlasiusContext)ctx; 227 const bool implicit = context->implicit; 228 const CeedScalar mu = context->newtonian_ctx.mu; 229 const CeedScalar cv = context->newtonian_ctx.cv; 230 const CeedScalar Rd = GasConstant(&context->newtonian_ctx); 231 const CeedScalar gamma = HeatCapacityRatio(&context->newtonian_ctx); 232 const CeedScalar T_inf = context->T_inf; 233 const CeedScalar P0 = context->P0; 234 const CeedScalar delta0 = context->delta0; 235 const CeedScalar U_inf = context->U_inf; 236 const bool weakT = context->weakT; 237 const CeedScalar rho_0 = P0 / (Rd * T_inf); 238 const CeedScalar x0 = U_inf*rho_0 / (mu*25/ (delta0*delta0)); 239 240 CeedPragmaSIMD 241 // Quadrature Point Loop 242 for (CeedInt i=0; i<Q; i++) { 243 // Setup 244 // -- Interp-to-Interp q_data 245 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 246 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 247 // We can effect this by swapping the sign on this weight 248 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 249 250 // Calculate inflow values 251 const CeedScalar x[3] = {X[0][i], X[1][i], X[2][i]}; 252 CeedScalar t12; 253 State s = BlasiusSolution(context, x, x0, 0, rho_0, &t12); 254 255 // enabling user to choose between weak T and weak rho inflow 256 CeedScalar drho, dE, dP; 257 if (weakT) { 258 // rho should be from the current solution 259 drho = dq[0][i]; 260 CeedScalar dE_internal = drho * cv * T_inf; 261 CeedScalar dE_kinetic = .5 * drho * Dot3(s.Y.velocity, s.Y.velocity); 262 dE = dE_internal + dE_kinetic; 263 dP = drho * Rd * T_inf; // interior rho with exterior T 264 } else { // rho specified, E_internal from solution 265 drho = 0; 266 dE = dq[4][i]; 267 dP = dE * (gamma - 1.); 268 } 269 const CeedScalar norm[3] = {q_data_sur[1][i], 270 q_data_sur[2][i], 271 q_data_sur[3][i] 272 }; 273 274 const CeedScalar u_normal = Dot3(norm, s.Y.velocity); 275 276 v[0][i] = - wdetJb * drho * u_normal; 277 for (int j=0; j<3; j++) 278 v[j+1][i] = -wdetJb * (drho * u_normal * s.Y.velocity[j] + norm[j] * dP); 279 v[4][i] = - wdetJb * u_normal * (dE + dP); 280 } // End Quadrature Point Loop 281 return 0; 282 } 283 284 #endif // blasius_h 285