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 <ceed.h> 16 #include <math.h> 17 #include "newtonian_state.h" 18 #include "newtonian_types.h" 19 #include "stabilization.h" 20 #include "utils.h" 21 22 // ***************************************************************************** 23 // This QFunction sets a "still" initial condition for generic Newtonian IG problems 24 // ***************************************************************************** 25 CEED_QFUNCTION(ICsNewtonianIG)(void *ctx, CeedInt Q, 26 const CeedScalar *const *in, CeedScalar *const *out) { 27 // Inputs 28 const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 29 30 // Outputs 31 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 32 33 // Context 34 const SetupContext context = (SetupContext)ctx; 35 const CeedScalar theta0 = context->theta0; 36 const CeedScalar P0 = context->P0; 37 const CeedScalar cv = context->cv; 38 const CeedScalar cp = context->cp; 39 const CeedScalar *g = context->g; 40 const CeedScalar Rd = cp - cv; 41 42 // Quadrature Point Loop 43 CeedPragmaSIMD 44 for (CeedInt i=0; i<Q; i++) { 45 CeedScalar q[5] = {0.}; 46 47 // Setup 48 // -- Coordinates 49 const CeedScalar x[3] = {X[0][i], X[1][i], X[2][i]}; 50 const CeedScalar e_potential = -Dot3(g, x); 51 52 // -- Density 53 const CeedScalar rho = P0 / (Rd*theta0); 54 55 // Initial Conditions 56 q[0] = rho; 57 q[1] = 0.0; 58 q[2] = 0.0; 59 q[3] = 0.0; 60 q[4] = rho * (cv*theta0 + e_potential); 61 62 for (CeedInt j=0; j<5; j++) 63 q0[j][i] = q[j]; 64 65 } // End of Quadrature Point Loop 66 return 0; 67 } 68 69 // ***************************************************************************** 70 // This QFunction sets a "still" initial condition for generic Newtonian IG 71 // problems in primitive variables 72 // ***************************************************************************** 73 CEED_QFUNCTION(ICsNewtonianIG_Prim)(void *ctx, CeedInt Q, 74 const CeedScalar *const *in, CeedScalar *const *out) { 75 // Outputs 76 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 77 78 // Context 79 const SetupContext context = (SetupContext)ctx; 80 const CeedScalar theta0 = context->theta0; 81 const CeedScalar P0 = context->P0; 82 83 // Quadrature Point Loop 84 CeedPragmaSIMD 85 for (CeedInt i=0; i<Q; i++) { 86 CeedScalar q[5] = {0.}; 87 88 // Initial Conditions 89 q[0] = P0; 90 q[1] = 0.0; 91 q[2] = 0.0; 92 q[3] = 0.0; 93 q[4] = theta0; 94 95 for (CeedInt j=0; j<5; j++) 96 q0[j][i] = q[j]; 97 98 } // End of Quadrature Point Loop 99 return 0; 100 } 101 102 // ***************************************************************************** 103 // This QFunction implements the following formulation of Navier-Stokes with 104 // explicit time stepping method 105 // 106 // This is 3D compressible Navier-Stokes in conservation form with state 107 // variables of density, momentum density, and total energy density. 108 // 109 // State Variables: q = ( rho, U1, U2, U3, E ) 110 // rho - Mass Density 111 // Ui - Momentum Density, Ui = rho ui 112 // E - Total Energy Density, E = rho (cv T + (u u)/2 + g z) 113 // 114 // Navier-Stokes Equations: 115 // drho/dt + div( U ) = 0 116 // dU/dt + div( rho (u x u) + P I3 ) + rho g khat = div( Fu ) 117 // dE/dt + div( (E + P) u ) = div( Fe ) 118 // 119 // Viscous Stress: 120 // Fu = mu (grad( u ) + grad( u )^T + lambda div ( u ) I3) 121 // 122 // Thermal Stress: 123 // Fe = u Fu + k grad( T ) 124 // Equation of State 125 // P = (gamma - 1) (E - rho (u u) / 2 - rho g z) 126 // 127 // Stabilization: 128 // Tau = diag(TauC, TauM, TauM, TauM, TauE) 129 // f1 = rho sqrt(ui uj gij) 130 // gij = dXi/dX * dXi/dX 131 // TauC = Cc f1 / (8 gii) 132 // TauM = min( 1 , 1 / f1 ) 133 // TauE = TauM / (Ce cv) 134 // 135 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 136 // 137 // Constants: 138 // lambda = - 2 / 3, From Stokes hypothesis 139 // mu , Dynamic viscosity 140 // k , Thermal conductivity 141 // cv , Specific heat, constant volume 142 // cp , Specific heat, constant pressure 143 // g , Gravity 144 // gamma = cp / cv, Specific heat ratio 145 // 146 // We require the product of the inverse of the Jacobian (dXdx_j,k) and 147 // its transpose (dXdx_k,j) to properly compute integrals of the form: 148 // int( gradv gradu ) 149 // 150 // ***************************************************************************** 151 CEED_QFUNCTION(RHSFunction_Newtonian)(void *ctx, CeedInt Q, 152 const CeedScalar *const *in, CeedScalar *const *out) { 153 // *INDENT-OFF* 154 // Inputs 155 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 156 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 157 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 158 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 159 // Outputs 160 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 161 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 162 // *INDENT-ON* 163 164 // Context 165 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 166 const CeedScalar *g = context->g; 167 const CeedScalar dt = context->dt; 168 169 CeedPragmaSIMD 170 // Quadrature Point Loop 171 for (CeedInt i=0; i<Q; i++) { 172 CeedScalar U[5]; 173 for (int j=0; j<5; j++) U[j] = q[j][i]; 174 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 175 State s = StateFromU(context, U, x_i); 176 177 // -- Interp-to-Interp q_data 178 const CeedScalar wdetJ = q_data[0][i]; 179 // -- Interp-to-Grad q_data 180 // ---- Inverse of change of coordinate matrix: X_i,j 181 // *INDENT-OFF* 182 const CeedScalar dXdx[3][3] = {{q_data[1][i], q_data[2][i], q_data[3][i]}, 183 {q_data[4][i], q_data[5][i], q_data[6][i]}, 184 {q_data[7][i], q_data[8][i], q_data[9][i]} 185 }; 186 // *INDENT-ON* 187 State grad_s[3]; 188 for (CeedInt j=0; j<3; j++) { 189 CeedScalar dx_i[3] = {0}, dU[5]; 190 for (CeedInt k=0; k<5; k++) 191 dU[k] = Grad_q[0][k][i] * dXdx[0][j] + 192 Grad_q[1][k][i] * dXdx[1][j] + 193 Grad_q[2][k][i] * dXdx[2][j]; 194 dx_i[j] = 1.; 195 grad_s[j] = StateFromU_fwd(context, s, dU, x_i, dx_i); 196 } 197 198 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 199 KMStrainRate(grad_s, strain_rate); 200 NewtonianStress(context, strain_rate, kmstress); 201 KMUnpack(kmstress, stress); 202 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 203 204 StateConservative F_inviscid[3]; 205 FluxInviscid(context, s, F_inviscid); 206 207 // Total flux 208 CeedScalar Flux[5][3]; 209 FluxTotal(F_inviscid, stress, Fe, Flux); 210 211 for (CeedInt j=0; j<3; j++) 212 for (CeedInt k=0; k<5; k++) 213 Grad_v[j][k][i] = wdetJ * (dXdx[j][0] * Flux[k][0] + 214 dXdx[j][1] * Flux[k][1] + 215 dXdx[j][2] * Flux[k][2]); 216 217 const CeedScalar body_force[5] = {0, s.U.density *g[0], s.U.density *g[1], s.U.density *g[2], 0}; 218 for (int j=0; j<5; j++) 219 v[j][i] = wdetJ * body_force[j]; 220 221 // -- Stabilization method: none (Galerkin), SU, or SUPG 222 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}; 223 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 224 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, x_i, stab); 225 226 for (CeedInt j=0; j<5; j++) 227 for (CeedInt k=0; k<3; k++) 228 Grad_v[k][j][i] -= wdetJ*(stab[j][0] * dXdx[k][0] + 229 stab[j][1] * dXdx[k][1] + 230 stab[j][2] * dXdx[k][2]); 231 232 } // End Quadrature Point Loop 233 234 // Return 235 return 0; 236 } 237 238 // ***************************************************************************** 239 // This QFunction implements the Navier-Stokes equations (mentioned above) with 240 // implicit time stepping method 241 // 242 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 243 // SUPG = Galerkin + grad(v) . ( Ai^T * Tau * (q_dot + Aj q,j - body force) ) 244 // (diffussive terms will be added later) 245 // 246 // ***************************************************************************** 247 CEED_QFUNCTION(IFunction_Newtonian)(void *ctx, CeedInt Q, 248 const CeedScalar *const *in, CeedScalar *const *out) { 249 // *INDENT-OFF* 250 // Inputs 251 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 252 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 253 (*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 254 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 255 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 256 // Outputs 257 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 258 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1], 259 (*jac_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[2]; 260 // *INDENT-ON* 261 // Context 262 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 263 const CeedScalar *g = context->g; 264 const CeedScalar dt = context->dt; 265 266 CeedPragmaSIMD 267 // Quadrature Point Loop 268 for (CeedInt i=0; i<Q; i++) { 269 CeedScalar U[5]; 270 for (CeedInt j=0; j<5; j++) U[j] = q[j][i]; 271 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 272 State s = StateFromU(context, U, x_i); 273 274 // -- Interp-to-Interp q_data 275 const CeedScalar wdetJ = q_data[0][i]; 276 // -- Interp-to-Grad q_data 277 // ---- Inverse of change of coordinate matrix: X_i,j 278 // *INDENT-OFF* 279 const CeedScalar dXdx[3][3] = {{q_data[1][i], q_data[2][i], q_data[3][i]}, 280 {q_data[4][i], q_data[5][i], q_data[6][i]}, 281 {q_data[7][i], q_data[8][i], q_data[9][i]} 282 }; 283 // *INDENT-ON* 284 State grad_s[3]; 285 for (CeedInt j=0; j<3; j++) { 286 CeedScalar dx_i[3] = {0}, dU[5]; 287 for (CeedInt k=0; k<5; k++) 288 dU[k] = Grad_q[0][k][i] * dXdx[0][j] + 289 Grad_q[1][k][i] * dXdx[1][j] + 290 Grad_q[2][k][i] * dXdx[2][j]; 291 dx_i[j] = 1.; 292 grad_s[j] = StateFromU_fwd(context, s, dU, x_i, dx_i); 293 } 294 295 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 296 KMStrainRate(grad_s, strain_rate); 297 NewtonianStress(context, strain_rate, kmstress); 298 KMUnpack(kmstress, stress); 299 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 300 301 StateConservative F_inviscid[3]; 302 FluxInviscid(context, s, F_inviscid); 303 304 // Total flux 305 CeedScalar Flux[5][3]; 306 FluxTotal(F_inviscid, stress, Fe, Flux); 307 308 for (CeedInt j=0; j<3; j++) 309 for (CeedInt k=0; k<5; k++) 310 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * Flux[k][0] + 311 dXdx[j][1] * Flux[k][1] + 312 dXdx[j][2] * Flux[k][2]); 313 314 const CeedScalar body_force[5] = {0, s.U.density *g[0], s.U.density *g[1], s.U.density *g[2], 0}; 315 for (CeedInt j=0; j<5; j++) 316 v[j][i] = wdetJ * (q_dot[j][i] - body_force[j]); 317 318 // -- Stabilization method: none (Galerkin), SU, or SUPG 319 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}; 320 for (CeedInt j=0; j<5; j++) U_dot[j] = q_dot[j][i]; 321 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 322 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, x_i, stab); 323 324 for (CeedInt j=0; j<5; j++) 325 for (CeedInt k=0; k<3; k++) 326 Grad_v[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 327 stab[j][1] * dXdx[k][1] + 328 stab[j][2] * dXdx[k][2]); 329 330 for (CeedInt j=0; j<5; j++) jac_data[j][i] = U[j]; 331 for (CeedInt j=0; j<6; j++) jac_data[5+j][i] = kmstress[j]; 332 for (CeedInt j=0; j<3; j++) jac_data[5+6+j][i] = Tau_d[j]; 333 334 } // End Quadrature Point Loop 335 336 // Return 337 return 0; 338 } 339 340 // ***************************************************************************** 341 // This QFunction implements the jacobean of the Navier-Stokes equations 342 // for implicit time stepping method. 343 // 344 // ***************************************************************************** 345 CEED_QFUNCTION(IJacobian_Newtonian)(void *ctx, CeedInt Q, 346 const CeedScalar *const *in, 347 CeedScalar *const *out) { 348 // *INDENT-OFF* 349 // Inputs 350 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 351 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 352 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 353 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 354 (*jac_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 355 // Outputs 356 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 357 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 358 // *INDENT-ON* 359 // Context 360 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 361 const CeedScalar *g = context->g; 362 363 CeedPragmaSIMD 364 // Quadrature Point Loop 365 for (CeedInt i=0; i<Q; i++) { 366 // -- Interp-to-Interp q_data 367 const CeedScalar wdetJ = q_data[0][i]; 368 // -- Interp-to-Grad q_data 369 // ---- Inverse of change of coordinate matrix: X_i,j 370 // *INDENT-OFF* 371 const CeedScalar dXdx[3][3] = {{q_data[1][i], q_data[2][i], q_data[3][i]}, 372 {q_data[4][i], q_data[5][i], q_data[6][i]}, 373 {q_data[7][i], q_data[8][i], q_data[9][i]} 374 }; 375 // *INDENT-ON* 376 377 CeedScalar U[5], kmstress[6], Tau_d[3] __attribute((unused)); 378 for (int j=0; j<5; j++) U[j] = jac_data[j][i]; 379 for (int j=0; j<6; j++) kmstress[j] = jac_data[5+j][i]; 380 for (int j=0; j<3; j++) Tau_d[j] = jac_data[5+6+j][i]; 381 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 382 State s = StateFromU(context, U, x_i); 383 384 CeedScalar dU[5], dx0[3] = {0}; 385 for (int j=0; j<5; j++) dU[j] = dq[j][i]; 386 State ds = StateFromU_fwd(context, s, dU, x_i, dx0); 387 388 State grad_ds[3]; 389 for (int j=0; j<3; j++) { 390 CeedScalar dUj[5]; 391 for (int k=0; k<5; k++) 392 dUj[k] = Grad_dq[0][k][i] * dXdx[0][j] + 393 Grad_dq[1][k][i] * dXdx[1][j] + 394 Grad_dq[2][k][i] * dXdx[2][j]; 395 grad_ds[j] = StateFromU_fwd(context, s, dUj, x_i, dx0); 396 } 397 398 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 399 KMStrainRate(grad_ds, dstrain_rate); 400 NewtonianStress(context, dstrain_rate, dkmstress); 401 KMUnpack(dkmstress, dstress); 402 KMUnpack(kmstress, stress); 403 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 404 405 StateConservative dF_inviscid[3]; 406 FluxInviscid_fwd(context, s, ds, dF_inviscid); 407 408 // Total flux 409 CeedScalar dFlux[5][3]; 410 FluxTotal(dF_inviscid, dstress, dFe, dFlux); 411 412 for (int j=0; j<3; j++) 413 for (int k=0; k<5; k++) 414 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * dFlux[k][0] + 415 dXdx[j][1] * dFlux[k][1] + 416 dXdx[j][2] * dFlux[k][2]); 417 418 const CeedScalar dbody_force[5] = {0, ds.U.density *g[0], ds.U.density *g[1], ds.U.density *g[2], 0}; 419 for (int j=0; j<5; j++) 420 v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]); 421 422 // -- Stabilization method: none (Galerkin), SU, or SUPG 423 CeedScalar dstab[5][3], U_dot[5] = {0}; 424 for (CeedInt j=0; j<5; j++) U_dot[j] = context->ijacobian_time_shift * dU[j]; 425 Stabilization(context, s, Tau_d, grad_ds, U_dot, dbody_force, x_i, dstab); 426 427 for (int j=0; j<5; j++) 428 for (int k=0; k<3; k++) 429 Grad_v[k][j][i] += wdetJ*(dstab[j][0] * dXdx[k][0] + 430 dstab[j][1] * dXdx[k][1] + 431 dstab[j][2] * dXdx[k][2]); 432 433 } // End Quadrature Point Loop 434 return 0; 435 } 436 437 // ***************************************************************************** 438 // Compute boundary integral (ie. for strongly set inflows) 439 // ***************************************************************************** 440 CEED_QFUNCTION_HELPER int BoundaryIntegral(void *ctx, CeedInt Q, 441 const CeedScalar *const *in, CeedScalar *const *out, 442 StateFromQi_t StateFromQi, StateFromQi_fwd_t StateFromQi_fwd) { 443 444 //*INDENT-OFF* 445 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 446 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 447 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 448 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 449 450 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA]) out[0], 451 (*jac_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA]) out[1]; 452 453 //*INDENT-ON* 454 455 const NewtonianIdealGasContext context = (NewtonianIdealGasContext) ctx; 456 const bool is_implicit = context->is_implicit; 457 458 CeedPragmaSIMD 459 for(CeedInt i=0; i<Q; i++) { 460 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 461 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 462 State s = StateFromQi(context, qi, x_i); 463 464 const CeedScalar wdetJb = (is_implicit ? -1. : 1.) * q_data_sur[0][i]; 465 // ---- Normal vector 466 const CeedScalar norm[3] = {q_data_sur[1][i], 467 q_data_sur[2][i], 468 q_data_sur[3][i] 469 }; 470 471 const CeedScalar dXdx[2][3] = { 472 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 473 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 474 }; 475 476 State grad_s[3]; 477 for (CeedInt j=0; j<3; j++) { 478 CeedScalar dx_i[3] = {0}, dqi[5]; 479 for (CeedInt k=0; k<5; k++) 480 dqi[k] = Grad_q[0][k][i] * dXdx[0][j] + 481 Grad_q[1][k][i] * dXdx[1][j]; 482 dx_i[j] = 1.; 483 grad_s[j] = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 484 } 485 486 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 487 KMStrainRate(grad_s, strain_rate); 488 NewtonianStress(context, strain_rate, kmstress); 489 KMUnpack(kmstress, stress); 490 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 491 492 StateConservative F_inviscid[3]; 493 FluxInviscid(context, s, F_inviscid); 494 495 CeedScalar Flux[5]; 496 FluxTotal_Boundary(F_inviscid, stress, Fe, norm, Flux); 497 498 for (CeedInt j=0; j<5; j++) v[j][i] = -wdetJb * Flux[j]; 499 500 for (int j=0; j<5; j++) jac_data_sur[j][i] = qi[j]; 501 for (int j=0; j<6; j++) jac_data_sur[5+j][i] = kmstress[j]; 502 } 503 return 0; 504 } 505 506 CEED_QFUNCTION(BoundaryIntegral_Conserv)(void *ctx, CeedInt Q, 507 const CeedScalar *const *in, CeedScalar *const *out) { 508 return BoundaryIntegral(ctx, Q, in, out, StateFromU, StateFromU_fwd); 509 } 510 511 CEED_QFUNCTION(BoundaryIntegral_Prim)(void *ctx, CeedInt Q, 512 const CeedScalar *const *in, CeedScalar *const *out) { 513 return BoundaryIntegral(ctx, Q, in, out, StateFromY, StateFromY_fwd); 514 } 515 516 // ***************************************************************************** 517 // Jacobian for "set nothing" boundary integral 518 // ***************************************************************************** 519 CEED_QFUNCTION_HELPER int BoundaryIntegral_Jacobian(void *ctx, CeedInt Q, 520 const CeedScalar *const *in, CeedScalar *const *out, 521 StateFromQi_t StateFromQi, StateFromQi_fwd_t StateFromQi_fwd) { 522 // *INDENT-OFF* 523 // Inputs 524 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 525 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 526 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 527 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 528 (*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 529 // Outputs 530 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 531 // *INDENT-ON* 532 533 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 534 const bool implicit = context->is_implicit; 535 536 CeedPragmaSIMD 537 // Quadrature Point Loop 538 for (CeedInt i=0; i<Q; i++) { 539 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 540 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 541 const CeedScalar norm[3] = {q_data_sur[1][i], 542 q_data_sur[2][i], 543 q_data_sur[3][i] 544 }; 545 const CeedScalar dXdx[2][3] = { 546 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 547 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 548 }; 549 550 CeedScalar qi[5], kmstress[6], dqi[5], dx_i[3] = {0.}; 551 for (int j=0; j<5; j++) qi[j] = jac_data_sur[j][i]; 552 for (int j=0; j<6; j++) kmstress[j] = jac_data_sur[5+j][i]; 553 for (int j=0; j<5; j++) dqi[j] = dq[j][i]; 554 555 State s = StateFromQi(context, qi, x_i); 556 State ds = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 557 558 State grad_ds[3]; 559 for (CeedInt j=0; j<3; j++) { 560 CeedScalar dx_i[3] = {0}, dqi_j[5]; 561 for (CeedInt k=0; k<5; k++) 562 dqi_j[k] = Grad_dq[0][k][i] * dXdx[0][j] + 563 Grad_dq[1][k][i] * dXdx[1][j]; 564 dx_i[j] = 1.; 565 grad_ds[j] = StateFromQi_fwd(context, s, dqi_j, x_i, dx_i); 566 } 567 568 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 569 KMStrainRate(grad_ds, dstrain_rate); 570 NewtonianStress(context, dstrain_rate, dkmstress); 571 KMUnpack(dkmstress, dstress); 572 KMUnpack(kmstress, stress); 573 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 574 575 StateConservative dF_inviscid[3]; 576 FluxInviscid_fwd(context, s, ds, dF_inviscid); 577 578 CeedScalar dFlux[5]; 579 FluxTotal_Boundary(dF_inviscid, dstress, dFe, norm, dFlux); 580 581 for (int j=0; j<5; j++) v[j][i] = -wdetJb * dFlux[j]; 582 } // End Quadrature Point Loop 583 return 0; 584 } 585 586 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Conserv)(void *ctx, CeedInt Q, 587 const CeedScalar *const *in, CeedScalar *const *out) { 588 return BoundaryIntegral_Jacobian(ctx, Q, in, out, StateFromU, StateFromU_fwd); 589 } 590 591 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Prim)(void *ctx, CeedInt Q, 592 const CeedScalar *const *in, CeedScalar *const *out) { 593 return BoundaryIntegral_Jacobian(ctx, Q, in, out, StateFromY, StateFromY_fwd); 594 } 595 596 // ***************************************************************************** 597 // Outflow boundary condition, weakly setting a constant pressure 598 // ***************************************************************************** 599 CEED_QFUNCTION_HELPER int PressureOutflow(void *ctx, CeedInt Q, 600 const CeedScalar *const *in, CeedScalar *const *out, 601 StateFromQi_t StateFromQi, StateFromQi_fwd_t StateFromQi_fwd) { 602 // *INDENT-OFF* 603 // Inputs 604 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 605 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 606 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 607 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 608 // Outputs 609 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 610 (*jac_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[1]; 611 // *INDENT-ON* 612 613 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 614 const bool implicit = context->is_implicit; 615 const CeedScalar P0 = context->P0; 616 617 CeedPragmaSIMD 618 // Quadrature Point Loop 619 for (CeedInt i=0; i<Q; i++) { 620 // Setup 621 // -- Interp in 622 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 623 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 624 State s = StateFromQi(context, qi, x_i); 625 s.Y.pressure = P0; 626 627 // -- Interp-to-Interp q_data 628 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 629 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 630 // We can effect this by swapping the sign on this weight 631 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 632 633 // ---- Normal vector 634 const CeedScalar norm[3] = {q_data_sur[1][i], q_data_sur[2][i], q_data_sur[3][i]}; 635 636 const CeedScalar dXdx[2][3] = { 637 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 638 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 639 }; 640 641 State grad_s[3]; 642 for (CeedInt j=0; j<3; j++) { 643 CeedScalar dx_i[3] = {0}, dqi[5]; 644 for (CeedInt k=0; k<5; k++) 645 dqi[k] = Grad_q[0][k][i] * dXdx[0][j] + 646 Grad_q[1][k][i] * dXdx[1][j]; 647 dx_i[j] = 1.; 648 grad_s[j] = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 649 } 650 651 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 652 KMStrainRate(grad_s, strain_rate); 653 NewtonianStress(context, strain_rate, kmstress); 654 KMUnpack(kmstress, stress); 655 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 656 657 StateConservative F_inviscid[3]; 658 FluxInviscid(context, s, F_inviscid); 659 660 CeedScalar Flux[5]; 661 FluxTotal_Boundary(F_inviscid, stress, Fe, norm, Flux); 662 663 for (CeedInt j=0; j<5; j++) v[j][i] = -wdetJb * Flux[j]; 664 665 // Save values for Jacobian 666 for (int j=0; j<5; j++) jac_data_sur[j][i] = qi[j]; 667 for (int j=0; j<6; j++) jac_data_sur[5+j][i] = kmstress[j]; 668 } // End Quadrature Point Loop 669 return 0; 670 } 671 672 CEED_QFUNCTION(PressureOutflow_Conserv)(void *ctx, CeedInt Q, 673 const CeedScalar *const *in, CeedScalar *const *out) { 674 return PressureOutflow(ctx, Q, in, out, StateFromU, StateFromU_fwd); 675 } 676 677 CEED_QFUNCTION(PressureOutflow_Prim)(void *ctx, CeedInt Q, 678 const CeedScalar *const *in, CeedScalar *const *out) { 679 return PressureOutflow(ctx, Q, in, out, StateFromY, StateFromY_fwd); 680 } 681 682 // ***************************************************************************** 683 // Jacobian for weak-pressure outflow boundary condition 684 // ***************************************************************************** 685 CEED_QFUNCTION_HELPER int PressureOutflow_Jacobian(void *ctx, CeedInt Q, 686 const CeedScalar *const *in, CeedScalar *const *out, 687 StateFromQi_t StateFromQi, StateFromQi_fwd_t StateFromQi_fwd) { 688 // *INDENT-OFF* 689 // Inputs 690 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 691 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 692 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 693 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 694 (*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 695 // Outputs 696 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 697 // *INDENT-ON* 698 699 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 700 const bool implicit = context->is_implicit; 701 702 CeedPragmaSIMD 703 // Quadrature Point Loop 704 for (CeedInt i=0; i<Q; i++) { 705 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 706 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 707 const CeedScalar norm[3] = {q_data_sur[1][i], q_data_sur[2][i], q_data_sur[3][i]}; 708 const CeedScalar dXdx[2][3] = { 709 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 710 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 711 }; 712 713 CeedScalar qi[5], kmstress[6], dqi[5], dx_i[3] = {0.}; 714 for (int j=0; j<5; j++) qi[j] = jac_data_sur[j][i]; 715 for (int j=0; j<6; j++) kmstress[j] = jac_data_sur[5+j][i]; 716 for (int j=0; j<5; j++) dqi[j] = dq[j][i]; 717 718 State s = StateFromQi(context, qi, x_i); 719 State ds = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 720 s.Y.pressure = context->P0; 721 ds.Y.pressure = 0.; 722 723 State grad_ds[3]; 724 for (CeedInt j=0; j<3; j++) { 725 CeedScalar dx_i[3] = {0}, dqi_j[5]; 726 for (CeedInt k=0; k<5; k++) 727 dqi_j[k] = Grad_dq[0][k][i] * dXdx[0][j] + 728 Grad_dq[1][k][i] * dXdx[1][j]; 729 dx_i[j] = 1.; 730 grad_ds[j] = StateFromQi_fwd(context, s, dqi_j, x_i, dx_i); 731 } 732 733 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 734 KMStrainRate(grad_ds, dstrain_rate); 735 NewtonianStress(context, dstrain_rate, dkmstress); 736 KMUnpack(dkmstress, dstress); 737 KMUnpack(kmstress, stress); 738 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 739 740 StateConservative dF_inviscid[3]; 741 FluxInviscid_fwd(context, s, ds, dF_inviscid); 742 743 CeedScalar dFlux[5]; 744 FluxTotal_Boundary(dF_inviscid, dstress, dFe, norm, dFlux); 745 746 for (int j=0; j<5; j++) v[j][i] = -wdetJb * dFlux[j]; 747 } // End Quadrature Point Loop 748 return 0; 749 } 750 751 CEED_QFUNCTION(PressureOutflow_Jacobian_Conserv)(void *ctx, CeedInt Q, 752 const CeedScalar *const *in, CeedScalar *const *out) { 753 return PressureOutflow_Jacobian(ctx, Q, in, out, StateFromU, StateFromU_fwd); 754 } 755 756 CEED_QFUNCTION(PressureOutflow_Jacobian_Prim)(void *ctx, CeedInt Q, 757 const CeedScalar *const *in, CeedScalar *const *out) { 758 return PressureOutflow_Jacobian(ctx, Q, in, out, StateFromY, StateFromY_fwd); 759 } 760 761 // ***************************************************************************** 762 // This QFunction implements the Navier-Stokes equations (mentioned above) in 763 // primitive variables and with implicit time stepping method 764 // 765 // ***************************************************************************** 766 CEED_QFUNCTION(IFunction_Newtonian_Prim)(void *ctx, CeedInt Q, 767 const CeedScalar *const *in, CeedScalar *const *out) { 768 // *INDENT-OFF* 769 // Inputs 770 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 771 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 772 (*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 773 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 774 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 775 // Outputs 776 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 777 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1], 778 (*jac_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[2]; 779 // *INDENT-ON* 780 // Context 781 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 782 const CeedScalar *g = context->g; 783 const CeedScalar dt = context->dt; 784 785 CeedPragmaSIMD 786 // Quadrature Point Loop 787 for (CeedInt i=0; i<Q; i++) { 788 CeedScalar Y[5]; 789 for (CeedInt j=0; j<5; j++) Y[j] = q[j][i]; 790 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 791 State s = StateFromY(context, Y, x_i); 792 793 // -- Interp-to-Interp q_data 794 const CeedScalar wdetJ = q_data[0][i]; 795 // -- Interp-to-Grad q_data 796 // ---- Inverse of change of coordinate matrix: X_i,j 797 // *INDENT-OFF* 798 const CeedScalar dXdx[3][3] = {{q_data[1][i], q_data[2][i], q_data[3][i]}, 799 {q_data[4][i], q_data[5][i], q_data[6][i]}, 800 {q_data[7][i], q_data[8][i], q_data[9][i]} 801 }; 802 // *INDENT-ON* 803 State grad_s[3]; 804 for (CeedInt j=0; j<3; j++) { 805 CeedScalar dx_i[3] = {0}, dY[5]; 806 for (CeedInt k=0; k<5; k++) 807 dY[k] = Grad_q[0][k][i] * dXdx[0][j] + 808 Grad_q[1][k][i] * dXdx[1][j] + 809 Grad_q[2][k][i] * dXdx[2][j]; 810 dx_i[j] = 1.; 811 grad_s[j] = StateFromY_fwd(context, s, dY, x_i, dx_i); 812 } 813 814 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 815 KMStrainRate(grad_s, strain_rate); 816 NewtonianStress(context, strain_rate, kmstress); 817 KMUnpack(kmstress, stress); 818 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 819 820 StateConservative F_inviscid[3]; 821 FluxInviscid(context, s, F_inviscid); 822 823 // Total flux 824 CeedScalar Flux[5][3]; 825 FluxTotal(F_inviscid, stress, Fe, Flux); 826 827 for (CeedInt j=0; j<3; j++) 828 for (CeedInt k=0; k<5; k++) 829 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * Flux[k][0] + 830 dXdx[j][1] * Flux[k][1] + 831 dXdx[j][2] * Flux[k][2]); 832 833 const CeedScalar body_force[5] = {0, s.U.density *g[0], s.U.density *g[1], s.U.density *g[2], 0}; 834 835 CeedScalar Y_dot[5], dx0[3] = {0}; 836 for (int j=0; j<5; j++) Y_dot[j] = q_dot[j][i]; 837 State s_dot = StateFromY_fwd(context, s, Y_dot, x_i, dx0); 838 839 CeedScalar U_dot[5] = {0.}; 840 UnpackState_U(s_dot.U, U_dot); 841 842 for (CeedInt j=0; j<5; j++) 843 v[j][i] = wdetJ * (U_dot[j] - body_force[j]); 844 845 // -- Stabilization method: none (Galerkin), SU, or SUPG 846 CeedScalar Tau_d[3], stab[5][3]; 847 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 848 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, x_i, stab); 849 850 for (CeedInt j=0; j<5; j++) 851 for (CeedInt k=0; k<3; k++) 852 Grad_v[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 853 stab[j][1] * dXdx[k][1] + 854 stab[j][2] * dXdx[k][2]); 855 856 for (CeedInt j=0; j<5; j++) jac_data[j][i] = Y[j]; 857 for (CeedInt j=0; j<6; j++) jac_data[5+j][i] = kmstress[j]; 858 for (CeedInt j=0; j<3; j++) jac_data[5+6+j][i] = Tau_d[j]; 859 860 } // End Quadrature Point Loop 861 862 // Return 863 return 0; 864 } 865 866 // ***************************************************************************** 867 // This QFunction implements the jacobean of the Navier-Stokes equations 868 // in primitive variables for implicit time stepping method. 869 // 870 // ***************************************************************************** 871 CEED_QFUNCTION(IJacobian_Newtonian_Prim)(void *ctx, CeedInt Q, 872 const CeedScalar *const *in, CeedScalar *const *out) { 873 // *INDENT-OFF* 874 // Inputs 875 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 876 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 877 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 878 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 879 (*jac_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 880 // Outputs 881 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 882 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 883 // *INDENT-ON* 884 // Context 885 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 886 const CeedScalar *g = context->g; 887 888 CeedPragmaSIMD 889 // Quadrature Point Loop 890 for (CeedInt i=0; i<Q; i++) { 891 // -- Interp-to-Interp q_data 892 const CeedScalar wdetJ = q_data[0][i]; 893 // -- Interp-to-Grad q_data 894 // ---- Inverse of change of coordinate matrix: X_i,j 895 // *INDENT-OFF* 896 const CeedScalar dXdx[3][3] = {{q_data[1][i], q_data[2][i], q_data[3][i]}, 897 {q_data[4][i], q_data[5][i], q_data[6][i]}, 898 {q_data[7][i], q_data[8][i], q_data[9][i]} 899 }; 900 // *INDENT-ON* 901 902 CeedScalar Y[5], kmstress[6], Tau_d[3] __attribute((unused)); 903 for (int j=0; j<5; j++) Y[j] = jac_data[j][i]; 904 for (int j=0; j<6; j++) kmstress[j] = jac_data[5+j][i]; 905 for (int j=0; j<3; j++) Tau_d[j] = jac_data[5+6+j][i]; 906 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 907 State s = StateFromY(context, Y, x_i); 908 909 CeedScalar dY[5], dx0[3] = {0}; 910 for (int j=0; j<5; j++) dY[j] = dq[j][i]; 911 State ds = StateFromY_fwd(context, s, dY, x_i, dx0); 912 913 State grad_ds[3]; 914 for (int j=0; j<3; j++) { 915 CeedScalar dYj[5]; 916 for (int k=0; k<5; k++) 917 dYj[k] = Grad_dq[0][k][i] * dXdx[0][j] + 918 Grad_dq[1][k][i] * dXdx[1][j] + 919 Grad_dq[2][k][i] * dXdx[2][j]; 920 grad_ds[j] = StateFromY_fwd(context, s, dYj, x_i, dx0); 921 } 922 923 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 924 KMStrainRate(grad_ds, dstrain_rate); 925 NewtonianStress(context, dstrain_rate, dkmstress); 926 KMUnpack(dkmstress, dstress); 927 KMUnpack(kmstress, stress); 928 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 929 930 StateConservative dF_inviscid[3]; 931 FluxInviscid_fwd(context, s, ds, dF_inviscid); 932 933 // Total flux 934 CeedScalar dFlux[5][3]; 935 FluxTotal(dF_inviscid, dstress, dFe, dFlux); 936 937 for (int j=0; j<3; j++) 938 for (int k=0; k<5; k++) 939 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * dFlux[k][0] + 940 dXdx[j][1] * dFlux[k][1] + 941 dXdx[j][2] * dFlux[k][2]); 942 943 const CeedScalar dbody_force[5] = {0, ds.U.density *g[0], ds.U.density *g[1], ds.U.density *g[2], 0}; 944 CeedScalar dU[5] = {0.}; 945 UnpackState_U(ds.U, dU); 946 947 for (int j=0; j<5; j++) 948 v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]); 949 950 // -- Stabilization method: none (Galerkin), SU, or SUPG 951 CeedScalar dstab[5][3], U_dot[5] = {0}; 952 for (CeedInt j=0; j<5; j++) U_dot[j] = context->ijacobian_time_shift * dU[j]; 953 Stabilization(context, s, Tau_d, grad_ds, U_dot, dbody_force, x_i, dstab); 954 955 for (int j=0; j<5; j++) 956 for (int k=0; k<3; k++) 957 Grad_v[k][j][i] += wdetJ*(dstab[j][0] * dXdx[k][0] + 958 dstab[j][1] * dXdx[k][1] + 959 dstab[j][2] * dXdx[k][2]); 960 961 } // End Quadrature Point Loop 962 return 0; 963 } 964 // ***************************************************************************** 965 966 #endif // newtonian_h 967