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