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 <math.h> 16 #include <ceed.h> 17 #include "newtonian_types.h" 18 19 typedef struct BlasiusContext_ *BlasiusContext; 20 struct BlasiusContext_ { 21 bool implicit; // !< Using implicit timesteping or not 22 bool weakT; // !< flag to set Temperature weakly at inflow 23 CeedScalar delta0; // !< Boundary layer height at inflow 24 CeedScalar Uinf; // !< Velocity at boundary layer edge 25 CeedScalar P0; // !< Pressure at outflow 26 CeedScalar theta0; // !< Temperature at inflow 27 struct NewtonianIdealGasContext_ newtonian_ctx; 28 }; 29 30 #ifndef M_PI 31 #define M_PI 3.14159265358979323846 32 #endif 33 34 void CEED_QFUNCTION_HELPER(BlasiusSolution)(const CeedScalar y, 35 const CeedScalar Uinf, const CeedScalar x0, const CeedScalar x, 36 const CeedScalar rho, CeedScalar *u, CeedScalar *v, CeedScalar *t12, 37 const NewtonianIdealGasContext newt_ctx) { 38 39 CeedInt nprofs = 50; 40 // *INDENT-OFF* 41 CeedScalar eta_table[] = { 42 0.000000000000000000e+00, 1.282051282051281937e-01, 2.564102564102563875e-01, 3.846153846153845812e-01, 5.128205128205127750e-01, 43 6.410256410256409687e-01, 7.692307692307691624e-01, 8.974358974358973562e-01, 1.025641025641025550e+00, 1.153846153846153744e+00, 44 1.282051282051281937e+00, 1.410256410256410131e+00, 1.538461538461538325e+00, 1.666666666666666519e+00, 1.794871794871794712e+00, 45 1.923076923076922906e+00, 2.051282051282051100e+00, 2.179487179487179294e+00, 2.307692307692307487e+00, 2.435897435897435681e+00, 46 2.564102564102563875e+00, 2.692307692307692069e+00, 2.820512820512820262e+00, 2.948717948717948456e+00, 3.076923076923076650e+00, 47 3.205128205128204844e+00, 3.333333333333333037e+00, 3.461538461538461231e+00, 3.589743589743589425e+00, 3.717948717948717618e+00, 48 3.846153846153845812e+00, 3.974358974358974006e+00, 4.102564102564102200e+00, 4.230769230769229949e+00, 4.358974358974358587e+00, 49 4.487179487179487225e+00, 4.615384615384614975e+00, 4.743589743589742724e+00, 4.871794871794871362e+00, 5.000000000000000000e+00, 50 5.500000000000000000e+00, 6.000000000000000000e+00, 6.500000000000000000e+00, 7.000000000000000000e+00, 7.500000000000000000e+00, 51 8.000000000000000000e+00, 8.500000000000000000e+00, 9.000000000000000000e+00, 9.500000000000000000e+00, 1.000000000000000000e+01}; 52 53 CeedScalar f_table[] = { 54 0.000000000000000000e+00, 2.728923405566200267e-03, 1.091524811461423369e-02, 2.455658828897525764e-02, 4.364674649279581820e-02, 55 6.817382707725749835e-02, 9.811838418932711248e-02, 1.334516294237205192e-01, 1.741337304561980659e-01, 2.201122374410622862e-01, 56 2.713206781625860375e-01, 3.276773654929600599e-01, 3.890844612583744255e-01, 4.554273387986328414e-01, 5.265742820946719416e-01, 57 6.023765522220410062e-01, 6.826688421431770237e-01, 7.672701287583111318e-01, 8.559849171804534418e-01, 9.486048570979430661e-01, 58 1.044910695686512625e+00, 1.144674516826549082e+00, 1.247662203367335465e+00, 1.353636048811749593e+00, 1.462357437868362364e+00, 59 1.573589512396551759e+00, 1.687099740622293842e+00, 1.802662313062363353e+00, 1.920060297987626230e+00, 2.039087501786055245e+00, 60 2.159549994377929050e+00, 2.281267275838891884e+00, 2.404073076539093190e+00, 2.527815798402052838e+00, 2.652358618452637540e+00, 61 2.777579287003750341e+00, 2.903369661199559637e+00, 3.029635020019957992e+00, 3.156293209307130088e+00, 3.283273665161465349e+00, 62 3.780571892998292771e+00, 4.279620922520262383e+00, 4.779322325882148448e+00, 5.279238811036782053e+00, 5.779218028455369804e+00, 63 6.279213431354994768e+00, 6.779212528163703233e+00, 7.279212370655419484e+00, 7.779212346288013613e+00, 8.279212342945751146e+00}; 64 65 CeedScalar fp_table[] = { 66 0.000000000000000000e+00, 4.257083277988830267e-02, 8.513297869782740501e-02, 1.276641169537044151e-01, 1.701271279078802878e-01, 67 2.124702831905590783e-01, 2.546276046951935212e-01, 2.965194442747576264e-01, 3.380533304776729975e-01, 3.791251204629754179e-01, 68 4.196204840172004791e-01, 4.594167322894788796e-01, 4.983849866855867838e-01, 5.363926638765821320e-01, 5.733062319885513514e-01, 69 6.089941719927144392e-01, 6.433300586189647507e-01, 6.761956584341198839e-01, 7.074839307288774970e-01, 7.371018110314454530e-01, 70 7.649726585225528064e-01, 7.910382579383948842e-01, 8.152602836158657773e-01, 8.376211573266827415e-01, 8.581242609418713307e-01, 71 8.767934976651666767e-01, 8.936722290953328374e-01, 9.088216471306606037e-01, 9.223186672607004422e-01, 9.342534510898168332e-01, 72 9.447266795705382414e-01, 9.538467037387058367e-01, 9.617266968332524035e-01, 9.684819213624265011e-01, 9.742272083384174719e-01, 73 9.790747253056680810e-01, 9.831320868743089747e-01, 9.865008381344084754e-01, 9.892753192614093249e-01, 9.915419001656551323e-01, 74 9.968788209317821503e-01, 9.989728724371175206e-01, 9.996990677381791812e-01, 9.999216041491896245e-01, 9.999818594083667023e-01, 75 9.999962745365539307e-01, 9.999993214550036980e-01, 9.999998904550418954e-01, 9.999999843329338001e-01, 9.999999980166356384e-01}; 76 77 CeedScalar fpp_table[] = { 78 3.320573362157903663e-01, 3.320379743512646420e-01, 3.319024760665882368e-01, 3.315350015070190337e-01, 3.308206767975666041e-01, 79 3.296466995822193158e-01, 3.279038639411161471e-01, 3.254884713737624113e-01, 3.223045750196085746e-01, 3.182664816607024272e-01, 80 3.133014118810801829e-01, 3.073521951089355775e-01, 3.003798556086043625e-01, 2.923659305537876785e-01, 2.833143548208253981e-01, 81 2.732527514995234941e-01, 2.622329840371728227e-01, 2.503308560706500874e-01, 2.376448876931176457e-01, 2.242941499773744018e-01, 82 2.104151994284793603e-01, 1.961582158440171031e-01, 1.816825052623964043e-01, 1.671515786102889534e-01, 1.527280512426029968e-01, 83 1.385686249977987894e-01, 1.248194106805364800e-01, 1.116118251613979206e-01, 9.905925581301598670e-02, 8.725462988794610575e-02, 84 7.626896310981794158e-02, 6.615089622448211415e-02, 5.692716644118058639e-02, 4.860390768479891377e-02, 4.116863313890323922e-02, 85 3.459272784597366285e-02, 2.883426862493499582e-02, 2.384099224121952881e-02, 1.955324839409207718e-02, 1.590679868531958210e-02, 86 6.578593141419011685e-03, 2.402039843751689954e-03, 7.741093231657678389e-04, 2.201689553063347941e-04, 5.526217815680267893e-05, 87 1.224092624232004387e-05, 2.392841910090350858e-06, 4.127879363882133676e-07, 6.284244603762621373e-08, 8.442944409712819646e-09}; 88 // *INDENT-ON* 89 90 CeedScalar nu = newt_ctx->mu / rho; 91 CeedScalar eta = y*sqrt(Uinf/(nu*(x0+x))); 92 CeedInt idx=-1; 93 94 for(CeedInt i=0; i<nprofs; i++) { 95 if (eta < eta_table[i]) { 96 idx = i; 97 break; 98 } 99 } 100 CeedScalar f, fp, fpp; 101 102 if (idx > 0) { // eta within the bounds of eta_table 103 CeedScalar coeff = (eta - eta_table[idx-1]) / (eta_table[idx] - eta_table[idx 104 -1]); 105 106 f = f_table[idx-1] + coeff*( f_table[idx] - f_table[idx-1] ); 107 fp = fp_table[idx-1] + coeff*( fp_table[idx] - fp_table[idx-1] ); 108 fpp = fpp_table[idx-1] + coeff*( fpp_table[idx] - fpp_table[idx-1] ); 109 } else { // eta outside bounds of eta_table 110 f = f_table[nprofs-1]; 111 fp = fp_table[nprofs-1]; 112 fpp = fpp_table[nprofs-1]; 113 eta = eta_table[nprofs-1]; 114 } 115 116 *u = Uinf*fp; 117 *t12 = rho*nu*Uinf*fpp*sqrt(Uinf/(nu*(x0+x))); 118 *v = 0.5*sqrt(nu*Uinf/(x0+x))*(eta*fp - f); 119 } 120 121 // ***************************************************************************** 122 // This QFunction sets a Blasius boundary layer for the initial condition 123 // ***************************************************************************** 124 CEED_QFUNCTION(ICsBlasius)(void *ctx, CeedInt Q, 125 const CeedScalar *const *in, CeedScalar *const *out) { 126 // Inputs 127 const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 128 129 // Outputs 130 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 131 132 const BlasiusContext context = (BlasiusContext)ctx; 133 const CeedScalar cv = context->newtonian_ctx.cv; 134 const CeedScalar cp = context->newtonian_ctx.cp; 135 const CeedScalar gamma = cp/cv; 136 const CeedScalar mu = context->newtonian_ctx.mu; 137 138 const CeedScalar theta0 = context->theta0; 139 const CeedScalar P0 = context->P0; 140 const CeedScalar delta0 = context->delta0; 141 const CeedScalar Uinf = context->Uinf; 142 143 const CeedScalar e_internal = cv * theta0; 144 const CeedScalar rho = P0 / ((gamma - 1) * e_internal); 145 const CeedScalar x0 = Uinf*rho / (mu*25/ (delta0*delta0) ); 146 CeedScalar u, v, t12; 147 148 // Quadrature Point Loop 149 CeedPragmaSIMD 150 for (CeedInt i=0; i<Q; i++) { 151 const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]}; 152 153 BlasiusSolution(x[1], Uinf, x0, x[0], rho, &u, &v, &t12, 154 &context->newtonian_ctx); 155 156 q0[0][i] = rho; 157 q0[1][i] = u * rho; 158 q0[2][i] = v * rho; 159 q0[3][i] = 0.; 160 q0[4][i] = rho * e_internal + 0.5*(u*u + v*v)*rho; 161 } // End of Quadrature Point Loop 162 return 0; 163 } 164 165 // ***************************************************************************** 166 CEED_QFUNCTION(Blasius_Inflow)(void *ctx, CeedInt Q, 167 const CeedScalar *const *in, 168 CeedScalar *const *out) { 169 // *INDENT-OFF* 170 // Inputs 171 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 172 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1], 173 (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 174 175 // Outputs 176 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 177 // *INDENT-ON* 178 const BlasiusContext context = (BlasiusContext)ctx; 179 const bool implicit = context->implicit; 180 const CeedScalar mu = context->newtonian_ctx.mu; 181 const CeedScalar cv = context->newtonian_ctx.cv; 182 const CeedScalar cp = context->newtonian_ctx.cp; 183 const CeedScalar Rd = cp - cv; 184 const CeedScalar gamma = cp/cv; 185 186 const CeedScalar theta0 = context->theta0; 187 const CeedScalar P0 = context->P0; 188 const CeedScalar delta0 = context->delta0; 189 const CeedScalar Uinf = context->Uinf; 190 const bool weakT = context->weakT; 191 const CeedScalar rho_0 = P0 / (Rd * theta0); 192 const CeedScalar x0 = Uinf*rho_0 / (mu*25/ (delta0*delta0) ); 193 194 CeedPragmaSIMD 195 // Quadrature Point Loop 196 for (CeedInt i=0; i<Q; i++) { 197 // Setup 198 // -- Interp-to-Interp q_data 199 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 200 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 201 // We can effect this by swapping the sign on this weight 202 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 203 204 // Calculate inflow values 205 const CeedScalar x[3] = {X[0][i], X[1][i], X[2][i]}; 206 CeedScalar velocity[3] = {0.}; 207 CeedScalar t12; 208 BlasiusSolution(x[1], Uinf, x0, x[0], rho_0, &velocity[0], &velocity[1], 209 &t12, &context->newtonian_ctx); 210 211 // enabling user to choose between weak T and weak rho inflow 212 CeedScalar rho,E_internal, P, E_kinetic; 213 if (weakT) { 214 // rho should be from the current solution 215 rho = q[0][i]; 216 // Temperature is being set weakly (theta0) and for constant cv this sets E_internal 217 E_internal = rho * cv * theta0; 218 // Find pressure using 219 P = rho*Rd*theta0; // interior rho with exterior T 220 E_kinetic = .5 * rho * Dot3(velocity, velocity); 221 } else { 222 // Fixing rho weakly on the inflow to a value consistent with theta0 and P0 223 rho = rho_0; 224 E_kinetic = .5 * rho * Dot3(velocity, velocity); 225 E_internal = q[4][i] - E_kinetic; // uses set rho and u but E from solution 226 P = E_internal * (gamma - 1.); 227 } 228 const CeedScalar E = E_internal + E_kinetic; 229 // ---- Normal vect 230 const CeedScalar norm[3] = {q_data_sur[1][i], 231 q_data_sur[2][i], 232 q_data_sur[3][i] 233 }; 234 235 // The Physics 236 // Zero v so all future terms can safely sum into it 237 for (int j=0; j<5; j++) v[j][i] = 0.; 238 239 const CeedScalar u_normal = Dot3(norm, velocity); 240 241 // The Physics 242 // -- Density 243 v[0][i] -= wdetJb * rho * u_normal; // interior rho 244 245 // -- Momentum 246 for (int j=0; j<3; j++) 247 v[j+1][i] -= wdetJb * (rho * u_normal * velocity[j] + // interior rho 248 norm[j] * P); // mixed P 249 v[2][i] -= wdetJb * t12 ; 250 251 // -- Total Energy Density 252 v[4][i] -= wdetJb * u_normal * (E + P); 253 v[4][i] -= wdetJb * t12 * velocity[1]; 254 255 } // End Quadrature Point Loop 256 return 0; 257 } 258 259 CEED_QFUNCTION(Blasius_Inflow_Jacobian)(void *ctx, CeedInt Q, 260 const CeedScalar *const *in, 261 CeedScalar *const *out) { 262 // *INDENT-OFF* 263 // Inputs 264 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 265 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1], 266 (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 267 268 // Outputs 269 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 270 // *INDENT-ON* 271 const BlasiusContext context = (BlasiusContext)ctx; 272 const bool implicit = context->implicit; 273 const CeedScalar mu = context->newtonian_ctx.mu; 274 const CeedScalar cv = context->newtonian_ctx.cv; 275 const CeedScalar cp = context->newtonian_ctx.cp; 276 const CeedScalar Rd = cp - cv; 277 const CeedScalar gamma = cp/cv; 278 279 const CeedScalar theta0 = context->theta0; 280 const CeedScalar P0 = context->P0; 281 const CeedScalar delta0 = context->delta0; 282 const CeedScalar Uinf = context->Uinf; 283 const bool weakT = context->weakT; 284 const CeedScalar rho_0 = P0 / (Rd * theta0); 285 const CeedScalar x0 = Uinf*rho_0 / (mu*25/ (delta0*delta0) ); 286 287 CeedPragmaSIMD 288 // Quadrature Point Loop 289 for (CeedInt i=0; i<Q; i++) { 290 // Setup 291 // -- Interp-to-Interp q_data 292 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 293 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 294 // We can effect this by swapping the sign on this weight 295 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 296 297 // Calculate inflow values 298 const CeedScalar x[3] = {X[0][i], X[1][i], X[2][i]}; 299 CeedScalar velocity[3] = {0}; 300 CeedScalar t12; 301 BlasiusSolution(x[1], Uinf, x0, x[0], rho_0, &velocity[0], &velocity[1], 302 &t12, &context->newtonian_ctx); 303 304 // enabling user to choose between weak T and weak rho inflow 305 CeedScalar drho, dE, dP; 306 if (weakT) { 307 // rho should be from the current solution 308 drho = dq[0][i]; 309 CeedScalar dE_internal = drho * cv * theta0; 310 CeedScalar dE_kinetic = .5 * drho * Dot3(velocity, velocity); 311 dE = dE_internal + dE_kinetic; 312 dP = drho * Rd * theta0; // interior rho with exterior T 313 } else { // rho specified, E_internal from solution 314 drho = 0; 315 dE = dq[4][i]; 316 dP = dE * (gamma - 1.); 317 } 318 const CeedScalar norm[3] = {q_data_sur[1][i], 319 q_data_sur[2][i], 320 q_data_sur[3][i] 321 }; 322 323 const CeedScalar u_normal = Dot3(norm, velocity); 324 325 v[0][i] = - wdetJb * drho * u_normal; 326 for (int j=0; j<3; j++) 327 v[j+1][i] = -wdetJb * (drho * u_normal * velocity[j] + norm[j] * dP); 328 v[4][i] = - wdetJb * u_normal * (dE + dP); 329 } // End Quadrature Point Loop 330 return 0; 331 } 332 333 // ***************************************************************************** 334 CEED_QFUNCTION(Blasius_Outflow)(void *ctx, CeedInt Q, 335 const CeedScalar *const *in, 336 CeedScalar *const *out) { 337 // *INDENT-OFF* 338 // Inputs 339 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 340 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1], 341 (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 342 // Outputs 343 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 344 (*jac_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[1]; 345 // *INDENT-ON* 346 347 const BlasiusContext context = (BlasiusContext)ctx; 348 const bool implicit = context->implicit; 349 const CeedScalar mu = context->newtonian_ctx.mu; 350 const CeedScalar cv = context->newtonian_ctx.cv; 351 const CeedScalar cp = context->newtonian_ctx.cp; 352 const CeedScalar Rd = cp - cv; 353 354 const CeedScalar theta0 = context->theta0; 355 const CeedScalar P0 = context->P0; 356 const CeedScalar rho_0 = P0 / (Rd*theta0); 357 const CeedScalar delta0 = context->delta0; 358 const CeedScalar Uinf = context->Uinf; 359 const CeedScalar x0 = Uinf*rho_0 / (mu*25/ (delta0*delta0) ); 360 361 CeedPragmaSIMD 362 // Quadrature Point Loop 363 for (CeedInt i=0; i<Q; i++) { 364 // Setup 365 // -- Interp in 366 const CeedScalar rho = q[0][i]; 367 const CeedScalar u[3] = {q[1][i] / rho, 368 q[2][i] / rho, 369 q[3][i] / rho 370 }; 371 const CeedScalar E = q[4][i]; 372 373 // -- Interp-to-Interp q_data 374 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 375 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 376 // We can effect this by swapping the sign on this weight 377 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 378 379 // ---- Normal vect 380 const CeedScalar norm[3] = {q_data_sur[1][i], 381 q_data_sur[2][i], 382 q_data_sur[3][i] 383 }; 384 385 // The Physics 386 // Zero v so all future terms can safely sum into it 387 for (int j=0; j<5; j++) v[j][i] = 0.; 388 389 // Implementing outflow condition 390 const CeedScalar P = P0; // pressure 391 const CeedScalar u_normal = Dot3(norm, u); // Normal velocity 392 393 // Calculate prescribed outflow traction values 394 const CeedScalar x[3] = {X[0][i], X[1][i], X[2][i]}; 395 CeedScalar velocity[3] = {0.}; 396 CeedScalar t12; 397 BlasiusSolution(x[1], Uinf, x0, x[0], rho_0, &velocity[0], &velocity[1], 398 &t12, &context->newtonian_ctx); 399 // The Physics 400 // -- Density 401 v[0][i] -= wdetJb * rho * u_normal; 402 403 // -- Momentum 404 for (int j=0; j<3; j++) 405 v[j+1][i] -= wdetJb *(rho * u_normal * u[j] + norm[j] * P); 406 v[2][i] += wdetJb * t12 ; 407 408 // -- Total Energy Density 409 v[4][i] -= wdetJb * u_normal * (E + P); 410 v[4][i] += wdetJb * t12 * velocity[1]; 411 412 // Save values for Jacobian 413 jac_data_sur[0][i] = rho; 414 jac_data_sur[1][i] = u[0]; 415 jac_data_sur[2][i] = u[1]; 416 jac_data_sur[3][i] = u[2]; 417 jac_data_sur[4][i] = E; 418 } // End Quadrature Point Loop 419 return 0; 420 } 421 422 CEED_QFUNCTION(Blasius_Outflow_Jacobian)(void *ctx, CeedInt Q, 423 const CeedScalar *const *in, 424 CeedScalar *const *out) { 425 // *INDENT-OFF* 426 // Inputs 427 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 428 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1], 429 (*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 430 // Outputs 431 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 432 // *INDENT-ON* 433 434 const BlasiusContext context = (BlasiusContext)ctx; 435 const bool implicit = context->implicit; 436 437 CeedPragmaSIMD 438 // Quadrature Point Loop 439 for (CeedInt i=0; i<Q; i++) { 440 const CeedScalar rho = jac_data_sur[0][i]; 441 const CeedScalar u[3] = {jac_data_sur[1][i], jac_data_sur[2][i], jac_data_sur[3][i]}; 442 const CeedScalar E = jac_data_sur[4][i]; 443 444 const CeedScalar drho = dq[0][i]; 445 const CeedScalar dmomentum[3] = {dq[1][i], dq[2][i], dq[3][i]}; 446 const CeedScalar dE = dq[4][i]; 447 448 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 449 const CeedScalar norm[3] = {q_data_sur[1][i], 450 q_data_sur[2][i], 451 q_data_sur[3][i] 452 }; 453 454 CeedScalar du[3]; 455 for (int j=0; j<3; j++) du[j] = (dmomentum[j] - u[j] * drho) / rho; 456 const CeedScalar u_normal = Dot3(norm, u); 457 const CeedScalar du_normal = Dot3(norm, du); 458 const CeedScalar dmomentum_normal = drho * u_normal + rho * du_normal; 459 const CeedScalar P = context->P0; 460 const CeedScalar dP = 0; 461 462 v[0][i] = -wdetJb * dmomentum_normal; 463 for (int j=0; j<3; j++) 464 v[j+1][i] = -wdetJb * (dmomentum_normal * u[j] + rho * u_normal * du[j]); 465 v[4][i] = -wdetJb * (du_normal * (E + P) + u_normal * (dE + dP)); 466 } // End Quadrature Point Loop 467 return 0; 468 } 469 470 #endif // blasius_h 471