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 #include "stabilization.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], 183 q_data[2][i], 184 q_data[3][i]}, 185 {q_data[4][i], 186 q_data[5][i], 187 q_data[6][i]}, 188 {q_data[7][i], 189 q_data[8][i], 190 q_data[9][i]} 191 }; 192 // *INDENT-ON* 193 State grad_s[3]; 194 for (CeedInt j=0; j<3; j++) { 195 CeedScalar dx_i[3] = {0}, dU[5]; 196 for (CeedInt k=0; k<5; k++) 197 dU[k] = Grad_q[0][k][i] * dXdx[0][j] + 198 Grad_q[1][k][i] * dXdx[1][j] + 199 Grad_q[2][k][i] * dXdx[2][j]; 200 dx_i[j] = 1.; 201 grad_s[j] = StateFromU_fwd(context, s, dU, x_i, dx_i); 202 } 203 204 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 205 KMStrainRate(grad_s, strain_rate); 206 NewtonianStress(context, strain_rate, kmstress); 207 KMUnpack(kmstress, stress); 208 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 209 210 StateConservative F_inviscid[3]; 211 FluxInviscid(context, s, F_inviscid); 212 213 // Total flux 214 CeedScalar Flux[5][3]; 215 FluxTotal(F_inviscid, stress, Fe, Flux); 216 217 for (CeedInt j=0; j<3; j++) 218 for (CeedInt k=0; k<5; k++) 219 Grad_v[j][k][i] = wdetJ * (dXdx[j][0] * Flux[k][0] + 220 dXdx[j][1] * Flux[k][1] + 221 dXdx[j][2] * Flux[k][2]); 222 223 const CeedScalar body_force[5] = {0, s.U.density *g[0], s.U.density *g[1], s.U.density *g[2], 0}; 224 for (int j=0; j<5; j++) 225 v[j][i] = wdetJ * body_force[j]; 226 227 // -- Stabilization method: none (Galerkin), SU, or SUPG 228 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}; 229 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 230 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, x_i, stab); 231 232 for (CeedInt j=0; j<5; j++) 233 for (CeedInt k=0; k<3; k++) 234 Grad_v[k][j][i] -= wdetJ*(stab[j][0] * dXdx[k][0] + 235 stab[j][1] * dXdx[k][1] + 236 stab[j][2] * dXdx[k][2]); 237 238 } // End Quadrature Point Loop 239 240 // Return 241 return 0; 242 } 243 244 // ***************************************************************************** 245 // This QFunction implements the Navier-Stokes equations (mentioned above) with 246 // implicit time stepping method 247 // 248 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 249 // SUPG = Galerkin + grad(v) . ( Ai^T * Tau * (q_dot + Aj q,j - body force) ) 250 // (diffussive terms will be added later) 251 // 252 // ***************************************************************************** 253 CEED_QFUNCTION(IFunction_Newtonian)(void *ctx, CeedInt Q, 254 const CeedScalar *const *in, CeedScalar *const *out) { 255 // *INDENT-OFF* 256 // Inputs 257 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 258 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 259 (*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 260 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 261 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 262 // Outputs 263 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 264 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1], 265 (*jac_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[2]; 266 // *INDENT-ON* 267 // Context 268 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 269 const CeedScalar *g = context->g; 270 const CeedScalar dt = context->dt; 271 272 CeedPragmaSIMD 273 // Quadrature Point Loop 274 for (CeedInt i=0; i<Q; i++) { 275 CeedScalar U[5]; 276 for (CeedInt j=0; j<5; j++) U[j] = q[j][i]; 277 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 278 State s = StateFromU(context, U, x_i); 279 280 // -- Interp-to-Interp q_data 281 const CeedScalar wdetJ = q_data[0][i]; 282 // -- Interp-to-Grad q_data 283 // ---- Inverse of change of coordinate matrix: X_i,j 284 // *INDENT-OFF* 285 const CeedScalar dXdx[3][3] = {{q_data[1][i], 286 q_data[2][i], 287 q_data[3][i]}, 288 {q_data[4][i], 289 q_data[5][i], 290 q_data[6][i]}, 291 {q_data[7][i], 292 q_data[8][i], 293 q_data[9][i]} 294 }; 295 // *INDENT-ON* 296 State grad_s[3]; 297 for (CeedInt j=0; j<3; j++) { 298 CeedScalar dx_i[3] = {0}, dU[5]; 299 for (CeedInt k=0; k<5; k++) 300 dU[k] = Grad_q[0][k][i] * dXdx[0][j] + 301 Grad_q[1][k][i] * dXdx[1][j] + 302 Grad_q[2][k][i] * dXdx[2][j]; 303 dx_i[j] = 1.; 304 grad_s[j] = StateFromU_fwd(context, s, dU, x_i, dx_i); 305 } 306 307 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 308 KMStrainRate(grad_s, strain_rate); 309 NewtonianStress(context, strain_rate, kmstress); 310 KMUnpack(kmstress, stress); 311 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 312 313 StateConservative F_inviscid[3]; 314 FluxInviscid(context, s, F_inviscid); 315 316 // Total flux 317 CeedScalar Flux[5][3]; 318 FluxTotal(F_inviscid, stress, Fe, Flux); 319 320 for (CeedInt j=0; j<3; j++) 321 for (CeedInt k=0; k<5; k++) 322 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * Flux[k][0] + 323 dXdx[j][1] * Flux[k][1] + 324 dXdx[j][2] * Flux[k][2]); 325 326 const CeedScalar body_force[5] = {0, s.U.density *g[0], s.U.density *g[1], s.U.density *g[2], 0}; 327 for (CeedInt j=0; j<5; j++) 328 v[j][i] = wdetJ * (q_dot[j][i] - body_force[j]); 329 330 // -- Stabilization method: none (Galerkin), SU, or SUPG 331 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}; 332 for (CeedInt j=0; j<5; j++) U_dot[j] = q_dot[j][i]; 333 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 334 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, x_i, stab); 335 336 for (CeedInt j=0; j<5; j++) 337 for (CeedInt k=0; k<3; k++) 338 Grad_v[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 339 stab[j][1] * dXdx[k][1] + 340 stab[j][2] * dXdx[k][2]); 341 342 for (CeedInt j=0; j<5; j++) jac_data[j][i] = U[j]; 343 for (CeedInt j=0; j<6; j++) jac_data[5+j][i] = kmstress[j]; 344 for (CeedInt j=0; j<3; j++) jac_data[5+6+j][i] = Tau_d[j]; 345 346 } // End Quadrature Point Loop 347 348 // Return 349 return 0; 350 } 351 352 // ***************************************************************************** 353 // This QFunction implements the jacobean of the Navier-Stokes equations 354 // for implicit time stepping method. 355 // 356 // ***************************************************************************** 357 CEED_QFUNCTION(IJacobian_Newtonian)(void *ctx, CeedInt Q, 358 const CeedScalar *const *in, 359 CeedScalar *const *out) { 360 // *INDENT-OFF* 361 // Inputs 362 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 363 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 364 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 365 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 366 (*jac_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 367 // Outputs 368 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 369 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 370 // *INDENT-ON* 371 // Context 372 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 373 const CeedScalar *g = context->g; 374 375 CeedPragmaSIMD 376 // Quadrature Point Loop 377 for (CeedInt i=0; i<Q; i++) { 378 // -- Interp-to-Interp q_data 379 const CeedScalar wdetJ = q_data[0][i]; 380 // -- Interp-to-Grad q_data 381 // ---- Inverse of change of coordinate matrix: X_i,j 382 // *INDENT-OFF* 383 const CeedScalar dXdx[3][3] = {{q_data[1][i], 384 q_data[2][i], 385 q_data[3][i]}, 386 {q_data[4][i], 387 q_data[5][i], 388 q_data[6][i]}, 389 {q_data[7][i], 390 q_data[8][i], 391 q_data[9][i]} 392 }; 393 // *INDENT-ON* 394 395 CeedScalar U[5], kmstress[6], Tau_d[3] __attribute((unused)); 396 for (int j=0; j<5; j++) U[j] = jac_data[j][i]; 397 for (int j=0; j<6; j++) kmstress[j] = jac_data[5+j][i]; 398 for (int j=0; j<3; j++) Tau_d[j] = jac_data[5+6+j][i]; 399 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 400 State s = StateFromU(context, U, x_i); 401 402 CeedScalar dU[5], dx0[3] = {0}; 403 for (int j=0; j<5; j++) dU[j] = dq[j][i]; 404 State ds = StateFromU_fwd(context, s, dU, x_i, dx0); 405 406 State grad_ds[3]; 407 for (int j=0; j<3; j++) { 408 CeedScalar dUj[5]; 409 for (int k=0; k<5; k++) 410 dUj[k] = Grad_dq[0][k][i] * dXdx[0][j] + 411 Grad_dq[1][k][i] * dXdx[1][j] + 412 Grad_dq[2][k][i] * dXdx[2][j]; 413 grad_ds[j] = StateFromU_fwd(context, s, dUj, x_i, dx0); 414 } 415 416 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 417 KMStrainRate(grad_ds, dstrain_rate); 418 NewtonianStress(context, dstrain_rate, dkmstress); 419 KMUnpack(dkmstress, dstress); 420 KMUnpack(kmstress, stress); 421 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 422 423 StateConservative dF_inviscid[3]; 424 FluxInviscid_fwd(context, s, ds, dF_inviscid); 425 426 // Total flux 427 CeedScalar dFlux[5][3]; 428 FluxTotal(dF_inviscid, dstress, dFe, dFlux); 429 430 for (int j=0; j<3; j++) 431 for (int k=0; k<5; k++) 432 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * dFlux[k][0] + 433 dXdx[j][1] * dFlux[k][1] + 434 dXdx[j][2] * dFlux[k][2]); 435 436 const CeedScalar dbody_force[5] = {0, ds.U.density *g[0], ds.U.density *g[1], ds.U.density *g[2], 0}; 437 for (int j=0; j<5; j++) 438 v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]); 439 440 // -- Stabilization method: none (Galerkin), SU, or SUPG 441 CeedScalar dstab[5][3], U_dot[5] = {0}; 442 for (CeedInt j=0; j<5; j++) U_dot[j] = context->ijacobian_time_shift * dU[j]; 443 Stabilization(context, s, Tau_d, grad_ds, U_dot, dbody_force, x_i, dstab); 444 445 for (int j=0; j<5; j++) 446 for (int k=0; k<3; k++) 447 Grad_v[k][j][i] += wdetJ*(dstab[j][0] * dXdx[k][0] + 448 dstab[j][1] * dXdx[k][1] + 449 dstab[j][2] * dXdx[k][2]); 450 451 } // End Quadrature Point Loop 452 return 0; 453 } 454 455 // ***************************************************************************** 456 // Compute boundary integral (ie. for strongly set inflows) 457 // ***************************************************************************** 458 CEED_QFUNCTION(BoundaryIntegral)(void *ctx, CeedInt Q, 459 const CeedScalar *const *in, 460 CeedScalar *const *out) { 461 462 //*INDENT-OFF* 463 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 464 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 465 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 466 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 467 468 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA]) out[0], 469 (*jac_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA]) out[1]; 470 471 //*INDENT-ON* 472 473 const NewtonianIdealGasContext context = (NewtonianIdealGasContext) ctx; 474 const bool is_implicit = context->is_implicit; 475 State (*StateFromQi)(NewtonianIdealGasContext gas, 476 const CeedScalar qi[5], const CeedScalar x[3]); 477 State (*StateFromQi_fwd)(NewtonianIdealGasContext gas, 478 State s, const CeedScalar dqi[5], 479 const CeedScalar x[3], const CeedScalar dx[3]); 480 StateFromQi = context->is_primitive ? &StateFromY : &StateFromU; 481 StateFromQi_fwd = context->is_primitive ? &StateFromY_fwd : &StateFromU_fwd; 482 483 484 CeedPragmaSIMD 485 for(CeedInt i=0; i<Q; i++) { 486 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 487 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 488 State s = StateFromQi(context, qi, x_i); 489 490 const CeedScalar wdetJb = (is_implicit ? -1. : 1.) * q_data_sur[0][i]; 491 // ---- Normal vect 492 const CeedScalar norm[3] = {q_data_sur[1][i], 493 q_data_sur[2][i], 494 q_data_sur[3][i] 495 }; 496 497 const CeedScalar dXdx[2][3] = { 498 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 499 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 500 }; 501 502 State grad_s[3]; 503 for (CeedInt j=0; j<3; j++) { 504 CeedScalar dx_i[3] = {0}, dqi[5]; 505 for (CeedInt k=0; k<5; k++) 506 dqi[k] = Grad_q[0][k][i] * dXdx[0][j] + 507 Grad_q[1][k][i] * dXdx[1][j]; 508 dx_i[j] = 1.; 509 grad_s[j] = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 510 } 511 512 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 513 KMStrainRate(grad_s, strain_rate); 514 NewtonianStress(context, strain_rate, kmstress); 515 KMUnpack(kmstress, stress); 516 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 517 518 StateConservative F_inviscid[3]; 519 FluxInviscid(context, s, F_inviscid); 520 521 CeedScalar Flux[5] = {0.}; 522 for (int j=0; j<3; j++) { 523 Flux[0] += F_inviscid[j].density * norm[j]; 524 for (int k=0; k<3; k++) 525 Flux[k+1] += (F_inviscid[j].momentum[k] - stress[k][j]) * norm[j]; 526 Flux[4] += (F_inviscid[j].E_total + Fe[j]) * norm[j]; 527 } 528 529 // -- Density 530 v[0][i] = -wdetJb * Flux[0]; 531 532 // -- Momentum 533 for (CeedInt j=0; j<3; j++) 534 v[j+1][i] = -wdetJb * Flux[j+1]; 535 536 // -- Total Energy Density 537 v[4][i] = -wdetJb * Flux[4]; 538 539 if (context->is_primitive) { 540 jac_data_sur[0][i] = s.Y.pressure; 541 for (int j=0; j<3; j++) jac_data_sur[j+1][i] = s.Y.velocity[j]; 542 jac_data_sur[4][i] = s.Y.temperature; 543 } else { 544 jac_data_sur[0][i] = s.U.density; 545 for (int j=0; j<3; j++) jac_data_sur[j+1][i] = s.U.momentum[j]; 546 jac_data_sur[4][i] = s.U.E_total; 547 } 548 for (int j=0; j<6; j++) jac_data_sur[5+j][i] = kmstress[j]; 549 } 550 return 0; 551 } 552 553 // ***************************************************************************** 554 // Jacobian for "set nothing" boundary integral 555 // ***************************************************************************** 556 CEED_QFUNCTION(BoundaryIntegral_Jacobian)(void *ctx, CeedInt Q, 557 const CeedScalar *const *in, 558 CeedScalar *const *out) { 559 // *INDENT-OFF* 560 // Inputs 561 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 562 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 563 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 564 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 565 (*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 566 // Outputs 567 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 568 // *INDENT-ON* 569 570 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 571 const bool implicit = context->is_implicit; 572 State (*StateFromQi)(NewtonianIdealGasContext gas, 573 const CeedScalar qi[5], const CeedScalar x[3]); 574 State (*StateFromQi_fwd)(NewtonianIdealGasContext gas, 575 State s, const CeedScalar dqi[5], 576 const CeedScalar x[3], const CeedScalar dx[3]); 577 StateFromQi = context->is_primitive ? &StateFromY : &StateFromU; 578 StateFromQi_fwd = context->is_primitive ? &StateFromY_fwd : &StateFromU_fwd; 579 580 CeedPragmaSIMD 581 // Quadrature Point Loop 582 for (CeedInt i=0; i<Q; i++) { 583 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 584 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 585 const CeedScalar norm[3] = {q_data_sur[1][i], 586 q_data_sur[2][i], 587 q_data_sur[3][i] 588 }; 589 const CeedScalar dXdx[2][3] = { 590 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 591 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 592 }; 593 594 CeedScalar qi[5], kmstress[6], dqi[5], dx_i[3] = {0.}; 595 for (int j=0; j<5; j++) qi[j] = jac_data_sur[j][i]; 596 for (int j=0; j<6; j++) kmstress[j] = jac_data_sur[5+j][i]; 597 for (int j=0; j<5; j++) dqi[j] = dq[j][i]; 598 599 State s = StateFromQi(context, qi, x_i); 600 State ds = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 601 602 State grad_ds[3]; 603 for (CeedInt j=0; j<3; j++) { 604 CeedScalar dx_i[3] = {0}, dqi_j[5]; 605 for (CeedInt k=0; k<5; k++) 606 dqi_j[k] = Grad_dq[0][k][i] * dXdx[0][j] + 607 Grad_dq[1][k][i] * dXdx[1][j]; 608 dx_i[j] = 1.; 609 grad_ds[j] = StateFromQi_fwd(context, s, dqi_j, x_i, dx_i); 610 } 611 612 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 613 KMStrainRate(grad_ds, dstrain_rate); 614 NewtonianStress(context, dstrain_rate, dkmstress); 615 KMUnpack(dkmstress, dstress); 616 KMUnpack(kmstress, stress); 617 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 618 619 StateConservative dF_inviscid[3]; 620 FluxInviscid_fwd(context, s, ds, dF_inviscid); 621 622 CeedScalar dFlux[5] = {0.}; 623 for (int j=0; j<3; j++) { 624 dFlux[0] += dF_inviscid[j].density * norm[j]; 625 for (int k=0; k<3; k++) 626 dFlux[k+1] += (dF_inviscid[j].momentum[k] - dstress[k][j]) * norm[j]; 627 dFlux[4] += (dF_inviscid[j].E_total + dFe[j]) * norm[j]; 628 } 629 630 for (int j=0; j<5; j++) 631 v[j][i] = -wdetJb * dFlux[j]; 632 } // End Quadrature Point Loop 633 return 0; 634 } 635 636 // ***************************************************************************** 637 // Outflow boundary condition, weakly setting a constant pressure 638 // ***************************************************************************** 639 CEED_QFUNCTION(PressureOutflow)(void *ctx, CeedInt Q, 640 const CeedScalar *const *in, 641 CeedScalar *const *out) { 642 // *INDENT-OFF* 643 // Inputs 644 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 645 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 646 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 647 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 648 // Outputs 649 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 650 (*jac_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[1]; 651 // *INDENT-ON* 652 653 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 654 const bool implicit = context->is_implicit; 655 const CeedScalar P0 = context->P0; 656 657 State (*StateFromQi)(NewtonianIdealGasContext gas, 658 const CeedScalar qi[5], const CeedScalar x[3]); 659 State (*StateFromQi_fwd)(NewtonianIdealGasContext gas, 660 State s, const CeedScalar dqi[5], 661 const CeedScalar x[3], const CeedScalar dx[3]); 662 StateFromQi = context->is_primitive ? &StateFromY : &StateFromU; 663 StateFromQi_fwd = context->is_primitive ? &StateFromY_fwd : &StateFromU_fwd; 664 665 CeedPragmaSIMD 666 // Quadrature Point Loop 667 for (CeedInt i=0; i<Q; i++) { 668 // Setup 669 // -- Interp in 670 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 671 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 672 State s = StateFromQi(context, qi, x_i); 673 s.Y.pressure = P0; 674 675 // -- Interp-to-Interp q_data 676 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 677 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 678 // We can effect this by swapping the sign on this weight 679 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 680 681 // ---- Normal vect 682 const CeedScalar norm[3] = {q_data_sur[1][i], 683 q_data_sur[2][i], 684 q_data_sur[3][i] 685 }; 686 687 const CeedScalar dXdx[2][3] = { 688 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 689 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 690 }; 691 692 State grad_s[3]; 693 for (CeedInt j=0; j<3; j++) { 694 CeedScalar dx_i[3] = {0}, dqi[5]; 695 for (CeedInt k=0; k<5; k++) 696 dqi[k] = Grad_q[0][k][i] * dXdx[0][j] + 697 Grad_q[1][k][i] * dXdx[1][j]; 698 dx_i[j] = 1.; 699 grad_s[j] = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 700 } 701 702 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 703 KMStrainRate(grad_s, strain_rate); 704 NewtonianStress(context, strain_rate, kmstress); 705 KMUnpack(kmstress, stress); 706 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 707 708 StateConservative F_inviscid[3]; 709 FluxInviscid(context, s, F_inviscid); 710 711 CeedScalar Flux[5] = {0.}; 712 for (int j=0; j<3; j++) { 713 Flux[0] += F_inviscid[j].density * norm[j]; 714 for (int k=0; k<3; k++) 715 Flux[k+1] += (F_inviscid[j].momentum[k] - stress[k][j]) * norm[j]; 716 Flux[4] += (F_inviscid[j].E_total + Fe[j])*norm[j]; 717 } 718 719 // -- Density 720 v[0][i] = -wdetJb * Flux[0]; 721 722 // -- Momentum 723 for (CeedInt j=0; j<3; j++) 724 v[j+1][i] = -wdetJb * Flux[j+1]; 725 726 // -- Total Energy Density 727 v[4][i] = -wdetJb * Flux[4]; 728 729 // Save values for Jacobian 730 if (context->is_primitive) { 731 jac_data_sur[0][i] = s.Y.pressure; 732 for (int j=0; j<3; j++) jac_data_sur[j+1][i] = s.Y.velocity[j]; 733 jac_data_sur[4][i] = s.Y.temperature; 734 } else { 735 jac_data_sur[0][i] = s.U.density; 736 for (int j=0; j<3; j++) jac_data_sur[j+1][i] = s.U.momentum[j]; 737 jac_data_sur[4][i] = s.U.E_total; 738 } 739 for (int j=0; j<6; j++) jac_data_sur[5+j][i] = kmstress[j]; 740 } // End Quadrature Point Loop 741 return 0; 742 } 743 744 // ***************************************************************************** 745 // Jacobian for weak-pressure outflow boundary condition 746 // ***************************************************************************** 747 CEED_QFUNCTION(PressureOutflow_Jacobian)(void *ctx, CeedInt Q, 748 const CeedScalar *const *in, CeedScalar *const *out) { 749 // *INDENT-OFF* 750 // Inputs 751 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 752 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 753 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 754 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 755 (*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 756 // Outputs 757 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 758 // *INDENT-ON* 759 760 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 761 const bool implicit = context->is_implicit; 762 763 State (*StateFromQi)(NewtonianIdealGasContext gas, 764 const CeedScalar qi[5], const CeedScalar x[3]); 765 State (*StateFromQi_fwd)(NewtonianIdealGasContext gas, 766 State s, const CeedScalar dQi[5], 767 const CeedScalar x[3], const CeedScalar dx[3]); 768 StateFromQi = context->is_primitive ? &StateFromY : &StateFromU; 769 StateFromQi_fwd = context->is_primitive ? &StateFromY_fwd : &StateFromU_fwd; 770 771 CeedPragmaSIMD 772 // Quadrature Point Loop 773 for (CeedInt i=0; i<Q; i++) { 774 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 775 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 776 const CeedScalar norm[3] = {q_data_sur[1][i], 777 q_data_sur[2][i], 778 q_data_sur[3][i] 779 }; 780 const CeedScalar dXdx[2][3] = { 781 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 782 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 783 }; 784 785 CeedScalar qi[5], kmstress[6], dqi[5], dx_i[3] = {0.}; 786 for (int j=0; j<5; j++) qi[j] = jac_data_sur[j][i]; 787 for (int j=0; j<6; j++) kmstress[j] = jac_data_sur[5+j][i]; 788 for (int j=0; j<5; j++) dqi[j] = dq[j][i]; 789 790 State s = StateFromQi(context, qi, x_i); 791 State ds = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 792 s.Y.pressure = context->P0; 793 ds.Y.pressure = 0.; 794 795 State grad_ds[3]; 796 for (CeedInt j=0; j<3; j++) { 797 CeedScalar dx_i[3] = {0}, dqi_j[5]; 798 for (CeedInt k=0; k<5; k++) 799 dqi_j[k] = Grad_dq[0][k][i] * dXdx[0][j] + 800 Grad_dq[1][k][i] * dXdx[1][j]; 801 dx_i[j] = 1.; 802 grad_ds[j] = StateFromQi_fwd(context, s, dqi_j, x_i, dx_i); 803 } 804 805 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 806 KMStrainRate(grad_ds, dstrain_rate); 807 NewtonianStress(context, dstrain_rate, dkmstress); 808 KMUnpack(dkmstress, dstress); 809 KMUnpack(kmstress, stress); 810 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 811 812 StateConservative dF_inviscid[3]; 813 FluxInviscid_fwd(context, s, ds, dF_inviscid); 814 815 CeedScalar dFlux[5] = {0.}; 816 for (int j=0; j<3; j++) { 817 dFlux[0] += dF_inviscid[j].density * norm[j]; 818 for (int k=0; k<3; k++) 819 dFlux[k+1] += (dF_inviscid[j].momentum[k] - dstress[k][j]) * norm[j]; 820 dFlux[4] += (dF_inviscid[j].E_total + dFe[j]) * norm[j]; 821 } 822 823 for (int j=0; j<5; j++) 824 v[j][i] = -wdetJb * dFlux[j]; 825 } // End Quadrature Point Loop 826 return 0; 827 } 828 829 // ***************************************************************************** 830 // This QFunction implements the Navier-Stokes equations (mentioned above) in 831 // primitive variables and with implicit time stepping method 832 // 833 // ***************************************************************************** 834 CEED_QFUNCTION(IFunction_Newtonian_Prim)(void *ctx, CeedInt Q, 835 const CeedScalar *const *in, CeedScalar *const *out) { 836 // *INDENT-OFF* 837 // Inputs 838 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 839 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 840 (*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 841 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 842 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 843 // Outputs 844 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 845 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1], 846 (*jac_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[2]; 847 // *INDENT-ON* 848 // Context 849 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 850 const CeedScalar *g = context->g; 851 const CeedScalar dt = context->dt; 852 853 CeedPragmaSIMD 854 // Quadrature Point Loop 855 for (CeedInt i=0; i<Q; i++) { 856 CeedScalar Y[5]; 857 for (CeedInt j=0; j<5; j++) Y[j] = q[j][i]; 858 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 859 State s = StateFromY(context, Y, x_i); 860 861 // -- Interp-to-Interp q_data 862 const CeedScalar wdetJ = q_data[0][i]; 863 // -- Interp-to-Grad q_data 864 // ---- Inverse of change of coordinate matrix: X_i,j 865 // *INDENT-OFF* 866 const CeedScalar dXdx[3][3] = {{q_data[1][i], 867 q_data[2][i], 868 q_data[3][i]}, 869 {q_data[4][i], 870 q_data[5][i], 871 q_data[6][i]}, 872 {q_data[7][i], 873 q_data[8][i], 874 q_data[9][i]} 875 }; 876 // *INDENT-ON* 877 State grad_s[3]; 878 for (CeedInt j=0; j<3; j++) { 879 CeedScalar dx_i[3] = {0}, dY[5]; 880 for (CeedInt k=0; k<5; k++) 881 dY[k] = Grad_q[0][k][i] * dXdx[0][j] + 882 Grad_q[1][k][i] * dXdx[1][j] + 883 Grad_q[2][k][i] * dXdx[2][j]; 884 dx_i[j] = 1.; 885 grad_s[j] = StateFromY_fwd(context, s, dY, x_i, dx_i); 886 } 887 888 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 889 KMStrainRate(grad_s, strain_rate); 890 NewtonianStress(context, strain_rate, kmstress); 891 KMUnpack(kmstress, stress); 892 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 893 894 StateConservative F_inviscid[3]; 895 FluxInviscid(context, s, F_inviscid); 896 897 // Total flux 898 CeedScalar Flux[5][3]; 899 FluxTotal(F_inviscid, stress, Fe, Flux); 900 901 for (CeedInt j=0; j<3; j++) 902 for (CeedInt k=0; k<5; k++) 903 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * Flux[k][0] + 904 dXdx[j][1] * Flux[k][1] + 905 dXdx[j][2] * Flux[k][2]); 906 907 const CeedScalar body_force[5] = {0, s.U.density *g[0], s.U.density *g[1], s.U.density *g[2], 0}; 908 909 CeedScalar Y_dot[5], dx0[3] = {0}; 910 for (int j=0; j<5; j++) Y_dot[j] = q_dot[j][i]; 911 State s_dot = StateFromY_fwd(context, s, Y_dot, x_i, dx0); 912 913 CeedScalar U_dot[5] = {0.}; 914 UnpackState_U(s_dot.U, U_dot); 915 916 for (CeedInt j=0; j<5; j++) 917 v[j][i] = wdetJ * (U_dot[j] - body_force[j]); 918 919 // -- Stabilization method: none (Galerkin), SU, or SUPG 920 CeedScalar Tau_d[3], stab[5][3]; 921 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 922 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, x_i, stab); 923 924 for (CeedInt j=0; j<5; j++) 925 for (CeedInt k=0; k<3; k++) 926 Grad_v[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 927 stab[j][1] * dXdx[k][1] + 928 stab[j][2] * dXdx[k][2]); 929 930 for (CeedInt j=0; j<5; j++) jac_data[j][i] = Y[j]; 931 for (CeedInt j=0; j<6; j++) jac_data[5+j][i] = kmstress[j]; 932 for (CeedInt j=0; j<3; j++) jac_data[5+6+j][i] = Tau_d[j]; 933 934 } // End Quadrature Point Loop 935 936 // Return 937 return 0; 938 } 939 940 // ***************************************************************************** 941 // This QFunction implements the jacobean of the Navier-Stokes equations 942 // in primitive variables for implicit time stepping method. 943 // 944 // ***************************************************************************** 945 CEED_QFUNCTION(IJacobian_Newtonian_Prim)(void *ctx, CeedInt Q, 946 const CeedScalar *const *in, CeedScalar *const *out) { 947 // *INDENT-OFF* 948 // Inputs 949 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 950 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 951 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 952 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 953 (*jac_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 954 // Outputs 955 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 956 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 957 // *INDENT-ON* 958 // Context 959 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 960 const CeedScalar *g = context->g; 961 962 CeedPragmaSIMD 963 // Quadrature Point Loop 964 for (CeedInt i=0; i<Q; i++) { 965 // -- Interp-to-Interp q_data 966 const CeedScalar wdetJ = q_data[0][i]; 967 // -- Interp-to-Grad q_data 968 // ---- Inverse of change of coordinate matrix: X_i,j 969 // *INDENT-OFF* 970 const CeedScalar dXdx[3][3] = {{q_data[1][i], 971 q_data[2][i], 972 q_data[3][i]}, 973 {q_data[4][i], 974 q_data[5][i], 975 q_data[6][i]}, 976 {q_data[7][i], 977 q_data[8][i], 978 q_data[9][i]} 979 }; 980 // *INDENT-ON* 981 982 CeedScalar Y[5], kmstress[6], Tau_d[3] __attribute((unused)); 983 for (int j=0; j<5; j++) Y[j] = jac_data[j][i]; 984 for (int j=0; j<6; j++) kmstress[j] = jac_data[5+j][i]; 985 for (int j=0; j<3; j++) Tau_d[j] = jac_data[5+6+j][i]; 986 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 987 State s = StateFromY(context, Y, x_i); 988 989 CeedScalar dY[5], dx0[3] = {0}; 990 for (int j=0; j<5; j++) dY[j] = dq[j][i]; 991 State ds = StateFromY_fwd(context, s, dY, x_i, dx0); 992 993 State grad_ds[3]; 994 for (int j=0; j<3; j++) { 995 CeedScalar dYj[5]; 996 for (int k=0; k<5; k++) 997 dYj[k] = Grad_dq[0][k][i] * dXdx[0][j] + 998 Grad_dq[1][k][i] * dXdx[1][j] + 999 Grad_dq[2][k][i] * dXdx[2][j]; 1000 grad_ds[j] = StateFromY_fwd(context, s, dYj, x_i, dx0); 1001 } 1002 1003 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 1004 KMStrainRate(grad_ds, dstrain_rate); 1005 NewtonianStress(context, dstrain_rate, dkmstress); 1006 KMUnpack(dkmstress, dstress); 1007 KMUnpack(kmstress, stress); 1008 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 1009 1010 StateConservative dF_inviscid[3]; 1011 FluxInviscid_fwd(context, s, ds, dF_inviscid); 1012 1013 // Total flux 1014 CeedScalar dFlux[5][3]; 1015 FluxTotal(dF_inviscid, dstress, dFe, dFlux); 1016 1017 for (int j=0; j<3; j++) 1018 for (int k=0; k<5; k++) 1019 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * dFlux[k][0] + 1020 dXdx[j][1] * dFlux[k][1] + 1021 dXdx[j][2] * dFlux[k][2]); 1022 1023 const CeedScalar dbody_force[5] = {0, ds.U.density *g[0], ds.U.density *g[1], ds.U.density *g[2], 0}; 1024 CeedScalar dU[5] = {0.}; 1025 UnpackState_U(ds.U, dU); 1026 1027 for (int j=0; j<5; j++) 1028 v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]); 1029 1030 // -- Stabilization method: none (Galerkin), SU, or SUPG 1031 CeedScalar dstab[5][3], U_dot[5] = {0}; 1032 for (CeedInt j=0; j<5; j++) U_dot[j] = context->ijacobian_time_shift * dU[j]; 1033 Stabilization(context, s, Tau_d, grad_ds, U_dot, dbody_force, x_i, dstab); 1034 1035 for (int j=0; j<5; j++) 1036 for (int k=0; k<3; k++) 1037 Grad_v[k][j][i] += wdetJ*(dstab[j][0] * dXdx[k][0] + 1038 dstab[j][1] * dXdx[k][1] + 1039 dstab[j][2] * dXdx[k][2]); 1040 1041 } // End Quadrature Point Loop 1042 return 0; 1043 } 1044 // ***************************************************************************** 1045 1046 #endif // newtonian_h 1047