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