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(BoundaryIntegral)(void *ctx, CeedInt Q, 441 const CeedScalar *const *in, 442 CeedScalar *const *out) { 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 State (*StateFromQi)(NewtonianIdealGasContext gas, 458 const CeedScalar qi[5], const CeedScalar x[3]); 459 State (*StateFromQi_fwd)(NewtonianIdealGasContext gas, 460 State s, const CeedScalar dqi[5], 461 const CeedScalar x[3], const CeedScalar dx[3]); 462 StateFromQi = context->use_primitive ? &StateFromY : &StateFromU; 463 StateFromQi_fwd = context->use_primitive ? &StateFromY_fwd : &StateFromU_fwd; 464 465 466 CeedPragmaSIMD 467 for(CeedInt i=0; i<Q; i++) { 468 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 469 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 470 State s = StateFromQi(context, qi, x_i); 471 472 const CeedScalar wdetJb = (is_implicit ? -1. : 1.) * q_data_sur[0][i]; 473 // ---- Normal vect 474 const CeedScalar norm[3] = {q_data_sur[1][i], 475 q_data_sur[2][i], 476 q_data_sur[3][i] 477 }; 478 479 const CeedScalar dXdx[2][3] = { 480 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 481 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 482 }; 483 484 State grad_s[3]; 485 for (CeedInt j=0; j<3; j++) { 486 CeedScalar dx_i[3] = {0}, dqi[5]; 487 for (CeedInt k=0; k<5; k++) 488 dqi[k] = Grad_q[0][k][i] * dXdx[0][j] + 489 Grad_q[1][k][i] * dXdx[1][j]; 490 dx_i[j] = 1.; 491 grad_s[j] = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 492 } 493 494 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 495 KMStrainRate(grad_s, strain_rate); 496 NewtonianStress(context, strain_rate, kmstress); 497 KMUnpack(kmstress, stress); 498 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 499 500 StateConservative F_inviscid[3]; 501 FluxInviscid(context, s, F_inviscid); 502 503 CeedScalar Flux[5] = {0.}; 504 for (int j=0; j<3; j++) { 505 Flux[0] += F_inviscid[j].density * norm[j]; 506 for (int k=0; k<3; k++) 507 Flux[k+1] += (F_inviscid[j].momentum[k] - stress[k][j]) * norm[j]; 508 Flux[4] += (F_inviscid[j].E_total + Fe[j]) * norm[j]; 509 } 510 511 // -- Density 512 v[0][i] = -wdetJb * Flux[0]; 513 514 // -- Momentum 515 for (CeedInt j=0; j<3; j++) 516 v[j+1][i] = -wdetJb * Flux[j+1]; 517 518 // -- Total Energy Density 519 v[4][i] = -wdetJb * Flux[4]; 520 521 if (context->use_primitive) { 522 jac_data_sur[0][i] = s.Y.pressure; 523 for (int j=0; j<3; j++) jac_data_sur[j+1][i] = s.Y.velocity[j]; 524 jac_data_sur[4][i] = s.Y.temperature; 525 } else { 526 jac_data_sur[0][i] = s.U.density; 527 for (int j=0; j<3; j++) jac_data_sur[j+1][i] = s.U.momentum[j]; 528 jac_data_sur[4][i] = s.U.E_total; 529 } 530 for (int j=0; j<6; j++) jac_data_sur[5+j][i] = kmstress[j]; 531 } 532 return 0; 533 } 534 535 // ***************************************************************************** 536 // Jacobian for "set nothing" boundary integral 537 // ***************************************************************************** 538 CEED_QFUNCTION(BoundaryIntegral_Jacobian)(void *ctx, CeedInt Q, 539 const CeedScalar *const *in, 540 CeedScalar *const *out) { 541 // *INDENT-OFF* 542 // Inputs 543 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 544 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 545 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 546 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 547 (*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 548 // Outputs 549 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 550 // *INDENT-ON* 551 552 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 553 const bool implicit = context->is_implicit; 554 State (*StateFromQi)(NewtonianIdealGasContext gas, 555 const CeedScalar qi[5], const CeedScalar x[3]); 556 State (*StateFromQi_fwd)(NewtonianIdealGasContext gas, 557 State s, const CeedScalar dqi[5], 558 const CeedScalar x[3], const CeedScalar dx[3]); 559 StateFromQi = context->use_primitive ? &StateFromY : &StateFromU; 560 StateFromQi_fwd = context->use_primitive ? &StateFromY_fwd : &StateFromU_fwd; 561 562 CeedPragmaSIMD 563 // Quadrature Point Loop 564 for (CeedInt i=0; i<Q; i++) { 565 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 566 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 567 const CeedScalar norm[3] = {q_data_sur[1][i], 568 q_data_sur[2][i], 569 q_data_sur[3][i] 570 }; 571 const CeedScalar dXdx[2][3] = { 572 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 573 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 574 }; 575 576 CeedScalar qi[5], kmstress[6], dqi[5], dx_i[3] = {0.}; 577 for (int j=0; j<5; j++) qi[j] = jac_data_sur[j][i]; 578 for (int j=0; j<6; j++) kmstress[j] = jac_data_sur[5+j][i]; 579 for (int j=0; j<5; j++) dqi[j] = dq[j][i]; 580 581 State s = StateFromQi(context, qi, x_i); 582 State ds = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 583 584 State grad_ds[3]; 585 for (CeedInt j=0; j<3; j++) { 586 CeedScalar dx_i[3] = {0}, dqi_j[5]; 587 for (CeedInt k=0; k<5; k++) 588 dqi_j[k] = Grad_dq[0][k][i] * dXdx[0][j] + 589 Grad_dq[1][k][i] * dXdx[1][j]; 590 dx_i[j] = 1.; 591 grad_ds[j] = StateFromQi_fwd(context, s, dqi_j, x_i, dx_i); 592 } 593 594 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 595 KMStrainRate(grad_ds, dstrain_rate); 596 NewtonianStress(context, dstrain_rate, dkmstress); 597 KMUnpack(dkmstress, dstress); 598 KMUnpack(kmstress, stress); 599 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 600 601 StateConservative dF_inviscid[3]; 602 FluxInviscid_fwd(context, s, ds, dF_inviscid); 603 604 CeedScalar dFlux[5] = {0.}; 605 for (int j=0; j<3; j++) { 606 dFlux[0] += dF_inviscid[j].density * norm[j]; 607 for (int k=0; k<3; k++) 608 dFlux[k+1] += (dF_inviscid[j].momentum[k] - dstress[k][j]) * norm[j]; 609 dFlux[4] += (dF_inviscid[j].E_total + dFe[j]) * norm[j]; 610 } 611 612 for (int j=0; j<5; j++) 613 v[j][i] = -wdetJb * dFlux[j]; 614 } // End Quadrature Point Loop 615 return 0; 616 } 617 618 // ***************************************************************************** 619 // Outflow boundary condition, weakly setting a constant pressure 620 // ***************************************************************************** 621 CEED_QFUNCTION(PressureOutflow)(void *ctx, CeedInt Q, 622 const CeedScalar *const *in, 623 CeedScalar *const *out) { 624 // *INDENT-OFF* 625 // Inputs 626 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 627 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 628 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 629 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 630 // Outputs 631 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 632 (*jac_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[1]; 633 // *INDENT-ON* 634 635 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 636 const bool implicit = context->is_implicit; 637 const CeedScalar P0 = context->P0; 638 639 State (*StateFromQi)(NewtonianIdealGasContext gas, 640 const CeedScalar qi[5], const CeedScalar x[3]); 641 State (*StateFromQi_fwd)(NewtonianIdealGasContext gas, 642 State s, const CeedScalar dqi[5], 643 const CeedScalar x[3], const CeedScalar dx[3]); 644 StateFromQi = context->use_primitive ? &StateFromY : &StateFromU; 645 StateFromQi_fwd = context->use_primitive ? &StateFromY_fwd : &StateFromU_fwd; 646 647 CeedPragmaSIMD 648 // Quadrature Point Loop 649 for (CeedInt i=0; i<Q; i++) { 650 // Setup 651 // -- Interp in 652 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 653 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 654 State s = StateFromQi(context, qi, x_i); 655 s.Y.pressure = P0; 656 657 // -- Interp-to-Interp q_data 658 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 659 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 660 // We can effect this by swapping the sign on this weight 661 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 662 663 // ---- Normal vect 664 const CeedScalar norm[3] = {q_data_sur[1][i], 665 q_data_sur[2][i], 666 q_data_sur[3][i] 667 }; 668 669 const CeedScalar dXdx[2][3] = { 670 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 671 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 672 }; 673 674 State grad_s[3]; 675 for (CeedInt j=0; j<3; j++) { 676 CeedScalar dx_i[3] = {0}, dqi[5]; 677 for (CeedInt k=0; k<5; k++) 678 dqi[k] = Grad_q[0][k][i] * dXdx[0][j] + 679 Grad_q[1][k][i] * dXdx[1][j]; 680 dx_i[j] = 1.; 681 grad_s[j] = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 682 } 683 684 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 685 KMStrainRate(grad_s, strain_rate); 686 NewtonianStress(context, strain_rate, kmstress); 687 KMUnpack(kmstress, stress); 688 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 689 690 StateConservative F_inviscid[3]; 691 FluxInviscid(context, s, F_inviscid); 692 693 CeedScalar Flux[5] = {0.}; 694 for (int j=0; j<3; j++) { 695 Flux[0] += F_inviscid[j].density * norm[j]; 696 for (int k=0; k<3; k++) 697 Flux[k+1] += (F_inviscid[j].momentum[k] - stress[k][j]) * norm[j]; 698 Flux[4] += (F_inviscid[j].E_total + Fe[j])*norm[j]; 699 } 700 701 // -- Density 702 v[0][i] = -wdetJb * Flux[0]; 703 704 // -- Momentum 705 for (CeedInt j=0; j<3; j++) 706 v[j+1][i] = -wdetJb * Flux[j+1]; 707 708 // -- Total Energy Density 709 v[4][i] = -wdetJb * Flux[4]; 710 711 // Save values for Jacobian 712 if (context->use_primitive) { 713 jac_data_sur[0][i] = s.Y.pressure; 714 for (int j=0; j<3; j++) jac_data_sur[j+1][i] = s.Y.velocity[j]; 715 jac_data_sur[4][i] = s.Y.temperature; 716 } else { 717 jac_data_sur[0][i] = s.U.density; 718 for (int j=0; j<3; j++) jac_data_sur[j+1][i] = s.U.momentum[j]; 719 jac_data_sur[4][i] = s.U.E_total; 720 } 721 for (int j=0; j<6; j++) jac_data_sur[5+j][i] = kmstress[j]; 722 } // End Quadrature Point Loop 723 return 0; 724 } 725 726 // ***************************************************************************** 727 // Jacobian for weak-pressure outflow boundary condition 728 // ***************************************************************************** 729 CEED_QFUNCTION(PressureOutflow_Jacobian)(void *ctx, CeedInt Q, 730 const CeedScalar *const *in, CeedScalar *const *out) { 731 // *INDENT-OFF* 732 // Inputs 733 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 734 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 735 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 736 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 737 (*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 738 // Outputs 739 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 740 // *INDENT-ON* 741 742 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 743 const bool implicit = context->is_implicit; 744 745 State (*StateFromQi)(NewtonianIdealGasContext gas, 746 const CeedScalar qi[5], const CeedScalar x[3]); 747 State (*StateFromQi_fwd)(NewtonianIdealGasContext gas, 748 State s, const CeedScalar dQi[5], 749 const CeedScalar x[3], const CeedScalar dx[3]); 750 StateFromQi = context->use_primitive ? &StateFromY : &StateFromU; 751 StateFromQi_fwd = context->use_primitive ? &StateFromY_fwd : &StateFromU_fwd; 752 753 CeedPragmaSIMD 754 // Quadrature Point Loop 755 for (CeedInt i=0; i<Q; i++) { 756 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 757 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 758 const CeedScalar norm[3] = {q_data_sur[1][i], 759 q_data_sur[2][i], 760 q_data_sur[3][i] 761 }; 762 const CeedScalar dXdx[2][3] = { 763 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 764 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 765 }; 766 767 CeedScalar qi[5], kmstress[6], dqi[5], dx_i[3] = {0.}; 768 for (int j=0; j<5; j++) qi[j] = jac_data_sur[j][i]; 769 for (int j=0; j<6; j++) kmstress[j] = jac_data_sur[5+j][i]; 770 for (int j=0; j<5; j++) dqi[j] = dq[j][i]; 771 772 State s = StateFromQi(context, qi, x_i); 773 State ds = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 774 s.Y.pressure = context->P0; 775 ds.Y.pressure = 0.; 776 777 State grad_ds[3]; 778 for (CeedInt j=0; j<3; j++) { 779 CeedScalar dx_i[3] = {0}, dqi_j[5]; 780 for (CeedInt k=0; k<5; k++) 781 dqi_j[k] = Grad_dq[0][k][i] * dXdx[0][j] + 782 Grad_dq[1][k][i] * dXdx[1][j]; 783 dx_i[j] = 1.; 784 grad_ds[j] = StateFromQi_fwd(context, s, dqi_j, x_i, dx_i); 785 } 786 787 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 788 KMStrainRate(grad_ds, dstrain_rate); 789 NewtonianStress(context, dstrain_rate, dkmstress); 790 KMUnpack(dkmstress, dstress); 791 KMUnpack(kmstress, stress); 792 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 793 794 StateConservative dF_inviscid[3]; 795 FluxInviscid_fwd(context, s, ds, dF_inviscid); 796 797 CeedScalar dFlux[5] = {0.}; 798 for (int j=0; j<3; j++) { 799 dFlux[0] += dF_inviscid[j].density * norm[j]; 800 for (int k=0; k<3; k++) 801 dFlux[k+1] += (dF_inviscid[j].momentum[k] - dstress[k][j]) * norm[j]; 802 dFlux[4] += (dF_inviscid[j].E_total + dFe[j]) * norm[j]; 803 } 804 805 for (int j=0; j<5; j++) 806 v[j][i] = -wdetJb * dFlux[j]; 807 } // End Quadrature Point Loop 808 return 0; 809 } 810 811 // ***************************************************************************** 812 // This QFunction implements the Navier-Stokes equations (mentioned above) in 813 // primitive variables and with implicit time stepping method 814 // 815 // ***************************************************************************** 816 CEED_QFUNCTION(IFunction_Newtonian_Prim)(void *ctx, CeedInt Q, 817 const CeedScalar *const *in, CeedScalar *const *out) { 818 // *INDENT-OFF* 819 // Inputs 820 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 821 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 822 (*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 823 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 824 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 825 // Outputs 826 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 827 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1], 828 (*jac_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[2]; 829 // *INDENT-ON* 830 // Context 831 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 832 const CeedScalar *g = context->g; 833 const CeedScalar dt = context->dt; 834 835 CeedPragmaSIMD 836 // Quadrature Point Loop 837 for (CeedInt i=0; i<Q; i++) { 838 CeedScalar Y[5]; 839 for (CeedInt j=0; j<5; j++) Y[j] = q[j][i]; 840 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 841 State s = StateFromY(context, Y, x_i); 842 843 // -- Interp-to-Interp q_data 844 const CeedScalar wdetJ = q_data[0][i]; 845 // -- Interp-to-Grad q_data 846 // ---- Inverse of change of coordinate matrix: X_i,j 847 // *INDENT-OFF* 848 const CeedScalar dXdx[3][3] = {{q_data[1][i], q_data[2][i], q_data[3][i]}, 849 {q_data[4][i], q_data[5][i], q_data[6][i]}, 850 {q_data[7][i], q_data[8][i], q_data[9][i]} 851 }; 852 // *INDENT-ON* 853 State grad_s[3]; 854 for (CeedInt j=0; j<3; j++) { 855 CeedScalar dx_i[3] = {0}, dY[5]; 856 for (CeedInt k=0; k<5; k++) 857 dY[k] = Grad_q[0][k][i] * dXdx[0][j] + 858 Grad_q[1][k][i] * dXdx[1][j] + 859 Grad_q[2][k][i] * dXdx[2][j]; 860 dx_i[j] = 1.; 861 grad_s[j] = StateFromY_fwd(context, s, dY, x_i, dx_i); 862 } 863 864 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 865 KMStrainRate(grad_s, strain_rate); 866 NewtonianStress(context, strain_rate, kmstress); 867 KMUnpack(kmstress, stress); 868 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 869 870 StateConservative F_inviscid[3]; 871 FluxInviscid(context, s, F_inviscid); 872 873 // Total flux 874 CeedScalar Flux[5][3]; 875 FluxTotal(F_inviscid, stress, Fe, Flux); 876 877 for (CeedInt j=0; j<3; j++) 878 for (CeedInt k=0; k<5; k++) 879 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * Flux[k][0] + 880 dXdx[j][1] * Flux[k][1] + 881 dXdx[j][2] * Flux[k][2]); 882 883 const CeedScalar body_force[5] = {0, s.U.density *g[0], s.U.density *g[1], s.U.density *g[2], 0}; 884 885 CeedScalar Y_dot[5], dx0[3] = {0}; 886 for (int j=0; j<5; j++) Y_dot[j] = q_dot[j][i]; 887 State s_dot = StateFromY_fwd(context, s, Y_dot, x_i, dx0); 888 889 CeedScalar U_dot[5] = {0.}; 890 UnpackState_U(s_dot.U, U_dot); 891 892 for (CeedInt j=0; j<5; j++) 893 v[j][i] = wdetJ * (U_dot[j] - body_force[j]); 894 895 // -- Stabilization method: none (Galerkin), SU, or SUPG 896 CeedScalar Tau_d[3], stab[5][3]; 897 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 898 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, x_i, stab); 899 900 for (CeedInt j=0; j<5; j++) 901 for (CeedInt k=0; k<3; k++) 902 Grad_v[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 903 stab[j][1] * dXdx[k][1] + 904 stab[j][2] * dXdx[k][2]); 905 906 for (CeedInt j=0; j<5; j++) jac_data[j][i] = Y[j]; 907 for (CeedInt j=0; j<6; j++) jac_data[5+j][i] = kmstress[j]; 908 for (CeedInt j=0; j<3; j++) jac_data[5+6+j][i] = Tau_d[j]; 909 910 } // End Quadrature Point Loop 911 912 // Return 913 return 0; 914 } 915 916 // ***************************************************************************** 917 // This QFunction implements the jacobean of the Navier-Stokes equations 918 // in primitive variables for implicit time stepping method. 919 // 920 // ***************************************************************************** 921 CEED_QFUNCTION(IJacobian_Newtonian_Prim)(void *ctx, CeedInt Q, 922 const CeedScalar *const *in, CeedScalar *const *out) { 923 // *INDENT-OFF* 924 // Inputs 925 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 926 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 927 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 928 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 929 (*jac_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 930 // Outputs 931 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 932 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 933 // *INDENT-ON* 934 // Context 935 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 936 const CeedScalar *g = context->g; 937 938 CeedPragmaSIMD 939 // Quadrature Point Loop 940 for (CeedInt i=0; i<Q; i++) { 941 // -- Interp-to-Interp q_data 942 const CeedScalar wdetJ = q_data[0][i]; 943 // -- Interp-to-Grad q_data 944 // ---- Inverse of change of coordinate matrix: X_i,j 945 // *INDENT-OFF* 946 const CeedScalar dXdx[3][3] = {{q_data[1][i], q_data[2][i], q_data[3][i]}, 947 {q_data[4][i], q_data[5][i], q_data[6][i]}, 948 {q_data[7][i], q_data[8][i], q_data[9][i]} 949 }; 950 // *INDENT-ON* 951 952 CeedScalar Y[5], kmstress[6], Tau_d[3] __attribute((unused)); 953 for (int j=0; j<5; j++) Y[j] = jac_data[j][i]; 954 for (int j=0; j<6; j++) kmstress[j] = jac_data[5+j][i]; 955 for (int j=0; j<3; j++) Tau_d[j] = jac_data[5+6+j][i]; 956 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 957 State s = StateFromY(context, Y, x_i); 958 959 CeedScalar dY[5], dx0[3] = {0}; 960 for (int j=0; j<5; j++) dY[j] = dq[j][i]; 961 State ds = StateFromY_fwd(context, s, dY, x_i, dx0); 962 963 State grad_ds[3]; 964 for (int j=0; j<3; j++) { 965 CeedScalar dYj[5]; 966 for (int k=0; k<5; k++) 967 dYj[k] = Grad_dq[0][k][i] * dXdx[0][j] + 968 Grad_dq[1][k][i] * dXdx[1][j] + 969 Grad_dq[2][k][i] * dXdx[2][j]; 970 grad_ds[j] = StateFromY_fwd(context, s, dYj, x_i, dx0); 971 } 972 973 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 974 KMStrainRate(grad_ds, dstrain_rate); 975 NewtonianStress(context, dstrain_rate, dkmstress); 976 KMUnpack(dkmstress, dstress); 977 KMUnpack(kmstress, stress); 978 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 979 980 StateConservative dF_inviscid[3]; 981 FluxInviscid_fwd(context, s, ds, dF_inviscid); 982 983 // Total flux 984 CeedScalar dFlux[5][3]; 985 FluxTotal(dF_inviscid, dstress, dFe, dFlux); 986 987 for (int j=0; j<3; j++) 988 for (int k=0; k<5; k++) 989 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * dFlux[k][0] + 990 dXdx[j][1] * dFlux[k][1] + 991 dXdx[j][2] * dFlux[k][2]); 992 993 const CeedScalar dbody_force[5] = {0, ds.U.density *g[0], ds.U.density *g[1], ds.U.density *g[2], 0}; 994 CeedScalar dU[5] = {0.}; 995 UnpackState_U(ds.U, dU); 996 997 for (int j=0; j<5; j++) 998 v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]); 999 1000 // -- Stabilization method: none (Galerkin), SU, or SUPG 1001 CeedScalar dstab[5][3], U_dot[5] = {0}; 1002 for (CeedInt j=0; j<5; j++) U_dot[j] = context->ijacobian_time_shift * dU[j]; 1003 Stabilization(context, s, Tau_d, grad_ds, U_dot, dbody_force, x_i, dstab); 1004 1005 for (int j=0; j<5; j++) 1006 for (int k=0; k<3; k++) 1007 Grad_v[k][j][i] += wdetJ*(dstab[j][0] * dXdx[k][0] + 1008 dstab[j][1] * dXdx[k][1] + 1009 dstab[j][2] * dXdx[k][2]); 1010 1011 } // End Quadrature Point Loop 1012 return 0; 1013 } 1014 // ***************************************************************************** 1015 1016 #endif // newtonian_h 1017