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 U[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 797 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 798 const State s = StateFromU(context, U, x_i); 799 800 const CeedScalar wdetJb = (is_implicit ? -1. : 1.) * q_data_sur[0][i]; 801 // ---- Normal vect 802 const CeedScalar norm[3] = {q_data_sur[1][i], 803 q_data_sur[2][i], 804 q_data_sur[3][i] 805 }; 806 807 const CeedScalar dXdx[2][3] = { 808 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 809 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 810 }; 811 812 State grad_s[3]; 813 for (CeedInt j=0; j<3; j++) { 814 CeedScalar dx_i[3] = {0}, dU[5]; 815 for (CeedInt k=0; k<5; k++) 816 dU[k] = Grad_q[0][k][i] * dXdx[0][j] + 817 Grad_q[1][k][i] * dXdx[1][j]; 818 dx_i[j] = 1.; 819 grad_s[j] = StateFromU_fwd(context, s, dU, x_i, dx_i); 820 } 821 822 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 823 KMStrainRate(grad_s, strain_rate); 824 NewtonianStress(context, strain_rate, kmstress); 825 KMUnpack(kmstress, stress); 826 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 827 828 StateConservative F_inviscid[3]; 829 FluxInviscid(context, s, F_inviscid); 830 831 CeedScalar Flux[5] = {0.}; 832 for (int j=0; j<3; j++) { 833 Flux[0] += F_inviscid[j].density * norm[j]; 834 for (int k=0; k<3; k++) 835 Flux[k+1] += (F_inviscid[j].momentum[k] - stress[k][j]) * norm[j]; 836 Flux[4] += (F_inviscid[j].E_total + Fe[j])*norm[j]; 837 } 838 839 // -- Density 840 v[0][i] = -wdetJb * Flux[0]; 841 842 // -- Momentum 843 for (CeedInt j=0; j<3; j++) 844 v[j+1][i] = -wdetJb * Flux[j+1]; 845 846 // -- Total Energy Density 847 v[4][i] = -wdetJb * Flux[4]; 848 849 jac_data_sur[0][i] = s.U.density; 850 jac_data_sur[1][i] = s.Y.velocity[0]; 851 jac_data_sur[2][i] = s.Y.velocity[1]; 852 jac_data_sur[3][i] = s.Y.velocity[2]; 853 jac_data_sur[4][i] = s.U.E_total; 854 for (int j=0; j<6; j++) jac_data_sur[5+j][i] = kmstress[j]; 855 } 856 return 0; 857 } 858 859 // Jacobian for "set nothing" boundary integral 860 CEED_QFUNCTION(BoundaryIntegral_Jacobian)(void *ctx, CeedInt Q, 861 const CeedScalar *const *in, 862 CeedScalar *const *out) { 863 // *INDENT-OFF* 864 // Inputs 865 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 866 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 867 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 868 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 869 (*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 870 // Outputs 871 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 872 // *INDENT-ON* 873 874 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 875 const bool implicit = context->is_implicit; 876 877 CeedPragmaSIMD 878 // Quadrature Point Loop 879 for (CeedInt i=0; i<Q; i++) { 880 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 881 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 882 const CeedScalar norm[3] = {q_data_sur[1][i], 883 q_data_sur[2][i], 884 q_data_sur[3][i] 885 }; 886 const CeedScalar dXdx[2][3] = { 887 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 888 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 889 }; 890 891 CeedScalar U[5], kmstress[6], dU[5], dx_i[3] = {0.}; 892 for (int j=0; j<5; j++) U[j] = jac_data_sur[j][i]; 893 for (int j=0; j<6; j++) kmstress[j] = jac_data_sur[5+j][i]; 894 for (int j=0; j<3; j++) U[j+1] *= U[0]; 895 for (int j=0; j<5; j++) dU[j] = dq[j][i]; 896 State s = StateFromU(context, U, x_i); 897 State ds = StateFromU_fwd(context, s, dU, x_i, dx_i); 898 899 State grad_ds[3]; 900 for (CeedInt j=0; j<3; j++) { 901 CeedScalar dx_i[3] = {0}, dUj[5]; 902 for (CeedInt k=0; k<5; k++) 903 dUj[k] = Grad_dq[0][k][i] * dXdx[0][j] + 904 Grad_dq[1][k][i] * dXdx[1][j]; 905 dx_i[j] = 1.; 906 grad_ds[j] = StateFromU_fwd(context, s, dUj, x_i, dx_i); 907 } 908 909 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 910 KMStrainRate(grad_ds, dstrain_rate); 911 NewtonianStress(context, dstrain_rate, dkmstress); 912 KMUnpack(dkmstress, dstress); 913 KMUnpack(kmstress, stress); 914 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 915 916 StateConservative dF_inviscid[3]; 917 FluxInviscid_fwd(context, s, ds, dF_inviscid); 918 919 CeedScalar dFlux[5] = {0.}; 920 for (int j=0; j<3; j++) { 921 dFlux[0] += dF_inviscid[j].density * norm[j]; 922 for (int k=0; k<3; k++) 923 dFlux[k+1] += (dF_inviscid[j].momentum[k] - dstress[k][j]) * norm[j]; 924 dFlux[4] += (dF_inviscid[j].E_total + dFe[j]) * norm[j]; 925 } 926 927 for (int j=0; j<5; j++) 928 v[j][i] = -wdetJb * dFlux[j]; 929 } // End Quadrature Point Loop 930 return 0; 931 } 932 933 // Outflow boundary condition, weakly setting a constant pressure 934 CEED_QFUNCTION(PressureOutflow)(void *ctx, CeedInt Q, 935 const CeedScalar *const *in, 936 CeedScalar *const *out) { 937 // *INDENT-OFF* 938 // Inputs 939 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 940 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 941 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 942 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 943 // Outputs 944 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 945 (*jac_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[1]; 946 // *INDENT-ON* 947 948 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 949 const bool implicit = context->is_implicit; 950 const CeedScalar P0 = context->P0; 951 952 CeedPragmaSIMD 953 // Quadrature Point Loop 954 for (CeedInt i=0; i<Q; i++) { 955 // Setup 956 // -- Interp in 957 const CeedScalar U[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 958 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 959 State s = StateFromU(context, U, x_i); 960 s.Y.pressure = P0; 961 962 // -- Interp-to-Interp q_data 963 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 964 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 965 // We can effect this by swapping the sign on this weight 966 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 967 968 // ---- Normal vect 969 const CeedScalar norm[3] = {q_data_sur[1][i], 970 q_data_sur[2][i], 971 q_data_sur[3][i] 972 }; 973 974 const CeedScalar dXdx[2][3] = { 975 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 976 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 977 }; 978 979 State grad_s[3]; 980 for (CeedInt j=0; j<3; j++) { 981 CeedScalar dx_i[3] = {0}, dU[5]; 982 for (CeedInt k=0; k<5; k++) 983 dU[k] = Grad_q[0][k][i] * dXdx[0][j] + 984 Grad_q[1][k][i] * dXdx[1][j]; 985 dx_i[j] = 1.; 986 grad_s[j] = StateFromU_fwd(context, s, dU, x_i, dx_i); 987 } 988 989 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 990 KMStrainRate(grad_s, strain_rate); 991 NewtonianStress(context, strain_rate, kmstress); 992 KMUnpack(kmstress, stress); 993 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 994 995 StateConservative F_inviscid[3]; 996 FluxInviscid(context, s, F_inviscid); 997 998 CeedScalar Flux[5] = {0.}; 999 for (int j=0; j<3; j++) { 1000 Flux[0] += F_inviscid[j].density * norm[j]; 1001 for (int k=0; k<3; k++) 1002 Flux[k+1] += (F_inviscid[j].momentum[k] - stress[k][j]) * norm[j]; 1003 Flux[4] += (F_inviscid[j].E_total + Fe[j])*norm[j]; 1004 } 1005 1006 // -- Density 1007 v[0][i] = -wdetJb * Flux[0]; 1008 1009 // -- Momentum 1010 for (CeedInt j=0; j<3; j++) 1011 v[j+1][i] = -wdetJb * Flux[j+1]; 1012 1013 // -- Total Energy Density 1014 v[4][i] = -wdetJb * Flux[4]; 1015 1016 // Save values for Jacobian 1017 jac_data_sur[0][i] = s.U.density; 1018 jac_data_sur[1][i] = s.Y.velocity[0]; 1019 jac_data_sur[2][i] = s.Y.velocity[1]; 1020 jac_data_sur[3][i] = s.Y.velocity[2]; 1021 jac_data_sur[4][i] = s.U.E_total; 1022 for (int j=0; j<6; j++) jac_data_sur[5+j][i] = kmstress[j]; 1023 } // End Quadrature Point Loop 1024 return 0; 1025 } 1026 1027 // Jacobian for weak-pressure outflow boundary condition 1028 CEED_QFUNCTION(PressureOutflow_Jacobian)(void *ctx, CeedInt Q, 1029 const CeedScalar *const *in, 1030 CeedScalar *const *out) { 1031 // *INDENT-OFF* 1032 // Inputs 1033 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 1034 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 1035 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 1036 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 1037 (*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 1038 // Outputs 1039 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 1040 // *INDENT-ON* 1041 1042 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 1043 const bool implicit = context->is_implicit; 1044 1045 CeedPragmaSIMD 1046 // Quadrature Point Loop 1047 for (CeedInt i=0; i<Q; i++) { 1048 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 1049 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 1050 const CeedScalar norm[3] = {q_data_sur[1][i], 1051 q_data_sur[2][i], 1052 q_data_sur[3][i] 1053 }; 1054 const CeedScalar dXdx[2][3] = { 1055 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 1056 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 1057 }; 1058 1059 CeedScalar U[5], kmstress[6], dU[5], dx_i[3] = {0.}; 1060 for (int j=0; j<5; j++) U[j] = jac_data_sur[j][i]; 1061 for (int j=0; j<6; j++) kmstress[j] = jac_data_sur[5+j][i]; 1062 for (int j=0; j<3; j++) U[j+1] *= U[0]; 1063 for (int j=0; j<5; j++) dU[j] = dq[j][i]; 1064 State s = StateFromU(context, U, x_i); 1065 State ds = StateFromU_fwd(context, s, dU, x_i, dx_i); 1066 s.Y.pressure = context->P0; 1067 ds.Y.pressure = 0.; 1068 1069 State grad_ds[3]; 1070 for (CeedInt j=0; j<3; j++) { 1071 CeedScalar dx_i[3] = {0}, dUj[5]; 1072 for (CeedInt k=0; k<5; k++) 1073 dUj[k] = Grad_dq[0][k][i] * dXdx[0][j] + 1074 Grad_dq[1][k][i] * dXdx[1][j]; 1075 dx_i[j] = 1.; 1076 grad_ds[j] = StateFromU_fwd(context, s, dUj, x_i, dx_i); 1077 } 1078 1079 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 1080 KMStrainRate(grad_ds, dstrain_rate); 1081 NewtonianStress(context, dstrain_rate, dkmstress); 1082 KMUnpack(dkmstress, dstress); 1083 KMUnpack(kmstress, stress); 1084 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 1085 1086 StateConservative dF_inviscid[3]; 1087 FluxInviscid_fwd(context, s, ds, dF_inviscid); 1088 1089 CeedScalar dFlux[5] = {0.}; 1090 for (int j=0; j<3; j++) { 1091 dFlux[0] += dF_inviscid[j].density * norm[j]; 1092 for (int k=0; k<3; k++) 1093 dFlux[k+1] += (dF_inviscid[j].momentum[k] - dstress[k][j]) * norm[j]; 1094 dFlux[4] += (dF_inviscid[j].E_total + dFe[j]) * norm[j]; 1095 } 1096 1097 for (int j=0; j<5; j++) 1098 v[j][i] = -wdetJb * dFlux[j]; 1099 } // End Quadrature Point Loop 1100 return 0; 1101 } 1102 1103 // ***************************************************************************** 1104 // This QFunction implements the Navier-Stokes equations (mentioned above) in 1105 // primitive variables and with implicit time stepping method 1106 // 1107 // ***************************************************************************** 1108 CEED_QFUNCTION(IFunction_Newtonian_Prim)(void *ctx, CeedInt Q, 1109 const CeedScalar *const *in, CeedScalar *const *out) { 1110 // *INDENT-OFF* 1111 // Inputs 1112 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 1113 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 1114 (*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 1115 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 1116 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 1117 // Outputs 1118 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 1119 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1], 1120 (*jac_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[2]; 1121 // *INDENT-ON* 1122 // Context 1123 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 1124 const CeedScalar mu = context->mu; 1125 const CeedScalar cv = context->cv; 1126 const CeedScalar cp = context->cp; 1127 const CeedScalar *g = context->g; 1128 const CeedScalar dt = context->dt; 1129 const CeedScalar gamma = cp / cv; 1130 const CeedScalar Rd = cp - cv; 1131 1132 CeedPragmaSIMD 1133 // Quadrature Point Loop 1134 for (CeedInt i=0; i<Q; i++) { 1135 CeedScalar Y[5]; 1136 for (CeedInt j=0; j<5; j++) Y[j] = q[j][i]; 1137 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 1138 State s = StateFromY(context, Y, x_i); 1139 1140 // -- Interp-to-Interp q_data 1141 const CeedScalar wdetJ = q_data[0][i]; 1142 // -- Interp-to-Grad q_data 1143 // ---- Inverse of change of coordinate matrix: X_i,j 1144 // *INDENT-OFF* 1145 const CeedScalar dXdx[3][3] = {{q_data[1][i], 1146 q_data[2][i], 1147 q_data[3][i]}, 1148 {q_data[4][i], 1149 q_data[5][i], 1150 q_data[6][i]}, 1151 {q_data[7][i], 1152 q_data[8][i], 1153 q_data[9][i]} 1154 }; 1155 // *INDENT-ON* 1156 State grad_s[3]; 1157 for (CeedInt j=0; j<3; j++) { 1158 CeedScalar dx_i[3] = {0}, dY[5]; 1159 for (CeedInt k=0; k<5; k++) 1160 dY[k] = Grad_q[0][k][i] * dXdx[0][j] + 1161 Grad_q[1][k][i] * dXdx[1][j] + 1162 Grad_q[2][k][i] * dXdx[2][j]; 1163 dx_i[j] = 1.; 1164 grad_s[j] = StateFromY_fwd(context, s, dY, x_i, dx_i); 1165 } 1166 1167 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 1168 KMStrainRate(grad_s, strain_rate); 1169 NewtonianStress(context, strain_rate, kmstress); 1170 KMUnpack(kmstress, stress); 1171 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 1172 1173 StateConservative F_inviscid[3]; 1174 FluxInviscid(context, s, F_inviscid); 1175 1176 // Total flux 1177 CeedScalar Flux[5][3]; 1178 for (CeedInt j=0; j<3; j++) { 1179 Flux[0][j] = F_inviscid[j].density; 1180 for (CeedInt k=0; k<3; k++) 1181 Flux[k+1][j] = F_inviscid[j].momentum[k] - stress[k][j]; 1182 Flux[4][j] = F_inviscid[j].E_total + Fe[j]; 1183 } 1184 1185 for (CeedInt j=0; j<3; j++) { 1186 for (CeedInt k=0; k<5; k++) { 1187 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * Flux[k][0] + 1188 dXdx[j][1] * Flux[k][1] + 1189 dXdx[j][2] * Flux[k][2]); 1190 } 1191 } 1192 1193 const CeedScalar body_force[5] = {0, s.U.density *g[0], s.U.density *g[1], s.U.density *g[2], 0}; 1194 1195 CeedScalar Y_dot[5], dx0[3] = {0}; 1196 for (int j=0; j<5; j++) Y_dot[j] = q_dot[j][i]; 1197 State s_dot = StateFromY_fwd(context, s, Y_dot, x_i, dx0); 1198 1199 CeedScalar U_dot[5] = {0.}; 1200 U_dot[0] = s_dot.U.density; 1201 for (CeedInt j=0; j<3; j++) 1202 U_dot[j+1] = s_dot.U.momentum[j]; 1203 U_dot[4] = s_dot.U.E_total; 1204 1205 for (CeedInt j=0; j<5; j++) 1206 v[j][i] = wdetJ * (U_dot[j] - body_force[j]); 1207 1208 // jacob_F_conv[3][5][5] = dF(convective)/dq at each direction 1209 CeedScalar jacob_F_conv[3][5][5] = {0}; 1210 computeFluxJacobian_NS(jacob_F_conv, s.U.density, s.Y.velocity, s.U.E_total, 1211 gamma, g, x_i); 1212 CeedScalar grad_U[5][3]; 1213 for (CeedInt j=0; j<3; j++) { 1214 grad_U[0][j] = grad_s[j].U.density; 1215 for (CeedInt k=0; k<3; k++) grad_U[k+1][j] = grad_s[j].U.momentum[k]; 1216 grad_U[4][j] = grad_s[j].U.E_total; 1217 } 1218 1219 // strong_conv = dF/dq * dq/dx (Strong convection) 1220 CeedScalar strong_conv[5] = {0}; 1221 for (CeedInt j=0; j<3; j++) 1222 for (CeedInt k=0; k<5; k++) 1223 for (CeedInt l=0; l<5; l++) 1224 strong_conv[k] += jacob_F_conv[j][k][l] * grad_U[l][j]; 1225 1226 // Strong residual 1227 CeedScalar strong_res[5]; 1228 for (CeedInt j=0; j<5; j++) 1229 strong_res[j] = U_dot[j] + strong_conv[j] - body_force[j]; 1230 1231 // -- Stabilization method: none, SU, or SUPG 1232 CeedScalar stab[5][3] = {{0.}}; 1233 CeedScalar tau_strong_res[5] = {0.}, tau_strong_res_conservative[5] = {0}; 1234 CeedScalar tau_strong_conv[5] = {0.}, tau_strong_conv_conservative[5] = {0}; 1235 CeedScalar Tau_d[3] = {0.}; 1236 switch (context->stabilization) { 1237 case STAB_NONE: // Galerkin 1238 break; 1239 case STAB_SU: // SU 1240 Tau_diagPrim(Tau_d, dXdx, s.Y.velocity, cv, context, mu, dt, s.U.density); 1241 tau_strong_conv[0] = Tau_d[0] * strong_conv[0]; 1242 tau_strong_conv[1] = Tau_d[1] * strong_conv[1]; 1243 tau_strong_conv[2] = Tau_d[1] * strong_conv[2]; 1244 tau_strong_conv[3] = Tau_d[1] * strong_conv[3]; 1245 tau_strong_conv[4] = Tau_d[2] * strong_conv[4]; 1246 PrimitiveToConservative_fwd(s.U.density, s.Y.velocity, s.U.E_total, Rd, cv, 1247 tau_strong_conv, tau_strong_conv_conservative); 1248 for (CeedInt j=0; j<3; j++) 1249 for (CeedInt k=0; k<5; k++) 1250 for (CeedInt l=0; l<5; l++) 1251 stab[k][j] += jacob_F_conv[j][k][l] * tau_strong_conv_conservative[l]; 1252 1253 for (CeedInt j=0; j<5; j++) 1254 for (CeedInt k=0; k<3; k++) 1255 Grad_v[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 1256 stab[j][1] * dXdx[k][1] + 1257 stab[j][2] * dXdx[k][2]); 1258 1259 break; 1260 case STAB_SUPG: // SUPG 1261 Tau_diagPrim(Tau_d, dXdx, s.Y.velocity, cv, context, mu, dt, s.U.density); 1262 tau_strong_res[0] = Tau_d[0] * strong_res[0]; 1263 tau_strong_res[1] = Tau_d[1] * strong_res[1]; 1264 tau_strong_res[2] = Tau_d[1] * strong_res[2]; 1265 tau_strong_res[3] = Tau_d[1] * strong_res[3]; 1266 tau_strong_res[4] = Tau_d[2] * strong_res[4]; 1267 1268 PrimitiveToConservative_fwd(s.U.density, s.Y.velocity, s.U.E_total, Rd, cv, 1269 tau_strong_res, tau_strong_res_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_res_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 break; 1281 } 1282 for (CeedInt j=0; j<5; j++) jac_data[j][i] = Y[j]; 1283 for (CeedInt j=0; j<6; j++) jac_data[5+j][i] = kmstress[j]; 1284 for (CeedInt j=0; j<3; j++) jac_data[5+6+j][i] = Tau_d[j]; 1285 1286 } // End Quadrature Point Loop 1287 1288 // Return 1289 return 0; 1290 } 1291 1292 // ***************************************************************************** 1293 // This QFunction implements the jacobean of the Navier-Stokes equations 1294 // in primitive variables for implicit time stepping method. 1295 // 1296 // ***************************************************************************** 1297 CEED_QFUNCTION(IJacobian_Newtonian_Prim)(void *ctx, CeedInt Q, 1298 const CeedScalar *const *in, CeedScalar *const *out) { 1299 // *INDENT-OFF* 1300 // Inputs 1301 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 1302 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 1303 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 1304 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 1305 (*jac_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 1306 // Outputs 1307 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 1308 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 1309 // *INDENT-ON* 1310 // Context 1311 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 1312 const CeedScalar *g = context->g; 1313 const CeedScalar cp = context->cp; 1314 const CeedScalar cv = context->cv; 1315 const CeedScalar Rd = cp - cv; 1316 const CeedScalar gamma = cp / cv; 1317 1318 CeedPragmaSIMD 1319 // Quadrature Point Loop 1320 for (CeedInt i=0; i<Q; i++) { 1321 // -- Interp-to-Interp q_data 1322 const CeedScalar wdetJ = q_data[0][i]; 1323 // -- Interp-to-Grad q_data 1324 // ---- Inverse of change of coordinate matrix: X_i,j 1325 // *INDENT-OFF* 1326 const CeedScalar dXdx[3][3] = {{q_data[1][i], 1327 q_data[2][i], 1328 q_data[3][i]}, 1329 {q_data[4][i], 1330 q_data[5][i], 1331 q_data[6][i]}, 1332 {q_data[7][i], 1333 q_data[8][i], 1334 q_data[9][i]} 1335 }; 1336 // *INDENT-ON* 1337 1338 CeedScalar Y[5], kmstress[6], Tau_d[3] __attribute((unused)); 1339 for (int j=0; j<5; j++) Y[j] = jac_data[j][i]; 1340 for (int j=0; j<6; j++) kmstress[j] = jac_data[5+j][i]; 1341 for (int j=0; j<3; j++) Tau_d[j] = jac_data[5+6+j][i]; 1342 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 1343 State s = StateFromY(context, Y, x_i); 1344 1345 CeedScalar dY[5], dx0[3] = {0}; 1346 for (int j=0; j<5; j++) dY[j] = dq[j][i]; 1347 State ds = StateFromY_fwd(context, s, dY, x_i, dx0); 1348 1349 State grad_ds[3]; 1350 for (int j=0; j<3; j++) { 1351 CeedScalar dYj[5]; 1352 for (int k=0; k<5; k++) 1353 dYj[k] = Grad_dq[0][k][i] * dXdx[0][j] + 1354 Grad_dq[1][k][i] * dXdx[1][j] + 1355 Grad_dq[2][k][i] * dXdx[2][j]; 1356 grad_ds[j] = StateFromY_fwd(context, s, dYj, x_i, dx0); 1357 } 1358 1359 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 1360 KMStrainRate(grad_ds, dstrain_rate); 1361 NewtonianStress(context, dstrain_rate, dkmstress); 1362 KMUnpack(dkmstress, dstress); 1363 KMUnpack(kmstress, stress); 1364 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 1365 1366 StateConservative dF_inviscid[3]; 1367 FluxInviscid_fwd(context, s, ds, dF_inviscid); 1368 1369 // Total flux 1370 CeedScalar dFlux[5][3]; 1371 for (int j=0; j<3; j++) { 1372 dFlux[0][j] = dF_inviscid[j].density; 1373 for (int k=0; k<3; k++) 1374 dFlux[k+1][j] = dF_inviscid[j].momentum[k] - dstress[k][j]; 1375 dFlux[4][j] = dF_inviscid[j].E_total + dFe[j]; 1376 } 1377 1378 for (int j=0; j<3; j++) { 1379 for (int k=0; k<5; k++) { 1380 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * dFlux[k][0] + 1381 dXdx[j][1] * dFlux[k][1] + 1382 dXdx[j][2] * dFlux[k][2]); 1383 } 1384 } 1385 1386 const CeedScalar dbody_force[5] = {0, 1387 ds.U.density *g[0], 1388 ds.U.density *g[1], 1389 ds.U.density *g[2], 1390 0 1391 }; 1392 CeedScalar dU[5] = {0.}; 1393 dU[0] = ds.U.density; 1394 for (CeedInt j=0; j<3; j++) 1395 dU[j+1] = ds.U.momentum[j]; 1396 dU[4] = ds.U.E_total; 1397 1398 for (int j=0; j<5; j++) 1399 v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]); 1400 1401 if (1) { 1402 CeedScalar jacob_F_conv[3][5][5] = {0}; 1403 computeFluxJacobian_NS(jacob_F_conv, s.U.density, s.Y.velocity, s.U.E_total, 1404 gamma, g, x_i); 1405 CeedScalar grad_dU[5][3]; 1406 for (int j=0; j<3; j++) { 1407 grad_dU[0][j] = grad_ds[j].U.density; 1408 for (int k=0; k<3; k++) grad_dU[k+1][j] = grad_ds[j].U.momentum[k]; 1409 grad_dU[4][j] = grad_ds[j].U.E_total; 1410 } 1411 CeedScalar dstrong_conv[5] = {0.}; 1412 for (int j=0; j<3; j++) 1413 for (int k=0; k<5; k++) 1414 for (int l=0; l<5; l++) 1415 dstrong_conv[k] += jacob_F_conv[j][k][l] * grad_dU[l][j]; 1416 1417 CeedScalar dstrong_res[5]; 1418 for (int j=0; j<5; j++) 1419 dstrong_res[j] = context->ijacobian_time_shift * dU[j] + 1420 dstrong_conv[j] - 1421 dbody_force[j]; 1422 1423 CeedScalar dtau_strong_res[5] = {0.}, 1424 dtau_strong_res_conservative[5] = {0.}; 1425 dtau_strong_res[0] = Tau_d[0] * dstrong_res[0]; 1426 dtau_strong_res[1] = Tau_d[1] * dstrong_res[1]; 1427 dtau_strong_res[2] = Tau_d[1] * dstrong_res[2]; 1428 dtau_strong_res[3] = Tau_d[1] * dstrong_res[3]; 1429 dtau_strong_res[4] = Tau_d[2] * dstrong_res[4]; 1430 PrimitiveToConservative_fwd(s.U.density, s.Y.velocity, s.U.E_total, Rd, cv, 1431 dtau_strong_res, dtau_strong_res_conservative); 1432 CeedScalar dstab[5][3] = {0}; 1433 for (int j=0; j<3; j++) 1434 for (int k=0; k<5; k++) 1435 for (int l=0; l<5; l++) 1436 dstab[k][j] += jacob_F_conv[j][k][l] * dtau_strong_res_conservative[l]; 1437 1438 for (int j=0; j<5; j++) 1439 for (int k=0; k<3; k++) 1440 Grad_v[k][j][i] += wdetJ*(dstab[j][0] * dXdx[k][0] + 1441 dstab[j][1] * dXdx[k][1] + 1442 dstab[j][2] * dXdx[k][2]); 1443 1444 } 1445 } // End Quadrature Point Loop 1446 return 0; 1447 } 1448 // ***************************************************************************** 1449 1450 #endif // newtonian_h 1451