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 newtonian_h 13 #define newtonian_h 14 15 #include <math.h> 16 #include <ceed.h> 17 #include "newtonian_types.h" 18 #include "newtonian_state.h" 19 #include "utils.h" 20 21 // ***************************************************************************** 22 // Helper function for computing flux Jacobian 23 // ***************************************************************************** 24 CEED_QFUNCTION_HELPER void computeFluxJacobian_NS(CeedScalar dF[3][5][5], 25 const CeedScalar rho, const CeedScalar u[3], const CeedScalar E, 26 const CeedScalar gamma, const CeedScalar g[3], const CeedScalar x[3]) { 27 CeedScalar u_sq = u[0]*u[0] + u[1]*u[1] + u[2]*u[2]; // Velocity square 28 CeedScalar e_potential = -(g[0]*x[0] + g[1]*x[1] + g[2]*x[2]); 29 for (CeedInt i=0; i<3; i++) { // Jacobian matrices for 3 directions 30 for (CeedInt j=0; j<3; j++) { // Rows of each Jacobian matrix 31 dF[i][j+1][0] = ((i==j) ? ((gamma-1.)*(u_sq/2. - e_potential)) : 0.) - 32 u[i]*u[j]; 33 for (CeedInt k=0; k<3; k++) { // Columns of each Jacobian matrix 34 dF[i][0][k+1] = ((i==k) ? 1. : 0.); 35 dF[i][j+1][k+1] = ((j==k) ? u[i] : 0.) + 36 ((i==k) ? u[j] : 0.) - 37 ((i==j) ? u[k] : 0.) * (gamma-1.); 38 dF[i][4][k+1] = ((i==k) ? (E*gamma/rho - (gamma-1.)*u_sq/2.) : 0.) - 39 (gamma-1.)*u[i]*u[k]; 40 } 41 dF[i][j+1][4] = ((i==j) ? (gamma-1.) : 0.); 42 } 43 dF[i][4][0] = u[i] * ((gamma-1.)*u_sq - E*gamma/rho); 44 dF[i][4][4] = u[i] * gamma; 45 } 46 } 47 48 // ***************************************************************************** 49 // Helper function for computing flux Jacobian of Primitive variables 50 // ***************************************************************************** 51 CEED_QFUNCTION_HELPER void computeFluxJacobian_NSp(CeedScalar dF[3][5][5], 52 const CeedScalar rho, const CeedScalar u[3], const CeedScalar E, 53 const CeedScalar Rd, const CeedScalar cv) { 54 CeedScalar u_sq = u[0]*u[0] + u[1]*u[1] + u[2]*u[2]; // Velocity square 55 // TODO Add in gravity's contribution 56 57 CeedScalar T = ( E / rho - u_sq / 2. ) / cv; 58 CeedScalar drdT = -rho / T; 59 CeedScalar drdP = 1. / ( Rd * T); 60 CeedScalar etot = E / rho ; 61 CeedScalar e2p = drdP * etot + 1. ; 62 CeedScalar e3p = ( E + rho * Rd * T ); 63 CeedScalar e4p = drdT * etot + rho * cv ; 64 65 for (CeedInt i=0; i<3; i++) { // Jacobian matrices for 3 directions 66 for (CeedInt j=0; j<3; j++) { // j counts F^{m_j} 67 // [row][col] of A_i 68 dF[i][j+1][0] = drdP * u[i] * u[j] + ((i==j) ? 1. : 0.); // F^{{m_j} wrt p 69 for (CeedInt k=0; k<3; k++) { // k counts the wrt vel_k 70 dF[i][0][k+1] = ((i==k) ? rho : 0.); // F^c wrt u_k 71 dF[i][j+1][k+1] = (((j==k) ? u[i] : 0.) + // F^m_j wrt u_k 72 ((i==k) ? u[j] : 0.) ) * rho; 73 dF[i][4][k+1] = rho * u[i] * u[k] 74 + ((i==k) ? e3p : 0.) ; // F^e wrt u_k 75 } 76 dF[i][j+1][4] = drdT * u[i] * u[j]; // F^{m_j} wrt T 77 } 78 dF[i][4][0] = u[i] * e2p; // F^e wrt p 79 dF[i][4][4] = u[i] * e4p; // F^e wrt T 80 dF[i][0][0] = u[i] * drdP; // F^c wrt p 81 dF[i][0][4] = u[i] * drdT; // F^c wrt T 82 } 83 } 84 85 CEED_QFUNCTION_HELPER void PrimitiveToConservative_fwd(const CeedScalar rho, 86 const CeedScalar u[3], const CeedScalar E, const CeedScalar Rd, 87 const CeedScalar cv, const CeedScalar dY[5], CeedScalar dU[5]) { 88 CeedScalar u_sq = u[0]*u[0] + u[1]*u[1] + u[2]*u[2]; 89 CeedScalar T = ( E / rho - u_sq / 2. ) / cv; 90 CeedScalar drdT = -rho / T; 91 CeedScalar drdP = 1. / ( Rd * T); 92 dU[0] = drdP * dY[0] + drdT * dY[4]; 93 CeedScalar de_kinetic = 0; 94 for (CeedInt i=0; i<3; i++) { 95 dU[1+i] = dU[0] * u[i] + rho * dY[1+i]; 96 de_kinetic += u[i] * dY[1+i]; 97 } 98 dU[4] = rho * cv * dY[4] + dU[0] * cv * T // internal energy: rho * e 99 + rho * de_kinetic + .5 * dU[0] * u_sq; // kinetic energy: .5 * rho * |u|^2 100 } 101 102 // ***************************************************************************** 103 // Helper function for computing Tau elements (stabilization constant) 104 // Model from: 105 // PHASTA 106 // 107 // Tau[i] = itau=0 which is diagonal-Shakib (3 values still but not spatial) 108 // 109 // Where NOT UPDATED YET 110 // ***************************************************************************** 111 CEED_QFUNCTION_HELPER void Tau_diagPrim(CeedScalar Tau_d[3], 112 const CeedScalar dXdx[3][3], const CeedScalar u[3], 113 const CeedScalar cv, const NewtonianIdealGasContext newt_ctx, 114 const CeedScalar mu, const CeedScalar dt, 115 const CeedScalar rho) { 116 // Context 117 const CeedScalar Ctau_t = newt_ctx->Ctau_t; 118 const CeedScalar Ctau_v = newt_ctx->Ctau_v; 119 const CeedScalar Ctau_C = newt_ctx->Ctau_C; 120 const CeedScalar Ctau_M = newt_ctx->Ctau_M; 121 const CeedScalar Ctau_E = newt_ctx->Ctau_E; 122 CeedScalar gijd[6]; 123 CeedScalar tau; 124 CeedScalar dts; 125 CeedScalar fact; 126 127 //*INDENT-OFF* 128 gijd[0] = dXdx[0][0] * dXdx[0][0] 129 + dXdx[1][0] * dXdx[1][0] 130 + dXdx[2][0] * dXdx[2][0]; 131 132 gijd[1] = dXdx[0][0] * dXdx[0][1] 133 + dXdx[1][0] * dXdx[1][1] 134 + dXdx[2][0] * dXdx[2][1]; 135 136 gijd[2] = dXdx[0][1] * dXdx[0][1] 137 + dXdx[1][1] * dXdx[1][1] 138 + dXdx[2][1] * dXdx[2][1]; 139 140 gijd[3] = dXdx[0][0] * dXdx[0][2] 141 + dXdx[1][0] * dXdx[1][2] 142 + dXdx[2][0] * dXdx[2][2]; 143 144 gijd[4] = dXdx[0][1] * dXdx[0][2] 145 + dXdx[1][1] * dXdx[1][2] 146 + dXdx[2][1] * dXdx[2][2]; 147 148 gijd[5] = dXdx[0][2] * dXdx[0][2] 149 + dXdx[1][2] * dXdx[1][2] 150 + dXdx[2][2] * dXdx[2][2]; 151 //*INDENT-ON* 152 153 dts = Ctau_t / dt ; 154 155 tau = rho*rho*((4. * dts * dts) 156 + u[0] * ( u[0] * gijd[0] + 2. * ( u[1] * gijd[1] + u[2] * gijd[3])) 157 + u[1] * ( u[1] * gijd[2] + 2. * u[2] * gijd[4]) 158 + u[2] * u[2] * gijd[5]) 159 + Ctau_v* mu * mu * 160 (gijd[0]*gijd[0] + gijd[2]*gijd[2] + gijd[5]*gijd[5] + 161 + 2. * (gijd[1]*gijd[1] + gijd[3]*gijd[3] + gijd[4]*gijd[4])); 162 163 fact=sqrt(tau); 164 165 Tau_d[0] = Ctau_C * fact / (rho*(gijd[0] + gijd[2] + gijd[5]))*0.125; 166 167 Tau_d[1] = Ctau_M / fact; 168 Tau_d[2] = Ctau_E / ( fact * cv ); 169 170 // consider putting back the way I initially had it Ctau_E * Tau_d[1] /cv 171 // to avoid a division if the compiler is smart enough to see that cv IS 172 // a constant that it could invert once for all elements 173 // but in that case energy tau is scaled by the product of Ctau_E * Ctau_M 174 // OR we could absorb cv into Ctau_E but this puts more burden on user to 175 // know how to change constants with a change of fluid or units. Same for 176 // Ctau_v * mu * mu IF AND ONLY IF we don't add viscosity law =f(T) 177 } 178 179 // ***************************************************************************** 180 // This QFunction sets a "still" initial condition for generic Newtonian IG problems 181 // ***************************************************************************** 182 CEED_QFUNCTION(ICsNewtonianIG)(void *ctx, CeedInt Q, 183 const CeedScalar *const *in, CeedScalar *const *out) { 184 // Inputs 185 const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 186 187 // Outputs 188 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 189 190 // Context 191 const SetupContext context = (SetupContext)ctx; 192 const CeedScalar theta0 = context->theta0; 193 const CeedScalar P0 = context->P0; 194 const CeedScalar cv = context->cv; 195 const CeedScalar cp = context->cp; 196 const CeedScalar *g = context->g; 197 const CeedScalar Rd = cp - cv; 198 199 // Quadrature Point Loop 200 CeedPragmaSIMD 201 for (CeedInt i=0; i<Q; i++) { 202 CeedScalar q[5] = {0.}; 203 204 // Setup 205 // -- Coordinates 206 const CeedScalar x[3] = {X[0][i], X[1][i], X[2][i]}; 207 const CeedScalar e_potential = -(g[0]*x[0] + g[1]*x[1] + g[2]*x[2]); 208 209 // -- Density 210 const CeedScalar rho = P0 / (Rd*theta0); 211 212 // Initial Conditions 213 q[0] = rho; 214 q[1] = 0.0; 215 q[2] = 0.0; 216 q[3] = 0.0; 217 q[4] = rho * (cv*theta0 + e_potential); 218 219 for (CeedInt j=0; j<5; j++) 220 q0[j][i] = q[j]; 221 } // End of Quadrature Point Loop 222 return 0; 223 } 224 225 // ***************************************************************************** 226 // This QFunction sets a "still" initial condition for generic Newtonian IG 227 // problems in primitive variables 228 // ***************************************************************************** 229 CEED_QFUNCTION(ICsNewtonianIG_Prim)(void *ctx, CeedInt Q, 230 const CeedScalar *const *in, CeedScalar *const *out) { 231 // Outputs 232 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 233 234 // Context 235 const SetupContext context = (SetupContext)ctx; 236 const CeedScalar theta0 = context->theta0; 237 const CeedScalar P0 = context->P0; 238 239 // Quadrature Point Loop 240 CeedPragmaSIMD 241 for (CeedInt i=0; i<Q; i++) { 242 CeedScalar q[5] = {0.}; 243 244 // Initial Conditions 245 q[0] = P0; 246 q[1] = 0.0; 247 q[2] = 0.0; 248 q[3] = 0.0; 249 q[4] = theta0; 250 251 for (CeedInt j=0; j<5; j++) 252 q0[j][i] = q[j]; 253 254 } // End of Quadrature Point Loop 255 return 0; 256 } 257 258 // ***************************************************************************** 259 // This QFunction implements the following formulation of Navier-Stokes with 260 // explicit time stepping method 261 // 262 // This is 3D compressible Navier-Stokes in conservation form with state 263 // variables of density, momentum density, and total energy density. 264 // 265 // State Variables: q = ( rho, U1, U2, U3, E ) 266 // rho - Mass Density 267 // Ui - Momentum Density, Ui = rho ui 268 // E - Total Energy Density, E = rho (cv T + (u u)/2 + g z) 269 // 270 // Navier-Stokes Equations: 271 // drho/dt + div( U ) = 0 272 // dU/dt + div( rho (u x u) + P I3 ) + rho g khat = div( Fu ) 273 // dE/dt + div( (E + P) u ) = div( Fe ) 274 // 275 // Viscous Stress: 276 // Fu = mu (grad( u ) + grad( u )^T + lambda div ( u ) I3) 277 // 278 // Thermal Stress: 279 // Fe = u Fu + k grad( T ) 280 // Equation of State 281 // P = (gamma - 1) (E - rho (u u) / 2 - rho g z) 282 // 283 // Stabilization: 284 // Tau = diag(TauC, TauM, TauM, TauM, TauE) 285 // f1 = rho sqrt(ui uj gij) 286 // gij = dXi/dX * dXi/dX 287 // TauC = Cc f1 / (8 gii) 288 // TauM = min( 1 , 1 / f1 ) 289 // TauE = TauM / (Ce cv) 290 // 291 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 292 // 293 // Constants: 294 // lambda = - 2 / 3, From Stokes hypothesis 295 // mu , Dynamic viscosity 296 // k , Thermal conductivity 297 // cv , Specific heat, constant volume 298 // cp , Specific heat, constant pressure 299 // g , Gravity 300 // gamma = cp / cv, Specific heat ratio 301 // 302 // We require the product of the inverse of the Jacobian (dXdx_j,k) and 303 // its transpose (dXdx_k,j) to properly compute integrals of the form: 304 // int( gradv gradu ) 305 // 306 // ***************************************************************************** 307 CEED_QFUNCTION(RHSFunction_Newtonian)(void *ctx, CeedInt Q, 308 const CeedScalar *const *in, CeedScalar *const *out) { 309 // *INDENT-OFF* 310 // Inputs 311 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 312 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 313 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 314 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 315 // Outputs 316 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 317 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 318 // *INDENT-ON* 319 320 // Context 321 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 322 const CeedScalar mu = context->mu; 323 const CeedScalar cv = context->cv; 324 const CeedScalar cp = context->cp; 325 const CeedScalar *g = context->g; 326 const CeedScalar dt = context->dt; 327 const CeedScalar gamma = cp / cv; 328 const CeedScalar Rd = cp - cv; 329 330 CeedPragmaSIMD 331 // Quadrature Point Loop 332 for (CeedInt i=0; i<Q; i++) { 333 CeedScalar U[5]; 334 for (int j=0; j<5; j++) U[j] = q[j][i]; 335 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 336 State s = StateFromU(context, U, x_i); 337 338 // -- Interp-to-Interp q_data 339 const CeedScalar wdetJ = q_data[0][i]; 340 // -- Interp-to-Grad q_data 341 // ---- Inverse of change of coordinate matrix: X_i,j 342 // *INDENT-OFF* 343 const CeedScalar dXdx[3][3] = {{q_data[1][i], 344 q_data[2][i], 345 q_data[3][i]}, 346 {q_data[4][i], 347 q_data[5][i], 348 q_data[6][i]}, 349 {q_data[7][i], 350 q_data[8][i], 351 q_data[9][i]} 352 }; 353 // *INDENT-ON* 354 State grad_s[3]; 355 for (CeedInt j=0; j<3; j++) { 356 CeedScalar dx_i[3] = {0}, dU[5]; 357 for (CeedInt k=0; k<5; k++) 358 dU[k] = Grad_q[0][k][i] * dXdx[0][j] + 359 Grad_q[1][k][i] * dXdx[1][j] + 360 Grad_q[2][k][i] * dXdx[2][j]; 361 dx_i[j] = 1.; 362 grad_s[j] = StateFromU_fwd(context, s, dU, x_i, dx_i); 363 } 364 365 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 366 KMStrainRate(grad_s, strain_rate); 367 NewtonianStress(context, strain_rate, kmstress); 368 KMUnpack(kmstress, stress); 369 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 370 371 StateConservative F_inviscid[3]; 372 FluxInviscid(context, s, F_inviscid); 373 374 // Total flux 375 CeedScalar Flux[5][3]; 376 for (CeedInt j=0; j<3; j++) { 377 Flux[0][j] = F_inviscid[j].density; 378 for (CeedInt k=0; k<3; k++) 379 Flux[k+1][j] = F_inviscid[j].momentum[k] - stress[k][j]; 380 Flux[4][j] = F_inviscid[j].E_total + Fe[j]; 381 } 382 383 for (CeedInt j=0; j<3; j++) { 384 for (CeedInt k=0; k<5; k++) { 385 Grad_v[j][k][i] = wdetJ * (dXdx[j][0] * Flux[k][0] + 386 dXdx[j][1] * Flux[k][1] + 387 dXdx[j][2] * Flux[k][2]); 388 } 389 } 390 391 const CeedScalar body_force[5] = {0, s.U.density *g[0], s.U.density *g[1], s.U.density *g[2], 0}; 392 for (int j=0; j<5; j++) 393 v[j][i] = wdetJ * body_force[j]; 394 395 // jacob_F_conv[3][5][5] = dF(convective)/dq at each direction 396 CeedScalar jacob_F_conv[3][5][5] = {0}; 397 computeFluxJacobian_NS(jacob_F_conv, s.U.density, s.Y.velocity, s.U.E_total, 398 gamma, g, x_i); 399 CeedScalar grad_U[5][3]; 400 for (CeedInt j=0; j<3; j++) { 401 grad_U[0][j] = grad_s[j].U.density; 402 for (CeedInt k=0; k<3; k++) grad_U[k+1][j] = grad_s[j].U.momentum[k]; 403 grad_U[4][j] = grad_s[j].U.E_total; 404 } 405 406 // strong_conv = dF/dq * dq/dx (Strong convection) 407 CeedScalar strong_conv[5] = {0}; 408 for (CeedInt j=0; j<3; j++) 409 for (CeedInt k=0; k<5; k++) 410 for (CeedInt l=0; l<5; l++) 411 strong_conv[k] += jacob_F_conv[j][k][l] * grad_U[l][j]; 412 413 // -- Stabilization method: none, SU, or SUPG 414 CeedScalar stab[5][3] = {{0.}}; 415 CeedScalar tau_strong_conv[5] = {0.}, tau_strong_conv_conservative[5] = {0}; 416 CeedScalar Tau_d[3] = {0.}; 417 switch (context->stabilization) { 418 case STAB_NONE: // Galerkin 419 break; 420 case STAB_SU: // SU 421 Tau_diagPrim(Tau_d, dXdx, s.Y.velocity, cv, context, mu, dt, s.U.density); 422 tau_strong_conv[0] = Tau_d[0] * strong_conv[0]; 423 tau_strong_conv[1] = Tau_d[1] * strong_conv[1]; 424 tau_strong_conv[2] = Tau_d[1] * strong_conv[2]; 425 tau_strong_conv[3] = Tau_d[1] * strong_conv[3]; 426 tau_strong_conv[4] = Tau_d[2] * strong_conv[4]; 427 PrimitiveToConservative_fwd(s.U.density, s.Y.velocity, s.U.E_total, Rd, cv, 428 tau_strong_conv, 429 tau_strong_conv_conservative); 430 for (CeedInt j=0; j<3; j++) 431 for (CeedInt k=0; k<5; k++) 432 for (CeedInt l=0; l<5; l++) 433 stab[k][j] += jacob_F_conv[j][k][l] * tau_strong_conv_conservative[l]; 434 435 for (CeedInt j=0; j<5; j++) 436 for (CeedInt k=0; k<3; k++) 437 Grad_v[k][j][i] -= wdetJ*(stab[j][0] * dXdx[k][0] + 438 stab[j][1] * dXdx[k][1] + 439 stab[j][2] * dXdx[k][2]); 440 break; 441 case STAB_SUPG: // SUPG is not implemented for explicit scheme 442 break; 443 } 444 445 } // End Quadrature Point Loop 446 447 // Return 448 return 0; 449 } 450 451 // ***************************************************************************** 452 // This QFunction implements the Navier-Stokes equations (mentioned above) with 453 // implicit time stepping method 454 // 455 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 456 // SUPG = Galerkin + grad(v) . ( Ai^T * Tau * (q_dot + Aj q,j - body force) ) 457 // (diffussive terms will be added later) 458 // 459 // ***************************************************************************** 460 CEED_QFUNCTION(IFunction_Newtonian)(void *ctx, CeedInt Q, 461 const CeedScalar *const *in, CeedScalar *const *out) { 462 // *INDENT-OFF* 463 // Inputs 464 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 465 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 466 (*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 467 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 468 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 469 // Outputs 470 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 471 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1], 472 (*jac_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[2]; 473 // *INDENT-ON* 474 // Context 475 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 476 const CeedScalar mu = context->mu; 477 const CeedScalar cv = context->cv; 478 const CeedScalar cp = context->cp; 479 const CeedScalar *g = context->g; 480 const CeedScalar dt = context->dt; 481 const CeedScalar gamma = cp / cv; 482 const CeedScalar Rd = cp - cv; 483 484 CeedPragmaSIMD 485 // Quadrature Point Loop 486 for (CeedInt i=0; i<Q; i++) { 487 CeedScalar U[5]; 488 for (CeedInt j=0; j<5; j++) U[j] = q[j][i]; 489 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 490 State s = StateFromU(context, U, x_i); 491 492 // -- Interp-to-Interp q_data 493 const CeedScalar wdetJ = q_data[0][i]; 494 // -- Interp-to-Grad q_data 495 // ---- Inverse of change of coordinate matrix: X_i,j 496 // *INDENT-OFF* 497 const CeedScalar dXdx[3][3] = {{q_data[1][i], 498 q_data[2][i], 499 q_data[3][i]}, 500 {q_data[4][i], 501 q_data[5][i], 502 q_data[6][i]}, 503 {q_data[7][i], 504 q_data[8][i], 505 q_data[9][i]} 506 }; 507 // *INDENT-ON* 508 State grad_s[3]; 509 for (CeedInt j=0; j<3; j++) { 510 CeedScalar dx_i[3] = {0}, dU[5]; 511 for (CeedInt k=0; k<5; k++) 512 dU[k] = Grad_q[0][k][i] * dXdx[0][j] + 513 Grad_q[1][k][i] * dXdx[1][j] + 514 Grad_q[2][k][i] * dXdx[2][j]; 515 dx_i[j] = 1.; 516 grad_s[j] = StateFromU_fwd(context, s, dU, x_i, dx_i); 517 } 518 519 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 520 KMStrainRate(grad_s, strain_rate); 521 NewtonianStress(context, strain_rate, kmstress); 522 KMUnpack(kmstress, stress); 523 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 524 525 StateConservative F_inviscid[3]; 526 FluxInviscid(context, s, F_inviscid); 527 528 529 // Total flux 530 CeedScalar Flux[5][3]; 531 for (CeedInt j=0; j<3; j++) { 532 Flux[0][j] = F_inviscid[j].density; 533 for (CeedInt k=0; k<3; k++) 534 Flux[k+1][j] = F_inviscid[j].momentum[k] - stress[k][j]; 535 Flux[4][j] = F_inviscid[j].E_total + Fe[j]; 536 } 537 538 for (CeedInt j=0; j<3; j++) { 539 for (CeedInt k=0; k<5; k++) { 540 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * Flux[k][0] + 541 dXdx[j][1] * Flux[k][1] + 542 dXdx[j][2] * Flux[k][2]); 543 } 544 } 545 546 const CeedScalar body_force[5] = {0, s.U.density *g[0], s.U.density *g[1], s.U.density *g[2], 0}; 547 for (CeedInt j=0; j<5; j++) 548 v[j][i] = wdetJ * (q_dot[j][i] - body_force[j]); 549 550 // jacob_F_conv[3][5][5] = dF(convective)/dq at each direction 551 CeedScalar jacob_F_conv[3][5][5] = {0}; 552 computeFluxJacobian_NS(jacob_F_conv, s.U.density, s.Y.velocity, s.U.E_total, 553 gamma, g, x_i); 554 CeedScalar grad_U[5][3]; 555 for (CeedInt j=0; j<3; j++) { 556 grad_U[0][j] = grad_s[j].U.density; 557 for (CeedInt k=0; k<3; k++) grad_U[k+1][j] = grad_s[j].U.momentum[k]; 558 grad_U[4][j] = grad_s[j].U.E_total; 559 } 560 561 // strong_conv = dF/dq * dq/dx (Strong convection) 562 CeedScalar strong_conv[5] = {0}; 563 for (CeedInt j=0; j<3; j++) 564 for (CeedInt k=0; k<5; k++) 565 for (CeedInt l=0; l<5; l++) 566 strong_conv[k] += jacob_F_conv[j][k][l] * grad_U[l][j]; 567 568 // Strong residual 569 CeedScalar strong_res[5]; 570 for (CeedInt j=0; j<5; j++) 571 strong_res[j] = q_dot[j][i] + strong_conv[j] - body_force[j]; 572 573 // -- Stabilization method: none, SU, or SUPG 574 CeedScalar stab[5][3] = {{0.}}; 575 CeedScalar tau_strong_res[5] = {0.}, tau_strong_res_conservative[5] = {0}; 576 CeedScalar tau_strong_conv[5] = {0.}, tau_strong_conv_conservative[5] = {0}; 577 CeedScalar Tau_d[3] = {0.}; 578 switch (context->stabilization) { 579 case STAB_NONE: // Galerkin 580 break; 581 case STAB_SU: // SU 582 Tau_diagPrim(Tau_d, dXdx, s.Y.velocity, cv, context, mu, dt, s.U.density); 583 tau_strong_conv[0] = Tau_d[0] * strong_conv[0]; 584 tau_strong_conv[1] = Tau_d[1] * strong_conv[1]; 585 tau_strong_conv[2] = Tau_d[1] * strong_conv[2]; 586 tau_strong_conv[3] = Tau_d[1] * strong_conv[3]; 587 tau_strong_conv[4] = Tau_d[2] * strong_conv[4]; 588 PrimitiveToConservative_fwd(s.U.density, s.Y.velocity, s.U.E_total, Rd, cv, 589 tau_strong_conv, tau_strong_conv_conservative); 590 for (CeedInt j=0; j<3; j++) 591 for (CeedInt k=0; k<5; k++) 592 for (CeedInt l=0; l<5; l++) 593 stab[k][j] += jacob_F_conv[j][k][l] * tau_strong_conv_conservative[l]; 594 595 for (CeedInt j=0; j<5; j++) 596 for (CeedInt k=0; k<3; k++) 597 Grad_v[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 598 stab[j][1] * dXdx[k][1] + 599 stab[j][2] * dXdx[k][2]); 600 601 break; 602 case STAB_SUPG: // SUPG 603 Tau_diagPrim(Tau_d, dXdx, s.Y.velocity, cv, context, mu, dt, s.U.density); 604 tau_strong_res[0] = Tau_d[0] * strong_res[0]; 605 tau_strong_res[1] = Tau_d[1] * strong_res[1]; 606 tau_strong_res[2] = Tau_d[1] * strong_res[2]; 607 tau_strong_res[3] = Tau_d[1] * strong_res[3]; 608 tau_strong_res[4] = Tau_d[2] * strong_res[4]; 609 610 PrimitiveToConservative_fwd(s.U.density, s.Y.velocity, s.U.E_total, Rd, cv, 611 tau_strong_res, tau_strong_res_conservative); 612 for (CeedInt j=0; j<3; j++) 613 for (CeedInt k=0; k<5; k++) 614 for (CeedInt l=0; l<5; l++) 615 stab[k][j] += jacob_F_conv[j][k][l] * tau_strong_res_conservative[l]; 616 617 for (CeedInt j=0; j<5; j++) 618 for (CeedInt k=0; k<3; k++) 619 Grad_v[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 620 stab[j][1] * dXdx[k][1] + 621 stab[j][2] * dXdx[k][2]); 622 break; 623 } 624 for (CeedInt j=0; j<5; j++) jac_data[j][i] = U[j]; 625 for (CeedInt j=0; j<6; j++) jac_data[5+j][i] = kmstress[j]; 626 for (CeedInt j=0; j<3; j++) jac_data[5+6+j][i] = Tau_d[j]; 627 628 } // End Quadrature Point Loop 629 630 // Return 631 return 0; 632 } 633 634 // ***************************************************************************** 635 // This QFunction implements the jacobean of the Navier-Stokes equations 636 // for implicit time stepping method. 637 // 638 // ***************************************************************************** 639 CEED_QFUNCTION(IJacobian_Newtonian)(void *ctx, CeedInt Q, 640 const CeedScalar *const *in, 641 CeedScalar *const *out) { 642 // *INDENT-OFF* 643 // Inputs 644 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 645 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 646 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 647 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 648 (*jac_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 649 // Outputs 650 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 651 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 652 // *INDENT-ON* 653 // Context 654 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 655 const CeedScalar *g = context->g; 656 const CeedScalar cp = context->cp; 657 const CeedScalar cv = context->cv; 658 const CeedScalar Rd = cp - cv; 659 const CeedScalar gamma = cp / cv; 660 661 CeedPragmaSIMD 662 // Quadrature Point Loop 663 for (CeedInt i=0; i<Q; i++) { 664 // -- Interp-to-Interp q_data 665 const CeedScalar wdetJ = q_data[0][i]; 666 // -- Interp-to-Grad q_data 667 // ---- Inverse of change of coordinate matrix: X_i,j 668 // *INDENT-OFF* 669 const CeedScalar dXdx[3][3] = {{q_data[1][i], 670 q_data[2][i], 671 q_data[3][i]}, 672 {q_data[4][i], 673 q_data[5][i], 674 q_data[6][i]}, 675 {q_data[7][i], 676 q_data[8][i], 677 q_data[9][i]} 678 }; 679 // *INDENT-ON* 680 681 CeedScalar U[5], kmstress[6], Tau_d[3] __attribute((unused)); 682 for (int j=0; j<5; j++) U[j] = jac_data[j][i]; 683 for (int j=0; j<6; j++) kmstress[j] = jac_data[5+j][i]; 684 for (int j=0; j<3; j++) Tau_d[j] = jac_data[5+6+j][i]; 685 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 686 State s = StateFromU(context, U, x_i); 687 688 CeedScalar dU[5], dx0[3] = {0}; 689 for (int j=0; j<5; j++) dU[j] = dq[j][i]; 690 State ds = StateFromU_fwd(context, s, dU, x_i, dx0); 691 692 State grad_ds[3]; 693 for (int j=0; j<3; j++) { 694 CeedScalar dUj[5]; 695 for (int k=0; k<5; k++) dUj[k] = Grad_dq[0][k][i] * dXdx[0][j] 696 + Grad_dq[1][k][i] * dXdx[1][j] 697 + Grad_dq[2][k][i] * dXdx[2][j]; 698 grad_ds[j] = StateFromU_fwd(context, s, dUj, x_i, dx0); 699 } 700 701 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 702 KMStrainRate(grad_ds, dstrain_rate); 703 NewtonianStress(context, dstrain_rate, dkmstress); 704 KMUnpack(dkmstress, dstress); 705 KMUnpack(kmstress, stress); 706 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 707 708 StateConservative dF_inviscid[3]; 709 FluxInviscid_fwd(context, s, ds, dF_inviscid); 710 711 // Total flux 712 CeedScalar dFlux[5][3]; 713 for (int j=0; j<3; j++) { 714 dFlux[0][j] = dF_inviscid[j].density; 715 for (int k=0; k<3; k++) 716 dFlux[k+1][j] = dF_inviscid[j].momentum[k] - dstress[k][j]; 717 dFlux[4][j] = dF_inviscid[j].E_total + dFe[j]; 718 } 719 720 for (int j=0; j<3; j++) { 721 for (int k=0; k<5; k++) { 722 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * dFlux[k][0] + 723 dXdx[j][1] * dFlux[k][1] + 724 dXdx[j][2] * dFlux[k][2]); 725 } 726 } 727 728 const CeedScalar dbody_force[5] = {0, ds.U.density *g[0], ds.U.density *g[1], ds.U.density *g[2], 0}; 729 for (int j=0; j<5; j++) 730 v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]); 731 732 if (1) { 733 CeedScalar jacob_F_conv[3][5][5] = {0}; 734 computeFluxJacobian_NS(jacob_F_conv, s.U.density, s.Y.velocity, s.U.E_total, 735 gamma, g, x_i); 736 CeedScalar grad_dU[5][3]; 737 for (int j=0; j<3; j++) { 738 grad_dU[0][j] = grad_ds[j].U.density; 739 for (int k=0; k<3; k++) grad_dU[k+1][j] = grad_ds[j].U.momentum[k]; 740 grad_dU[4][j] = grad_ds[j].U.E_total; 741 } 742 CeedScalar dstrong_conv[5] = {0}; 743 for (int j=0; j<3; j++) 744 for (int k=0; k<5; k++) 745 for (int l=0; l<5; l++) 746 dstrong_conv[k] += jacob_F_conv[j][k][l] * grad_dU[l][j]; 747 CeedScalar dstrong_res[5]; 748 for (int j=0; j<5; j++) 749 dstrong_res[j] = context->ijacobian_time_shift * dU[j] + dstrong_conv[j] - 750 dbody_force[j]; 751 CeedScalar dtau_strong_res[5] = {0.}, dtau_strong_res_conservative[5] = {0}; 752 dtau_strong_res[0] = Tau_d[0] * dstrong_res[0]; 753 dtau_strong_res[1] = Tau_d[1] * dstrong_res[1]; 754 dtau_strong_res[2] = Tau_d[1] * dstrong_res[2]; 755 dtau_strong_res[3] = Tau_d[1] * dstrong_res[3]; 756 dtau_strong_res[4] = Tau_d[2] * dstrong_res[4]; 757 PrimitiveToConservative_fwd(s.U.density, s.Y.velocity, s.U.E_total, Rd, cv, 758 dtau_strong_res, dtau_strong_res_conservative); 759 CeedScalar dstab[5][3] = {0}; 760 for (int j=0; j<3; j++) 761 for (int k=0; k<5; k++) 762 for (int l=0; l<5; l++) 763 dstab[k][j] += jacob_F_conv[j][k][l] * dtau_strong_res_conservative[l]; 764 for (int j=0; j<5; j++) 765 for (int k=0; k<3; k++) 766 Grad_v[k][j][i] += wdetJ*(dstab[j][0] * dXdx[k][0] + 767 dstab[j][1] * dXdx[k][1] + 768 dstab[j][2] * dXdx[k][2]); 769 770 } 771 } // End Quadrature Point Loop 772 return 0; 773 } 774 775 // Compute boundary integral (ie. for strongly set inflows) 776 CEED_QFUNCTION(BoundaryIntegral)(void *ctx, CeedInt Q, 777 const CeedScalar *const *in, 778 CeedScalar *const *out) { 779 780 //*INDENT-OFF* 781 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 782 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 783 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 784 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 785 786 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA]) out[0], 787 (*jac_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA]) out[1]; 788 789 //*INDENT-ON* 790 791 const NewtonianIdealGasContext context = (NewtonianIdealGasContext) ctx; 792 const bool is_implicit = context->is_implicit; 793 794 CeedPragmaSIMD 795 for(CeedInt i=0; i<Q; i++) { 796 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 797 const CeedScalar solution_state[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 798 State s = context->is_primitive ? StateFromY(context, solution_state, x_i) 799 : StateFromU(context, solution_state, x_i); 800 801 const CeedScalar wdetJb = (is_implicit ? -1. : 1.) * q_data_sur[0][i]; 802 // ---- Normal vect 803 const CeedScalar norm[3] = {q_data_sur[1][i], 804 q_data_sur[2][i], 805 q_data_sur[3][i] 806 }; 807 808 const CeedScalar dXdx[2][3] = { 809 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 810 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 811 }; 812 813 State grad_s[3]; 814 for (CeedInt j=0; j<3; j++) { 815 CeedScalar dx_i[3] = {0}, dU[5]; 816 for (CeedInt k=0; k<5; k++) 817 dU[k] = Grad_q[0][k][i] * dXdx[0][j] + 818 Grad_q[1][k][i] * dXdx[1][j]; 819 dx_i[j] = 1.; 820 grad_s[j] = StateFromU_fwd(context, s, dU, x_i, dx_i); 821 } 822 823 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 824 KMStrainRate(grad_s, strain_rate); 825 NewtonianStress(context, strain_rate, kmstress); 826 KMUnpack(kmstress, stress); 827 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 828 829 StateConservative F_inviscid[3]; 830 FluxInviscid(context, s, F_inviscid); 831 832 CeedScalar Flux[5] = {0.}; 833 for (int j=0; j<3; j++) { 834 Flux[0] += F_inviscid[j].density * norm[j]; 835 for (int k=0; k<3; k++) 836 Flux[k+1] += (F_inviscid[j].momentum[k] - stress[k][j]) * norm[j]; 837 Flux[4] += (F_inviscid[j].E_total + Fe[j])*norm[j]; 838 } 839 840 // -- Density 841 v[0][i] = -wdetJb * Flux[0]; 842 843 // -- Momentum 844 for (CeedInt j=0; j<3; j++) 845 v[j+1][i] = -wdetJb * Flux[j+1]; 846 847 // -- Total Energy Density 848 v[4][i] = -wdetJb * Flux[4]; 849 850 if (context->is_primitive) { 851 jac_data_sur[0][i] = s.Y.pressure; 852 for (int j=0; j<3; j++) jac_data_sur[j+1][i] = s.Y.velocity[j]; 853 jac_data_sur[4][i] = s.Y.temperature; 854 } else { 855 jac_data_sur[0][i] = s.U.density; 856 for (int j=0; j<3; j++) jac_data_sur[j+1][i] = s.U.momentum[j]; 857 jac_data_sur[4][i] = s.U.E_total; 858 } 859 for (int j=0; j<6; j++) jac_data_sur[5+j][i] = kmstress[j]; 860 } 861 return 0; 862 } 863 864 // Jacobian for "set nothing" boundary integral 865 CEED_QFUNCTION(BoundaryIntegral_Jacobian)(void *ctx, CeedInt Q, 866 const CeedScalar *const *in, 867 CeedScalar *const *out) { 868 // *INDENT-OFF* 869 // Inputs 870 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 871 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 872 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 873 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 874 (*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 875 // Outputs 876 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 877 // *INDENT-ON* 878 879 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 880 const bool implicit = context->is_implicit; 881 882 CeedPragmaSIMD 883 // Quadrature Point Loop 884 for (CeedInt i=0; i<Q; i++) { 885 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 886 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 887 const CeedScalar norm[3] = {q_data_sur[1][i], 888 q_data_sur[2][i], 889 q_data_sur[3][i] 890 }; 891 const CeedScalar dXdx[2][3] = { 892 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 893 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 894 }; 895 896 CeedScalar state[5], kmstress[6], dstate[5], dx_i[3] = {0.}; 897 for (int j=0; j<5; j++) state[j] = jac_data_sur[j][i]; 898 for (int j=0; j<6; j++) kmstress[j] = jac_data_sur[5+j][i]; 899 for (int j=0; j<5; j++) dstate[j] = dq[j][i]; 900 901 State s, ds; 902 if (context->is_primitive) { 903 s = StateFromY(context, state, x_i); 904 ds = StateFromY_fwd(context, s, dstate, x_i, dx_i); 905 } else { 906 s = StateFromU(context, state, x_i); 907 ds = StateFromU_fwd(context, s, dstate, x_i, dx_i); 908 } 909 910 State grad_ds[3]; 911 for (CeedInt j=0; j<3; j++) { 912 CeedScalar dx_i[3] = {0}, dUj[5]; 913 for (CeedInt k=0; k<5; k++) 914 dUj[k] = Grad_dq[0][k][i] * dXdx[0][j] + 915 Grad_dq[1][k][i] * dXdx[1][j]; 916 dx_i[j] = 1.; 917 grad_ds[j] = StateFromU_fwd(context, s, dUj, x_i, dx_i); 918 } 919 920 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 921 KMStrainRate(grad_ds, dstrain_rate); 922 NewtonianStress(context, dstrain_rate, dkmstress); 923 KMUnpack(dkmstress, dstress); 924 KMUnpack(kmstress, stress); 925 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 926 927 StateConservative dF_inviscid[3]; 928 FluxInviscid_fwd(context, s, ds, dF_inviscid); 929 930 CeedScalar dFlux[5] = {0.}; 931 for (int j=0; j<3; j++) { 932 dFlux[0] += dF_inviscid[j].density * norm[j]; 933 for (int k=0; k<3; k++) 934 dFlux[k+1] += (dF_inviscid[j].momentum[k] - dstress[k][j]) * norm[j]; 935 dFlux[4] += (dF_inviscid[j].E_total + dFe[j]) * norm[j]; 936 } 937 938 for (int j=0; j<5; j++) 939 v[j][i] = -wdetJb * dFlux[j]; 940 } // End Quadrature Point Loop 941 return 0; 942 } 943 944 // Outflow boundary condition, weakly setting a constant pressure 945 CEED_QFUNCTION(PressureOutflow)(void *ctx, CeedInt Q, 946 const CeedScalar *const *in, 947 CeedScalar *const *out) { 948 // *INDENT-OFF* 949 // Inputs 950 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 951 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 952 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 953 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 954 // Outputs 955 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 956 (*jac_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[1]; 957 // *INDENT-ON* 958 959 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 960 const bool implicit = context->is_implicit; 961 const CeedScalar P0 = context->P0; 962 963 CeedPragmaSIMD 964 // Quadrature Point Loop 965 for (CeedInt i=0; i<Q; i++) { 966 // Setup 967 // -- Interp in 968 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 969 const CeedScalar solution_state[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 970 State s = context->is_primitive ? StateFromY(context, solution_state, x_i) 971 : StateFromU(context, solution_state, x_i); 972 s.Y.pressure = P0; 973 974 // -- Interp-to-Interp q_data 975 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 976 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 977 // We can effect this by swapping the sign on this weight 978 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 979 980 // ---- Normal vect 981 const CeedScalar norm[3] = {q_data_sur[1][i], 982 q_data_sur[2][i], 983 q_data_sur[3][i] 984 }; 985 986 const CeedScalar dXdx[2][3] = { 987 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 988 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 989 }; 990 991 State grad_s[3]; 992 for (CeedInt j=0; j<3; j++) { 993 CeedScalar dx_i[3] = {0}, dU[5]; 994 for (CeedInt k=0; k<5; k++) 995 dU[k] = Grad_q[0][k][i] * dXdx[0][j] + 996 Grad_q[1][k][i] * dXdx[1][j]; 997 dx_i[j] = 1.; 998 grad_s[j] = StateFromU_fwd(context, s, dU, x_i, dx_i); 999 } 1000 1001 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 1002 KMStrainRate(grad_s, strain_rate); 1003 NewtonianStress(context, strain_rate, kmstress); 1004 KMUnpack(kmstress, stress); 1005 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 1006 1007 StateConservative F_inviscid[3]; 1008 FluxInviscid(context, s, F_inviscid); 1009 1010 CeedScalar Flux[5] = {0.}; 1011 for (int j=0; j<3; j++) { 1012 Flux[0] += F_inviscid[j].density * norm[j]; 1013 for (int k=0; k<3; k++) 1014 Flux[k+1] += (F_inviscid[j].momentum[k] - stress[k][j]) * norm[j]; 1015 Flux[4] += (F_inviscid[j].E_total + Fe[j])*norm[j]; 1016 } 1017 1018 // -- Density 1019 v[0][i] = -wdetJb * Flux[0]; 1020 1021 // -- Momentum 1022 for (CeedInt j=0; j<3; j++) 1023 v[j+1][i] = -wdetJb * Flux[j+1]; 1024 1025 // -- Total Energy Density 1026 v[4][i] = -wdetJb * Flux[4]; 1027 1028 // Save values for Jacobian 1029 if (context->is_primitive) { 1030 jac_data_sur[0][i] = s.Y.pressure; 1031 for (int j=0; j<3; j++) jac_data_sur[j+1][i] = s.Y.velocity[j]; 1032 jac_data_sur[4][i] = s.Y.temperature; 1033 } else { 1034 jac_data_sur[0][i] = s.U.density; 1035 for (int j=0; j<3; j++) jac_data_sur[j+1][i] = s.U.momentum[j]; 1036 jac_data_sur[4][i] = s.U.E_total; 1037 } 1038 for (int j=0; j<6; j++) jac_data_sur[5+j][i] = kmstress[j]; 1039 } // End Quadrature Point Loop 1040 return 0; 1041 } 1042 1043 // Jacobian for weak-pressure outflow boundary condition 1044 CEED_QFUNCTION(PressureOutflow_Jacobian)(void *ctx, CeedInt Q, 1045 const CeedScalar *const *in, 1046 CeedScalar *const *out) { 1047 // *INDENT-OFF* 1048 // Inputs 1049 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 1050 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 1051 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 1052 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 1053 (*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 1054 // Outputs 1055 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 1056 // *INDENT-ON* 1057 1058 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 1059 const bool implicit = context->is_implicit; 1060 1061 CeedPragmaSIMD 1062 // Quadrature Point Loop 1063 for (CeedInt i=0; i<Q; i++) { 1064 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 1065 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 1066 const CeedScalar norm[3] = {q_data_sur[1][i], 1067 q_data_sur[2][i], 1068 q_data_sur[3][i] 1069 }; 1070 const CeedScalar dXdx[2][3] = { 1071 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 1072 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 1073 }; 1074 1075 CeedScalar state[5], kmstress[6], dstate[5], dx_i[3] = {0.}; 1076 for (int j=0; j<5; j++) state[j] = jac_data_sur[j][i]; 1077 for (int j=0; j<6; j++) kmstress[j] = jac_data_sur[5+j][i]; 1078 for (int j=0; j<5; j++) dstate[j] = dq[j][i]; 1079 1080 State s, ds; 1081 if (context->is_primitive) { 1082 s = StateFromY(context, state, x_i); 1083 ds = StateFromY_fwd(context, s, dstate, x_i, dx_i); 1084 } else { 1085 s = StateFromU(context, state, x_i); 1086 ds = StateFromU_fwd(context, s, dstate, x_i, dx_i); 1087 } 1088 s.Y.pressure = context->P0; 1089 ds.Y.pressure = 0.; 1090 1091 State grad_ds[3]; 1092 for (CeedInt j=0; j<3; j++) { 1093 CeedScalar dx_i[3] = {0}, dUj[5]; 1094 for (CeedInt k=0; k<5; k++) 1095 dUj[k] = Grad_dq[0][k][i] * dXdx[0][j] + 1096 Grad_dq[1][k][i] * dXdx[1][j]; 1097 dx_i[j] = 1.; 1098 grad_ds[j] = StateFromU_fwd(context, s, dUj, x_i, dx_i); 1099 } 1100 1101 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 1102 KMStrainRate(grad_ds, dstrain_rate); 1103 NewtonianStress(context, dstrain_rate, dkmstress); 1104 KMUnpack(dkmstress, dstress); 1105 KMUnpack(kmstress, stress); 1106 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 1107 1108 StateConservative dF_inviscid[3]; 1109 FluxInviscid_fwd(context, s, ds, dF_inviscid); 1110 1111 CeedScalar dFlux[5] = {0.}; 1112 for (int j=0; j<3; j++) { 1113 dFlux[0] += dF_inviscid[j].density * norm[j]; 1114 for (int k=0; k<3; k++) 1115 dFlux[k+1] += (dF_inviscid[j].momentum[k] - dstress[k][j]) * norm[j]; 1116 dFlux[4] += (dF_inviscid[j].E_total + dFe[j]) * norm[j]; 1117 } 1118 1119 for (int j=0; j<5; j++) 1120 v[j][i] = -wdetJb * dFlux[j]; 1121 } // End Quadrature Point Loop 1122 return 0; 1123 } 1124 1125 // ***************************************************************************** 1126 // This QFunction implements the Navier-Stokes equations (mentioned above) in 1127 // primitive variables and with implicit time stepping method 1128 // 1129 // ***************************************************************************** 1130 CEED_QFUNCTION(IFunction_Newtonian_Prim)(void *ctx, CeedInt Q, 1131 const CeedScalar *const *in, CeedScalar *const *out) { 1132 // *INDENT-OFF* 1133 // Inputs 1134 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 1135 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 1136 (*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 1137 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 1138 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 1139 // Outputs 1140 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 1141 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1], 1142 (*jac_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[2]; 1143 // *INDENT-ON* 1144 // Context 1145 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 1146 const CeedScalar mu = context->mu; 1147 const CeedScalar cv = context->cv; 1148 const CeedScalar cp = context->cp; 1149 const CeedScalar *g = context->g; 1150 const CeedScalar dt = context->dt; 1151 const CeedScalar gamma = cp / cv; 1152 const CeedScalar Rd = cp - cv; 1153 1154 CeedPragmaSIMD 1155 // Quadrature Point Loop 1156 for (CeedInt i=0; i<Q; i++) { 1157 CeedScalar Y[5]; 1158 for (CeedInt j=0; j<5; j++) Y[j] = q[j][i]; 1159 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 1160 State s = StateFromY(context, Y, x_i); 1161 1162 // -- Interp-to-Interp q_data 1163 const CeedScalar wdetJ = q_data[0][i]; 1164 // -- Interp-to-Grad q_data 1165 // ---- Inverse of change of coordinate matrix: X_i,j 1166 // *INDENT-OFF* 1167 const CeedScalar dXdx[3][3] = {{q_data[1][i], 1168 q_data[2][i], 1169 q_data[3][i]}, 1170 {q_data[4][i], 1171 q_data[5][i], 1172 q_data[6][i]}, 1173 {q_data[7][i], 1174 q_data[8][i], 1175 q_data[9][i]} 1176 }; 1177 // *INDENT-ON* 1178 State grad_s[3]; 1179 for (CeedInt j=0; j<3; j++) { 1180 CeedScalar dx_i[3] = {0}, dY[5]; 1181 for (CeedInt k=0; k<5; k++) 1182 dY[k] = Grad_q[0][k][i] * dXdx[0][j] + 1183 Grad_q[1][k][i] * dXdx[1][j] + 1184 Grad_q[2][k][i] * dXdx[2][j]; 1185 dx_i[j] = 1.; 1186 grad_s[j] = StateFromY_fwd(context, s, dY, x_i, dx_i); 1187 } 1188 1189 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 1190 KMStrainRate(grad_s, strain_rate); 1191 NewtonianStress(context, strain_rate, kmstress); 1192 KMUnpack(kmstress, stress); 1193 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 1194 1195 StateConservative F_inviscid[3]; 1196 FluxInviscid(context, s, F_inviscid); 1197 1198 // Total flux 1199 CeedScalar Flux[5][3]; 1200 for (CeedInt j=0; j<3; j++) { 1201 Flux[0][j] = F_inviscid[j].density; 1202 for (CeedInt k=0; k<3; k++) 1203 Flux[k+1][j] = F_inviscid[j].momentum[k] - stress[k][j]; 1204 Flux[4][j] = F_inviscid[j].E_total + Fe[j]; 1205 } 1206 1207 for (CeedInt j=0; j<3; j++) { 1208 for (CeedInt k=0; k<5; k++) { 1209 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * Flux[k][0] + 1210 dXdx[j][1] * Flux[k][1] + 1211 dXdx[j][2] * Flux[k][2]); 1212 } 1213 } 1214 1215 const CeedScalar body_force[5] = {0, s.U.density *g[0], s.U.density *g[1], s.U.density *g[2], 0}; 1216 1217 CeedScalar Y_dot[5], dx0[3] = {0}; 1218 for (int j=0; j<5; j++) Y_dot[j] = q_dot[j][i]; 1219 State s_dot = StateFromY_fwd(context, s, Y_dot, x_i, dx0); 1220 1221 CeedScalar U_dot[5] = {0.}; 1222 U_dot[0] = s_dot.U.density; 1223 for (CeedInt j=0; j<3; j++) 1224 U_dot[j+1] = s_dot.U.momentum[j]; 1225 U_dot[4] = s_dot.U.E_total; 1226 1227 for (CeedInt j=0; j<5; j++) 1228 v[j][i] = wdetJ * (U_dot[j] - body_force[j]); 1229 1230 // jacob_F_conv[3][5][5] = dF(convective)/dq at each direction 1231 CeedScalar jacob_F_conv[3][5][5] = {0}; 1232 computeFluxJacobian_NS(jacob_F_conv, s.U.density, s.Y.velocity, s.U.E_total, 1233 gamma, g, x_i); 1234 CeedScalar grad_U[5][3]; 1235 for (CeedInt j=0; j<3; j++) { 1236 grad_U[0][j] = grad_s[j].U.density; 1237 for (CeedInt k=0; k<3; k++) grad_U[k+1][j] = grad_s[j].U.momentum[k]; 1238 grad_U[4][j] = grad_s[j].U.E_total; 1239 } 1240 1241 // strong_conv = dF/dq * dq/dx (Strong convection) 1242 CeedScalar strong_conv[5] = {0}; 1243 for (CeedInt j=0; j<3; j++) 1244 for (CeedInt k=0; k<5; k++) 1245 for (CeedInt l=0; l<5; l++) 1246 strong_conv[k] += jacob_F_conv[j][k][l] * grad_U[l][j]; 1247 1248 // Strong residual 1249 CeedScalar strong_res[5]; 1250 for (CeedInt j=0; j<5; j++) 1251 strong_res[j] = U_dot[j] + strong_conv[j] - body_force[j]; 1252 1253 // -- Stabilization method: none, SU, or SUPG 1254 CeedScalar stab[5][3] = {{0.}}; 1255 CeedScalar tau_strong_res[5] = {0.}, tau_strong_res_conservative[5] = {0}; 1256 CeedScalar tau_strong_conv[5] = {0.}, tau_strong_conv_conservative[5] = {0}; 1257 CeedScalar Tau_d[3] = {0.}; 1258 switch (context->stabilization) { 1259 case STAB_NONE: // Galerkin 1260 break; 1261 case STAB_SU: // SU 1262 Tau_diagPrim(Tau_d, dXdx, s.Y.velocity, cv, context, mu, dt, s.U.density); 1263 tau_strong_conv[0] = Tau_d[0] * strong_conv[0]; 1264 tau_strong_conv[1] = Tau_d[1] * strong_conv[1]; 1265 tau_strong_conv[2] = Tau_d[1] * strong_conv[2]; 1266 tau_strong_conv[3] = Tau_d[1] * strong_conv[3]; 1267 tau_strong_conv[4] = Tau_d[2] * strong_conv[4]; 1268 PrimitiveToConservative_fwd(s.U.density, s.Y.velocity, s.U.E_total, Rd, cv, 1269 tau_strong_conv, tau_strong_conv_conservative); 1270 for (CeedInt j=0; j<3; j++) 1271 for (CeedInt k=0; k<5; k++) 1272 for (CeedInt l=0; l<5; l++) 1273 stab[k][j] += jacob_F_conv[j][k][l] * tau_strong_conv_conservative[l]; 1274 1275 for (CeedInt j=0; j<5; j++) 1276 for (CeedInt k=0; k<3; k++) 1277 Grad_v[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 1278 stab[j][1] * dXdx[k][1] + 1279 stab[j][2] * dXdx[k][2]); 1280 1281 break; 1282 case STAB_SUPG: // SUPG 1283 Tau_diagPrim(Tau_d, dXdx, s.Y.velocity, cv, context, mu, dt, s.U.density); 1284 tau_strong_res[0] = Tau_d[0] * strong_res[0]; 1285 tau_strong_res[1] = Tau_d[1] * strong_res[1]; 1286 tau_strong_res[2] = Tau_d[1] * strong_res[2]; 1287 tau_strong_res[3] = Tau_d[1] * strong_res[3]; 1288 tau_strong_res[4] = Tau_d[2] * strong_res[4]; 1289 1290 PrimitiveToConservative_fwd(s.U.density, s.Y.velocity, s.U.E_total, Rd, cv, 1291 tau_strong_res, tau_strong_res_conservative); 1292 for (CeedInt j=0; j<3; j++) 1293 for (CeedInt k=0; k<5; k++) 1294 for (CeedInt l=0; l<5; l++) 1295 stab[k][j] += jacob_F_conv[j][k][l] * tau_strong_res_conservative[l]; 1296 1297 for (CeedInt j=0; j<5; j++) 1298 for (CeedInt k=0; k<3; k++) 1299 Grad_v[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 1300 stab[j][1] * dXdx[k][1] + 1301 stab[j][2] * dXdx[k][2]); 1302 break; 1303 } 1304 for (CeedInt j=0; j<5; j++) jac_data[j][i] = Y[j]; 1305 for (CeedInt j=0; j<6; j++) jac_data[5+j][i] = kmstress[j]; 1306 for (CeedInt j=0; j<3; j++) jac_data[5+6+j][i] = Tau_d[j]; 1307 1308 } // End Quadrature Point Loop 1309 1310 // Return 1311 return 0; 1312 } 1313 1314 // ***************************************************************************** 1315 // This QFunction implements the jacobean of the Navier-Stokes equations 1316 // in primitive variables for implicit time stepping method. 1317 // 1318 // ***************************************************************************** 1319 CEED_QFUNCTION(IJacobian_Newtonian_Prim)(void *ctx, CeedInt Q, 1320 const CeedScalar *const *in, CeedScalar *const *out) { 1321 // *INDENT-OFF* 1322 // Inputs 1323 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 1324 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 1325 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 1326 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 1327 (*jac_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 1328 // Outputs 1329 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 1330 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 1331 // *INDENT-ON* 1332 // Context 1333 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 1334 const CeedScalar *g = context->g; 1335 const CeedScalar cp = context->cp; 1336 const CeedScalar cv = context->cv; 1337 const CeedScalar Rd = cp - cv; 1338 const CeedScalar gamma = cp / cv; 1339 1340 CeedPragmaSIMD 1341 // Quadrature Point Loop 1342 for (CeedInt i=0; i<Q; i++) { 1343 // -- Interp-to-Interp q_data 1344 const CeedScalar wdetJ = q_data[0][i]; 1345 // -- Interp-to-Grad q_data 1346 // ---- Inverse of change of coordinate matrix: X_i,j 1347 // *INDENT-OFF* 1348 const CeedScalar dXdx[3][3] = {{q_data[1][i], 1349 q_data[2][i], 1350 q_data[3][i]}, 1351 {q_data[4][i], 1352 q_data[5][i], 1353 q_data[6][i]}, 1354 {q_data[7][i], 1355 q_data[8][i], 1356 q_data[9][i]} 1357 }; 1358 // *INDENT-ON* 1359 1360 CeedScalar Y[5], kmstress[6], Tau_d[3] __attribute((unused)); 1361 for (int j=0; j<5; j++) Y[j] = jac_data[j][i]; 1362 for (int j=0; j<6; j++) kmstress[j] = jac_data[5+j][i]; 1363 for (int j=0; j<3; j++) Tau_d[j] = jac_data[5+6+j][i]; 1364 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 1365 State s = StateFromY(context, Y, x_i); 1366 1367 CeedScalar dY[5], dx0[3] = {0}; 1368 for (int j=0; j<5; j++) dY[j] = dq[j][i]; 1369 State ds = StateFromY_fwd(context, s, dY, x_i, dx0); 1370 1371 State grad_ds[3]; 1372 for (int j=0; j<3; j++) { 1373 CeedScalar dYj[5]; 1374 for (int k=0; k<5; k++) 1375 dYj[k] = Grad_dq[0][k][i] * dXdx[0][j] + 1376 Grad_dq[1][k][i] * dXdx[1][j] + 1377 Grad_dq[2][k][i] * dXdx[2][j]; 1378 grad_ds[j] = StateFromY_fwd(context, s, dYj, x_i, dx0); 1379 } 1380 1381 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 1382 KMStrainRate(grad_ds, dstrain_rate); 1383 NewtonianStress(context, dstrain_rate, dkmstress); 1384 KMUnpack(dkmstress, dstress); 1385 KMUnpack(kmstress, stress); 1386 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 1387 1388 StateConservative dF_inviscid[3]; 1389 FluxInviscid_fwd(context, s, ds, dF_inviscid); 1390 1391 // Total flux 1392 CeedScalar dFlux[5][3]; 1393 for (int j=0; j<3; j++) { 1394 dFlux[0][j] = dF_inviscid[j].density; 1395 for (int k=0; k<3; k++) 1396 dFlux[k+1][j] = dF_inviscid[j].momentum[k] - dstress[k][j]; 1397 dFlux[4][j] = dF_inviscid[j].E_total + dFe[j]; 1398 } 1399 1400 for (int j=0; j<3; j++) { 1401 for (int k=0; k<5; k++) { 1402 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * dFlux[k][0] + 1403 dXdx[j][1] * dFlux[k][1] + 1404 dXdx[j][2] * dFlux[k][2]); 1405 } 1406 } 1407 1408 const CeedScalar dbody_force[5] = {0, 1409 ds.U.density *g[0], 1410 ds.U.density *g[1], 1411 ds.U.density *g[2], 1412 0 1413 }; 1414 CeedScalar dU[5] = {0.}; 1415 dU[0] = ds.U.density; 1416 for (CeedInt j=0; j<3; j++) 1417 dU[j+1] = ds.U.momentum[j]; 1418 dU[4] = ds.U.E_total; 1419 1420 for (int j=0; j<5; j++) 1421 v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]); 1422 1423 if (1) { 1424 CeedScalar jacob_F_conv[3][5][5] = {0}; 1425 computeFluxJacobian_NS(jacob_F_conv, s.U.density, s.Y.velocity, s.U.E_total, 1426 gamma, g, x_i); 1427 CeedScalar grad_dU[5][3]; 1428 for (int j=0; j<3; j++) { 1429 grad_dU[0][j] = grad_ds[j].U.density; 1430 for (int k=0; k<3; k++) grad_dU[k+1][j] = grad_ds[j].U.momentum[k]; 1431 grad_dU[4][j] = grad_ds[j].U.E_total; 1432 } 1433 CeedScalar dstrong_conv[5] = {0.}; 1434 for (int j=0; j<3; j++) 1435 for (int k=0; k<5; k++) 1436 for (int l=0; l<5; l++) 1437 dstrong_conv[k] += jacob_F_conv[j][k][l] * grad_dU[l][j]; 1438 1439 CeedScalar dstrong_res[5]; 1440 for (int j=0; j<5; j++) 1441 dstrong_res[j] = context->ijacobian_time_shift * dU[j] + 1442 dstrong_conv[j] - 1443 dbody_force[j]; 1444 1445 CeedScalar dtau_strong_res[5] = {0.}, 1446 dtau_strong_res_conservative[5] = {0.}; 1447 dtau_strong_res[0] = Tau_d[0] * dstrong_res[0]; 1448 dtau_strong_res[1] = Tau_d[1] * dstrong_res[1]; 1449 dtau_strong_res[2] = Tau_d[1] * dstrong_res[2]; 1450 dtau_strong_res[3] = Tau_d[1] * dstrong_res[3]; 1451 dtau_strong_res[4] = Tau_d[2] * dstrong_res[4]; 1452 PrimitiveToConservative_fwd(s.U.density, s.Y.velocity, s.U.E_total, Rd, cv, 1453 dtau_strong_res, dtau_strong_res_conservative); 1454 CeedScalar dstab[5][3] = {0}; 1455 for (int j=0; j<3; j++) 1456 for (int k=0; k<5; k++) 1457 for (int l=0; l<5; l++) 1458 dstab[k][j] += jacob_F_conv[j][k][l] * dtau_strong_res_conservative[l]; 1459 1460 for (int j=0; j<5; j++) 1461 for (int k=0; k<3; k++) 1462 Grad_v[k][j][i] += wdetJ*(dstab[j][0] * dXdx[k][0] + 1463 dstab[j][1] * dXdx[k][1] + 1464 dstab[j][2] * dXdx[k][2]); 1465 1466 } 1467 } // End Quadrature Point Loop 1468 return 0; 1469 } 1470 // ***************************************************************************** 1471 1472 #endif // newtonian_h 1473