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 typedef struct BlasiusContext_ *BlasiusContext; 21 struct BlasiusContext_ { 22 bool implicit; // !< Using implicit timesteping or not 23 bool weakT; // !< flag to set Temperature weakly at inflow 24 CeedScalar delta0; // !< Boundary layer height at inflow 25 CeedScalar U_inf; // !< Velocity at boundary layer edge 26 CeedScalar T_inf; // !< Temperature at boundary layer edge 27 CeedScalar T_wall; // !< Temperature at the wall 28 CeedScalar P0; // !< Pressure at outflow 29 CeedScalar x_inflow; // !< Location of inflow in x 30 CeedScalar n_cheb; // !< Number of Chebyshev terms 31 CeedScalar *X; // !< Chebyshev polynomial coordinate vector 32 CeedScalar eta_max; // !< Maximum eta in the domain 33 CeedScalar *Tf_cheb; // !< Chebyshev coefficient for f 34 CeedScalar *Th_cheb; // !< Chebyshev coefficient for h 35 struct NewtonianIdealGasContext_ newtonian_ctx; 36 }; 37 38 // ***************************************************************************** 39 // This helper function evaluates Chebyshev polynomials with a set of 40 // coefficients with all their derivatives represented as a recurrence table. 41 // ***************************************************************************** 42 CEED_QFUNCTION_HELPER void ChebyshevEval(int N, const double *Tf, double x, 43 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}, {0, 1, 4*x}, {0, 0, 4}, {0, 0, 0} 48 }; 49 for (int i=0; i<4; i++) { 50 // i-th derivative of f 51 f[i] = table[i][0] * Tf[0] + table[i][1] * Tf[1] + table[i][2] * Tf[2]; 52 } 53 for (int i=3; i<N; i++) { 54 // T_n(x) = 2xT_{n-1}(x) - T_{n-2}(x) 55 table[0][i%3] = 2 * x * table[0][(i-1) % 3] - table[0][(i-2)%3]; 56 // Differentiate Chebyshev polynomials with the recurrence relation 57 for (int j=1; j<4; j++) { 58 // T'_{n}(x)/n = 2T_{n-1}(x) + T'_{n-2}(x)/n-2 59 table[j][i%3] = i * (2 * table[j-1][(i-1) % 3] + table[j][(i-2)%3] / (i-2)); 60 } 61 for (int j=0; j<4; j++) { 62 f[j] += table[j][i%3] * Tf[i]; 63 } 64 } 65 for (int i=1; i<4; i++) { 66 // Transform derivatives from Chebyshev [-1, 1] to [0, eta_max]. 67 for (int j=0; j<i; j++) f[i] *= dX_deta; 68 } 69 } 70 71 // ***************************************************************************** 72 // This helper function computes the Blasius boundary layer solution. 73 // ***************************************************************************** 74 State CEED_QFUNCTION_HELPER(BlasiusSolution)(const BlasiusContext blasius, 75 const CeedScalar x[3], const CeedScalar x0, const CeedScalar x_inflow, 76 const CeedScalar rho, CeedScalar *t12) { 77 CeedInt N = blasius->n_cheb; 78 CeedScalar nu = blasius->newtonian_ctx.mu / rho; 79 CeedScalar eta = x[1]*sqrt(blasius->U_inf/(nu*(x0+x[0]-x_inflow))); 80 CeedScalar X = 2 * (eta / blasius->eta_max) - 1.; 81 CeedScalar U_inf = blasius->U_inf; 82 CeedScalar Rd = GasConstant(&blasius->newtonian_ctx); 83 84 CeedScalar f[4], h[4]; 85 ChebyshevEval(N, blasius->Tf_cheb, X, blasius->eta_max, f); 86 ChebyshevEval(N-1, blasius->Th_cheb, X, blasius->eta_max, h); 87 88 *t12 = rho*nu*U_inf*f[2]*sqrt(U_inf/(nu*(x0+x[0]-x_inflow))); 89 90 CeedScalar Y[5]; 91 Y[1] = U_inf * f[1]; 92 Y[2] = 0.5*sqrt(nu*U_inf/(x0+x[0]-x_inflow))*(eta*f[1] - f[0]); 93 Y[3] = 0.; 94 Y[4] = blasius->T_inf * h[0]; 95 Y[0] = rho * Rd * Y[4]; 96 return StateFromY(&blasius->newtonian_ctx, Y, x); 97 } 98 99 // ***************************************************************************** 100 // This QFunction sets a Blasius boundary layer for the initial condition 101 // ***************************************************************************** 102 CEED_QFUNCTION(ICsBlasius)(void *ctx, CeedInt Q, 103 const CeedScalar *const *in, CeedScalar *const *out) { 104 // Inputs 105 const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 106 107 // Outputs 108 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 109 110 const BlasiusContext context = (BlasiusContext)ctx; 111 const CeedScalar cv = context->newtonian_ctx.cv; 112 const CeedScalar mu = context->newtonian_ctx.mu; 113 const CeedScalar T_inf = context->T_inf; 114 const CeedScalar P0 = context->P0; 115 const CeedScalar delta0 = context->delta0; 116 const CeedScalar U_inf = context->U_inf; 117 const CeedScalar x_inflow = context->x_inflow; 118 const CeedScalar gamma = HeatCapacityRatio(&context->newtonian_ctx); 119 const CeedScalar e_internal = cv * T_inf; 120 const CeedScalar rho = P0 / ((gamma - 1) * e_internal); 121 const CeedScalar x0 = U_inf*rho / (mu*25/(delta0*delta0)); 122 CeedScalar t12; 123 124 // Quadrature Point Loop 125 CeedPragmaSIMD 126 for (CeedInt i=0; i<Q; i++) { 127 const CeedScalar x[3] = {X[0][i], X[1][i], 0.}; 128 State s = BlasiusSolution(context, x, x0, x_inflow, rho, &t12); 129 CeedScalar q[5] = {0}; 130 UnpackState_U(s.U, q); 131 for (CeedInt j=0; j<5; j++) q0[j][i] = q[j]; 132 133 } // End of Quadrature Point Loop 134 return 0; 135 } 136 137 // ***************************************************************************** 138 CEED_QFUNCTION(Blasius_Inflow)(void *ctx, CeedInt Q, 139 const CeedScalar *const *in, 140 CeedScalar *const *out) { 141 // *INDENT-OFF* 142 // Inputs 143 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 144 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 145 (*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 // *INDENT-ON* 150 const BlasiusContext context = (BlasiusContext)ctx; 151 const bool implicit = context->implicit; 152 const CeedScalar mu = context->newtonian_ctx.mu; 153 const CeedScalar cv = context->newtonian_ctx.cv; 154 const CeedScalar Rd = GasConstant(&context->newtonian_ctx); 155 const CeedScalar gamma = HeatCapacityRatio(&context->newtonian_ctx); 156 const CeedScalar T_inf = context->T_inf; 157 const CeedScalar P0 = context->P0; 158 const CeedScalar delta0 = context->delta0; 159 const CeedScalar U_inf = context->U_inf; 160 const CeedScalar x_inflow = context->x_inflow; 161 const bool weakT = context->weakT; 162 const CeedScalar rho_0 = P0 / (Rd * T_inf); 163 const CeedScalar x0 = U_inf*rho_0 / (mu*25/ Square(delta0)); 164 165 CeedPragmaSIMD 166 // Quadrature Point Loop 167 for (CeedInt i=0; i<Q; i++) { 168 // Setup 169 // -- Interp-to-Interp q_data 170 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 171 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 172 // We can effect this by swapping the sign on this weight 173 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 174 175 // Calculate inflow values 176 const CeedScalar x[3] = {X[0][i], X[1][i], 0.}; 177 CeedScalar t12; 178 State s = BlasiusSolution(context, x, x0, x_inflow, rho_0, &t12); 179 180 // enabling user to choose between weak T and weak rho inflow 181 CeedScalar rho,E_internal, P, E_kinetic; 182 if (weakT) { 183 // rho should be from the current solution 184 rho = q[0][i]; 185 // Temperature is being set weakly (T_inf) and for constant cv this sets E_internal 186 E_internal = rho * cv * T_inf; 187 // Find pressure using 188 P = rho*Rd*T_inf; // interior rho with exterior T 189 E_kinetic = .5 * rho * Dot3(s.Y.velocity, s.Y.velocity); 190 } else { 191 // Fixing rho weakly on the inflow to a value consistent with T_inf and P0 192 rho = rho_0; 193 E_kinetic = .5 * rho * Dot3(s.Y.velocity, s.Y.velocity); 194 E_internal = q[4][i] - E_kinetic; // uses set rho and u but E from solution 195 P = E_internal * (gamma - 1.); 196 } 197 const CeedScalar E = E_internal + E_kinetic; 198 // ---- Normal vect 199 const CeedScalar norm[3] = {q_data_sur[1][i], 200 q_data_sur[2][i], 201 q_data_sur[3][i] 202 }; 203 204 // The Physics 205 // Zero v so all future terms can safely sum into it 206 for (CeedInt j=0; j<5; j++) v[j][i] = 0.; 207 208 const CeedScalar u_normal = Dot3(norm, s.Y.velocity); 209 const CeedScalar viscous_flux[3] = {-t12 *norm[1], -t12 *norm[0], 0}; 210 211 // The Physics 212 // -- Density 213 v[0][i] -= wdetJb * rho * u_normal; // interior rho 214 215 // -- Momentum 216 for (CeedInt j=0; j<3; j++) 217 v[j+1][i] -= wdetJb * (rho * u_normal * s.Y.velocity[j] // interior rho 218 + norm[j] * P // mixed P 219 + viscous_flux[j]); 220 221 // -- Total Energy Density 222 v[4][i] -= wdetJb * (u_normal * (E + P) + Dot3(viscous_flux, s.Y.velocity)); 223 224 } // End Quadrature Point Loop 225 return 0; 226 } 227 228 // ***************************************************************************** 229 CEED_QFUNCTION(Blasius_Inflow_Jacobian)(void *ctx, CeedInt Q, 230 const CeedScalar *const *in, 231 CeedScalar *const *out) { 232 // *INDENT-OFF* 233 // Inputs 234 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 235 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 236 (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 237 238 // Outputs 239 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 240 // *INDENT-ON* 241 const BlasiusContext context = (BlasiusContext)ctx; 242 const bool implicit = context->implicit; 243 const CeedScalar mu = context->newtonian_ctx.mu; 244 const CeedScalar cv = context->newtonian_ctx.cv; 245 const CeedScalar Rd = GasConstant(&context->newtonian_ctx); 246 const CeedScalar gamma = HeatCapacityRatio(&context->newtonian_ctx); 247 const CeedScalar T_inf = context->T_inf; 248 const CeedScalar P0 = context->P0; 249 const CeedScalar delta0 = context->delta0; 250 const CeedScalar U_inf = context->U_inf; 251 const bool weakT = context->weakT; 252 const CeedScalar rho_0 = P0 / (Rd * T_inf); 253 const CeedScalar x0 = U_inf*rho_0 / (mu*25/ (delta0*delta0)); 254 255 CeedPragmaSIMD 256 // Quadrature Point Loop 257 for (CeedInt i=0; i<Q; i++) { 258 // Setup 259 // -- Interp-to-Interp q_data 260 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 261 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 262 // We can effect this by swapping the sign on this weight 263 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 264 265 // Calculate inflow values 266 const CeedScalar x[3] = {X[0][i], X[1][i], 0.}; 267 CeedScalar t12; 268 State s = BlasiusSolution(context, x, x0, 0, rho_0, &t12); 269 270 // enabling user to choose between weak T and weak rho inflow 271 CeedScalar drho, dE, dP; 272 if (weakT) { 273 // rho should be from the current solution 274 drho = dq[0][i]; 275 CeedScalar dE_internal = drho * cv * T_inf; 276 CeedScalar dE_kinetic = .5 * drho * Dot3(s.Y.velocity, s.Y.velocity); 277 dE = dE_internal + dE_kinetic; 278 dP = drho * Rd * T_inf; // interior rho with exterior T 279 } else { // rho specified, E_internal from solution 280 drho = 0; 281 dE = dq[4][i]; 282 dP = dE * (gamma - 1.); 283 } 284 const CeedScalar norm[3] = {q_data_sur[1][i], 285 q_data_sur[2][i], 286 q_data_sur[3][i] 287 }; 288 289 const CeedScalar u_normal = Dot3(norm, s.Y.velocity); 290 291 v[0][i] = - wdetJb * drho * u_normal; 292 for (int j=0; j<3; j++) 293 v[j+1][i] = -wdetJb * (drho * u_normal * s.Y.velocity[j] + norm[j] * dP); 294 v[4][i] = - wdetJb * u_normal * (dE + dP); 295 } // End Quadrature Point Loop 296 return 0; 297 } 298 299 #endif // blasius_h 300