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 20 #ifndef M_PI 21 #define M_PI 3.14159265358979323846 22 #endif 23 24 // ***************************************************************************** 25 // Helper function for computing flux Jacobian 26 // ***************************************************************************** 27 CEED_QFUNCTION_HELPER void computeFluxJacobian_NS(CeedScalar dF[3][5][5], 28 const CeedScalar rho, const CeedScalar u[3], const CeedScalar E, 29 const CeedScalar gamma, const CeedScalar g[3], const CeedScalar x[3]) { 30 CeedScalar u_sq = u[0]*u[0] + u[1]*u[1] + u[2]*u[2]; // Velocity square 31 CeedScalar e_potential = -(g[0]*x[0] + g[1]*x[1] + g[2]*x[2]); 32 for (CeedInt i=0; i<3; i++) { // Jacobian matrices for 3 directions 33 for (CeedInt j=0; j<3; j++) { // Rows of each Jacobian matrix 34 dF[i][j+1][0] = ((i==j) ? ((gamma-1.)*(u_sq/2. - e_potential)) : 0.) - 35 u[i]*u[j]; 36 for (CeedInt k=0; k<3; k++) { // Columns of each Jacobian matrix 37 dF[i][0][k+1] = ((i==k) ? 1. : 0.); 38 dF[i][j+1][k+1] = ((j==k) ? u[i] : 0.) + 39 ((i==k) ? u[j] : 0.) - 40 ((i==j) ? u[k] : 0.) * (gamma-1.); 41 dF[i][4][k+1] = ((i==k) ? (E*gamma/rho - (gamma-1.)*u_sq/2.) : 0.) - 42 (gamma-1.)*u[i]*u[k]; 43 } 44 dF[i][j+1][4] = ((i==j) ? (gamma-1.) : 0.); 45 } 46 dF[i][4][0] = u[i] * ((gamma-1.)*u_sq - E*gamma/rho); 47 dF[i][4][4] = u[i] * gamma; 48 } 49 } 50 51 // ***************************************************************************** 52 // Helper function for computing flux Jacobian of Primitive variables 53 // ***************************************************************************** 54 CEED_QFUNCTION_HELPER void computeFluxJacobian_NSp(CeedScalar dF[3][5][5], 55 const CeedScalar rho, const CeedScalar u[3], const CeedScalar E, 56 const CeedScalar Rd, const CeedScalar cv) { 57 CeedScalar u_sq = u[0]*u[0] + u[1]*u[1] + u[2]*u[2]; // Velocity square 58 // TODO Add in gravity's contribution 59 60 CeedScalar T = ( E / rho - u_sq / 2. ) / cv; 61 CeedScalar drdT = -rho / T; 62 CeedScalar drdP = 1. / ( Rd * T); 63 CeedScalar etot = E / rho ; 64 CeedScalar e2p = drdP * etot + 1. ; 65 CeedScalar e3p = ( E + rho * Rd * T ); 66 CeedScalar e4p = drdT * etot + rho * cv ; 67 68 for (CeedInt i=0; i<3; i++) { // Jacobian matrices for 3 directions 69 for (CeedInt j=0; j<3; j++) { // j counts F^{m_j} 70 // [row][col] of A_i 71 dF[i][j+1][0] = drdP * u[i] * u[j] + ((i==j) ? 1. : 0.); // F^{{m_j} wrt p 72 for (CeedInt k=0; k<3; k++) { // k counts the wrt vel_k 73 dF[i][0][k+1] = ((i==k) ? rho : 0.); // F^c wrt u_k 74 dF[i][j+1][k+1] = (((j==k) ? u[i] : 0.) + // F^m_j wrt u_k 75 ((i==k) ? u[j] : 0.) ) * rho; 76 dF[i][4][k+1] = rho * u[i] * u[k] 77 + ((i==k) ? e3p : 0.) ; // F^e wrt u_k 78 } 79 dF[i][j+1][4] = drdT * u[i] * u[j]; // F^{m_j} wrt T 80 } 81 dF[i][4][0] = u[i] * e2p; // F^e wrt p 82 dF[i][4][4] = u[i] * e4p; // F^e wrt T 83 dF[i][0][0] = u[i] * drdP; // F^c wrt p 84 dF[i][0][4] = u[i] * drdT; // F^c wrt T 85 } 86 } 87 88 CEED_QFUNCTION_HELPER void PrimitiveToConservative_fwd(const CeedScalar rho, 89 const CeedScalar u[3], const CeedScalar E, const CeedScalar Rd, 90 const CeedScalar cv, const CeedScalar dY[5], CeedScalar dU[5]) { 91 CeedScalar u_sq = u[0]*u[0] + u[1]*u[1] + u[2]*u[2]; 92 CeedScalar T = ( E / rho - u_sq / 2. ) / cv; 93 CeedScalar drdT = -rho / T; 94 CeedScalar drdP = 1. / ( Rd * T); 95 dU[0] = drdP * dY[0] + drdT * dY[4]; 96 CeedScalar de_kinetic = 0; 97 for (CeedInt i=0; i<3; i++) { 98 dU[1+i] = dU[0] * u[i] + rho * dY[1+i]; 99 de_kinetic += u[i] * dY[1+i]; 100 } 101 dU[4] = rho * cv * dY[4] + dU[0] * cv * T // internal energy: rho * e 102 + rho * de_kinetic + .5 * dU[0] * u_sq; // kinetic energy: .5 * rho * |u|^2 103 } 104 105 // ***************************************************************************** 106 // Helper function for computing Tau elements (stabilization constant) 107 // Model from: 108 // PHASTA 109 // 110 // Tau[i] = itau=0 which is diagonal-Shakib (3 values still but not spatial) 111 // 112 // Where NOT UPDATED YET 113 // ***************************************************************************** 114 CEED_QFUNCTION_HELPER void Tau_diagPrim(CeedScalar Tau_d[3], 115 const CeedScalar dXdx[3][3], const CeedScalar u[3], 116 const CeedScalar cv, const NewtonianIdealGasContext newt_ctx, 117 const CeedScalar mu, const CeedScalar dt, 118 const CeedScalar rho) { 119 // Context 120 const CeedScalar Ctau_t = newt_ctx->Ctau_t; 121 const CeedScalar Ctau_v = newt_ctx->Ctau_v; 122 const CeedScalar Ctau_C = newt_ctx->Ctau_C; 123 const CeedScalar Ctau_M = newt_ctx->Ctau_M; 124 const CeedScalar Ctau_E = newt_ctx->Ctau_E; 125 CeedScalar gijd[6]; 126 CeedScalar tau; 127 CeedScalar dts; 128 CeedScalar fact; 129 130 //*INDENT-OFF* 131 gijd[0] = dXdx[0][0] * dXdx[0][0] 132 + dXdx[1][0] * dXdx[1][0] 133 + dXdx[2][0] * dXdx[2][0]; 134 135 gijd[1] = dXdx[0][0] * dXdx[0][1] 136 + dXdx[1][0] * dXdx[1][1] 137 + dXdx[2][0] * dXdx[2][1]; 138 139 gijd[2] = dXdx[0][1] * dXdx[0][1] 140 + dXdx[1][1] * dXdx[1][1] 141 + dXdx[2][1] * dXdx[2][1]; 142 143 gijd[3] = dXdx[0][0] * dXdx[0][2] 144 + dXdx[1][0] * dXdx[1][2] 145 + dXdx[2][0] * dXdx[2][2]; 146 147 gijd[4] = dXdx[0][1] * dXdx[0][2] 148 + dXdx[1][1] * dXdx[1][2] 149 + dXdx[2][1] * dXdx[2][2]; 150 151 gijd[5] = dXdx[0][2] * dXdx[0][2] 152 + dXdx[1][2] * dXdx[1][2] 153 + dXdx[2][2] * dXdx[2][2]; 154 //*INDENT-ON* 155 156 dts = Ctau_t / dt ; 157 158 tau = rho*rho*((4. * dts * dts) 159 + u[0] * ( u[0] * gijd[0] + 2. * ( u[1] * gijd[1] + u[2] * gijd[3])) 160 + u[1] * ( u[1] * gijd[2] + 2. * u[2] * gijd[4]) 161 + u[2] * u[2] * gijd[5]) 162 + Ctau_v* mu * mu * 163 (gijd[0]*gijd[0] + gijd[2]*gijd[2] + gijd[5]*gijd[5] + 164 + 2. * (gijd[1]*gijd[1] + gijd[3]*gijd[3] + gijd[4]*gijd[4])); 165 166 fact=sqrt(tau); 167 168 Tau_d[0] = Ctau_C * fact / (rho*(gijd[0] + gijd[2] + gijd[5]))*0.125; 169 170 Tau_d[1] = Ctau_M / fact; 171 Tau_d[2] = Ctau_E / ( fact * cv ); 172 173 // consider putting back the way I initially had it Ctau_E * Tau_d[1] /cv 174 // to avoid a division if the compiler is smart enough to see that cv IS 175 // a constant that it could invert once for all elements 176 // but in that case energy tau is scaled by the product of Ctau_E * Ctau_M 177 // OR we could absorb cv into Ctau_E but this puts more burden on user to 178 // know how to change constants with a change of fluid or units. Same for 179 // Ctau_v * mu * mu IF AND ONLY IF we don't add viscosity law =f(T) 180 } 181 182 // ***************************************************************************** 183 // This QFunction sets a "still" initial condition for generic Newtonian IG problems 184 // ***************************************************************************** 185 CEED_QFUNCTION(ICsNewtonianIG)(void *ctx, CeedInt Q, 186 const CeedScalar *const *in, CeedScalar *const *out) { 187 // Inputs 188 const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 189 190 // Outputs 191 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 192 193 // Context 194 const SetupContext context = (SetupContext)ctx; 195 const CeedScalar theta0 = context->theta0; 196 const CeedScalar P0 = context->P0; 197 const CeedScalar cv = context->cv; 198 const CeedScalar cp = context->cp; 199 const CeedScalar *g = context->g; 200 const CeedScalar Rd = cp - cv; 201 202 // Quadrature Point Loop 203 CeedPragmaSIMD 204 for (CeedInt i=0; i<Q; i++) { 205 CeedScalar q[5] = {0.}; 206 207 // Setup 208 // -- Coordinates 209 const CeedScalar x[3] = {X[0][i], X[1][i], X[2][i]}; 210 const CeedScalar e_potential = -(g[0]*x[0] + g[1]*x[1] + g[2]*x[2]); 211 212 // -- Density 213 const CeedScalar rho = P0 / (Rd*theta0); 214 215 // Initial Conditions 216 q[0] = rho; 217 q[1] = 0.0; 218 q[2] = 0.0; 219 q[3] = 0.0; 220 q[4] = rho * (cv*theta0 + e_potential); 221 222 for (CeedInt j=0; j<5; j++) 223 q0[j][i] = q[j]; 224 } // End of Quadrature Point Loop 225 return 0; 226 } 227 228 // ***************************************************************************** 229 // This QFunction implements the following formulation of Navier-Stokes with 230 // explicit time stepping method 231 // 232 // This is 3D compressible Navier-Stokes in conservation form with state 233 // variables of density, momentum density, and total energy density. 234 // 235 // State Variables: q = ( rho, U1, U2, U3, E ) 236 // rho - Mass Density 237 // Ui - Momentum Density, Ui = rho ui 238 // E - Total Energy Density, E = rho (cv T + (u u)/2 + g z) 239 // 240 // Navier-Stokes Equations: 241 // drho/dt + div( U ) = 0 242 // dU/dt + div( rho (u x u) + P I3 ) + rho g khat = div( Fu ) 243 // dE/dt + div( (E + P) u ) = div( Fe ) 244 // 245 // Viscous Stress: 246 // Fu = mu (grad( u ) + grad( u )^T + lambda div ( u ) I3) 247 // 248 // Thermal Stress: 249 // Fe = u Fu + k grad( T ) 250 // Equation of State 251 // P = (gamma - 1) (E - rho (u u) / 2 - rho g z) 252 // 253 // Stabilization: 254 // Tau = diag(TauC, TauM, TauM, TauM, TauE) 255 // f1 = rho sqrt(ui uj gij) 256 // gij = dXi/dX * dXi/dX 257 // TauC = Cc f1 / (8 gii) 258 // TauM = min( 1 , 1 / f1 ) 259 // TauE = TauM / (Ce cv) 260 // 261 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 262 // 263 // Constants: 264 // lambda = - 2 / 3, From Stokes hypothesis 265 // mu , Dynamic viscosity 266 // k , Thermal conductivity 267 // cv , Specific heat, constant volume 268 // cp , Specific heat, constant pressure 269 // g , Gravity 270 // gamma = cp / cv, Specific heat ratio 271 // 272 // We require the product of the inverse of the Jacobian (dXdx_j,k) and 273 // its transpose (dXdx_k,j) to properly compute integrals of the form: 274 // int( gradv gradu ) 275 // 276 // ***************************************************************************** 277 CEED_QFUNCTION(RHSFunction_Newtonian)(void *ctx, CeedInt Q, 278 const CeedScalar *const *in, CeedScalar *const *out) { 279 // *INDENT-OFF* 280 // Inputs 281 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 282 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 283 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 284 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 285 // Outputs 286 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 287 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 288 // *INDENT-ON* 289 290 // Context 291 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 292 const CeedScalar mu = context->mu; 293 const CeedScalar cv = context->cv; 294 const CeedScalar cp = context->cp; 295 const CeedScalar *g = context->g; 296 const CeedScalar dt = context->dt; 297 const CeedScalar gamma = cp / cv; 298 const CeedScalar Rd = cp - cv; 299 300 CeedPragmaSIMD 301 // Quadrature Point Loop 302 for (CeedInt i=0; i<Q; i++) { 303 CeedScalar U[5]; 304 for (int j=0; j<5; j++) U[j] = q[j][i]; 305 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 306 State s = StateFromU(context, U, x_i); 307 308 // -- Interp-to-Interp q_data 309 const CeedScalar wdetJ = q_data[0][i]; 310 // -- Interp-to-Grad q_data 311 // ---- Inverse of change of coordinate matrix: X_i,j 312 // *INDENT-OFF* 313 const CeedScalar dXdx[3][3] = {{q_data[1][i], 314 q_data[2][i], 315 q_data[3][i]}, 316 {q_data[4][i], 317 q_data[5][i], 318 q_data[6][i]}, 319 {q_data[7][i], 320 q_data[8][i], 321 q_data[9][i]} 322 }; 323 // *INDENT-ON* 324 325 State grad_s[3]; 326 for (CeedInt j=0; j<3; j++) { 327 CeedScalar dx_i[3] = {0}, dU[5]; 328 for (CeedInt k=0; k<5; k++) 329 dU[k] = Grad_q[0][k][i] * dXdx[0][j] + 330 Grad_q[1][k][i] * dXdx[1][j] + 331 Grad_q[2][k][i] * dXdx[2][j]; 332 dx_i[j] = 1.; 333 grad_s[j] = StateFromU_fwd(context, s, dU, x_i, dx_i); 334 } 335 336 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 337 KMStrainRate(grad_s, strain_rate); 338 NewtonianStress(context, strain_rate, kmstress); 339 KMUnpack(kmstress, stress); 340 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 341 342 StateConservative F_inviscid[3]; 343 FluxInviscid(context, s, F_inviscid); 344 345 // Total flux 346 CeedScalar Flux[5][3]; 347 for (CeedInt j=0; j<3; j++) { 348 Flux[0][j] = F_inviscid[j].density; 349 for (CeedInt k=0; k<3; k++) 350 Flux[k+1][j] = F_inviscid[j].momentum[k] - stress[k][j]; 351 Flux[4][j] = F_inviscid[j].E_total + Fe[j]; 352 } 353 354 for (CeedInt j=0; j<3; j++) { 355 for (CeedInt k=0; k<5; k++) { 356 Grad_v[j][k][i] = wdetJ * (dXdx[j][0] * Flux[k][0] + 357 dXdx[j][1] * Flux[k][1] + 358 dXdx[j][2] * Flux[k][2]); 359 } 360 } 361 362 const CeedScalar body_force[5] = {0, s.U.density *g[0], s.U.density *g[1], s.U.density *g[2], 0}; 363 for (int j=0; j<5; j++) 364 v[j][i] = wdetJ * body_force[j]; 365 366 // jacob_F_conv[3][5][5] = dF(convective)/dq at each direction 367 CeedScalar jacob_F_conv[3][5][5] = {0}; 368 computeFluxJacobian_NS(jacob_F_conv, s.U.density, s.Y.velocity, s.U.E_total, 369 gamma, g, x_i); 370 CeedScalar grad_U[5][3]; 371 for (CeedInt j=0; j<3; j++) { 372 grad_U[0][j] = grad_s[j].U.density; 373 for (CeedInt k=0; k<3; k++) grad_U[k+1][j] = grad_s[j].U.momentum[k]; 374 grad_U[4][j] = grad_s[j].U.E_total; 375 } 376 377 // strong_conv = dF/dq * dq/dx (Strong convection) 378 CeedScalar strong_conv[5] = {0}; 379 for (CeedInt j=0; j<3; j++) 380 for (CeedInt k=0; k<5; k++) 381 for (CeedInt l=0; l<5; l++) 382 strong_conv[k] += jacob_F_conv[j][k][l] * grad_U[l][j]; 383 384 // -- Stabilization method: none, SU, or SUPG 385 CeedScalar stab[5][3] = {{0.}}; 386 CeedScalar tau_strong_conv[5] = {0.}, tau_strong_conv_conservative[5] = {0}; 387 CeedScalar Tau_d[3] = {0.}; 388 switch (context->stabilization) { 389 case STAB_NONE: // Galerkin 390 break; 391 case STAB_SU: // SU 392 Tau_diagPrim(Tau_d, dXdx, s.Y.velocity, cv, context, mu, dt, s.U.density); 393 tau_strong_conv[0] = Tau_d[0] * strong_conv[0]; 394 tau_strong_conv[1] = Tau_d[1] * strong_conv[1]; 395 tau_strong_conv[2] = Tau_d[1] * strong_conv[2]; 396 tau_strong_conv[3] = Tau_d[1] * strong_conv[3]; 397 tau_strong_conv[4] = Tau_d[2] * strong_conv[4]; 398 PrimitiveToConservative_fwd(s.U.density, s.Y.velocity, s.U.E_total, Rd, cv, 399 tau_strong_conv, 400 tau_strong_conv_conservative); 401 for (CeedInt j=0; j<3; j++) 402 for (CeedInt k=0; k<5; k++) 403 for (CeedInt l=0; l<5; l++) 404 stab[k][j] += jacob_F_conv[j][k][l] * tau_strong_conv_conservative[l]; 405 406 for (CeedInt j=0; j<5; j++) 407 for (CeedInt k=0; k<3; k++) 408 Grad_v[k][j][i] -= wdetJ*(stab[j][0] * dXdx[k][0] + 409 stab[j][1] * dXdx[k][1] + 410 stab[j][2] * dXdx[k][2]); 411 break; 412 case STAB_SUPG: // SUPG is not implemented for explicit scheme 413 break; 414 } 415 416 } // End Quadrature Point Loop 417 418 // Return 419 return 0; 420 } 421 422 // ***************************************************************************** 423 // This QFunction implements the Navier-Stokes equations (mentioned above) with 424 // implicit time stepping method 425 // 426 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 427 // SUPG = Galerkin + grad(v) . ( Ai^T * Tau * (q_dot + Aj q,j - body force) ) 428 // (diffussive terms will be added later) 429 // 430 // ***************************************************************************** 431 CEED_QFUNCTION(IFunction_Newtonian)(void *ctx, CeedInt Q, 432 const CeedScalar *const *in, 433 CeedScalar *const *out) { 434 // *INDENT-OFF* 435 // Inputs 436 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 437 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 438 (*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 439 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 440 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 441 // Outputs 442 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 443 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1], 444 (*jac_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[2]; 445 // *INDENT-ON* 446 // Context 447 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 448 const CeedScalar mu = context->mu; 449 const CeedScalar cv = context->cv; 450 const CeedScalar cp = context->cp; 451 const CeedScalar *g = context->g; 452 const CeedScalar dt = context->dt; 453 const CeedScalar gamma = cp / cv; 454 const CeedScalar Rd = cp-cv; 455 456 CeedPragmaSIMD 457 // Quadrature Point Loop 458 for (CeedInt i=0; i<Q; i++) { 459 CeedScalar U[5]; 460 for (CeedInt j=0; j<5; j++) U[j] = q[j][i]; 461 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 462 State s = StateFromU(context, U, x_i); 463 464 // -- Interp-to-Interp q_data 465 const CeedScalar wdetJ = q_data[0][i]; 466 // -- Interp-to-Grad q_data 467 // ---- Inverse of change of coordinate matrix: X_i,j 468 // *INDENT-OFF* 469 const CeedScalar dXdx[3][3] = {{q_data[1][i], 470 q_data[2][i], 471 q_data[3][i]}, 472 {q_data[4][i], 473 q_data[5][i], 474 q_data[6][i]}, 475 {q_data[7][i], 476 q_data[8][i], 477 q_data[9][i]} 478 }; 479 // *INDENT-ON* 480 State grad_s[3]; 481 for (CeedInt j=0; j<3; j++) { 482 CeedScalar dx_i[3] = {0}, dU[5]; 483 for (CeedInt k=0; k<5; k++) 484 dU[k] = Grad_q[0][k][i] * dXdx[0][j] + 485 Grad_q[1][k][i] * dXdx[1][j] + 486 Grad_q[2][k][i] * dXdx[2][j]; 487 dx_i[j] = 1.; 488 grad_s[j] = StateFromU_fwd(context, s, dU, x_i, dx_i); 489 } 490 491 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 492 KMStrainRate(grad_s, strain_rate); 493 NewtonianStress(context, strain_rate, kmstress); 494 KMUnpack(kmstress, stress); 495 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 496 497 StateConservative F_inviscid[3]; 498 FluxInviscid(context, s, F_inviscid); 499 500 501 // Total flux 502 CeedScalar Flux[5][3]; 503 for (CeedInt j=0; j<3; j++) { 504 Flux[0][j] = F_inviscid[j].density; 505 for (CeedInt k=0; k<3; k++) 506 Flux[k+1][j] = F_inviscid[j].momentum[k] - stress[k][j]; 507 Flux[4][j] = F_inviscid[j].E_total + Fe[j]; 508 } 509 510 for (CeedInt j=0; j<3; j++) { 511 for (CeedInt k=0; k<5; k++) { 512 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * Flux[k][0] + 513 dXdx[j][1] * Flux[k][1] + 514 dXdx[j][2] * Flux[k][2]); 515 } 516 } 517 518 const CeedScalar body_force[5] = {0, s.U.density *g[0], s.U.density *g[1], s.U.density *g[2], 0}; 519 for (CeedInt j=0; j<5; j++) 520 v[j][i] = wdetJ * (q_dot[j][i] - body_force[j]); 521 522 // jacob_F_conv[3][5][5] = dF(convective)/dq at each direction 523 CeedScalar jacob_F_conv[3][5][5] = {0}; 524 computeFluxJacobian_NS(jacob_F_conv, s.U.density, s.Y.velocity, s.U.E_total, 525 gamma, g, x_i); 526 CeedScalar grad_U[5][3]; 527 for (CeedInt j=0; j<3; j++) { 528 grad_U[0][j] = grad_s[j].U.density; 529 for (CeedInt k=0; k<3; k++) grad_U[k+1][j] = grad_s[j].U.momentum[k]; 530 grad_U[4][j] = grad_s[j].U.E_total; 531 } 532 533 // strong_conv = dF/dq * dq/dx (Strong convection) 534 CeedScalar strong_conv[5] = {0}; 535 for (CeedInt j=0; j<3; j++) 536 for (CeedInt k=0; k<5; k++) 537 for (CeedInt l=0; l<5; l++) 538 strong_conv[k] += jacob_F_conv[j][k][l] * grad_U[l][j]; 539 540 // Strong residual 541 CeedScalar strong_res[5]; 542 for (CeedInt j=0; j<5; j++) 543 strong_res[j] = q_dot[j][i] + strong_conv[j] - body_force[j]; 544 545 // -- Stabilization method: none, SU, or SUPG 546 CeedScalar stab[5][3] = {{0.}}; 547 CeedScalar tau_strong_res[5] = {0.}, tau_strong_res_conservative[5] = {0}; 548 CeedScalar tau_strong_conv[5] = {0.}, tau_strong_conv_conservative[5] = {0}; 549 CeedScalar Tau_d[3] = {0.}; 550 switch (context->stabilization) { 551 case STAB_NONE: // Galerkin 552 break; 553 case STAB_SU: // SU 554 Tau_diagPrim(Tau_d, dXdx, s.Y.velocity, cv, context, mu, dt, s.U.density); 555 tau_strong_conv[0] = Tau_d[0] * strong_conv[0]; 556 tau_strong_conv[1] = Tau_d[1] * strong_conv[1]; 557 tau_strong_conv[2] = Tau_d[1] * strong_conv[2]; 558 tau_strong_conv[3] = Tau_d[1] * strong_conv[3]; 559 tau_strong_conv[4] = Tau_d[2] * strong_conv[4]; 560 PrimitiveToConservative_fwd(s.U.density, s.Y.velocity, s.U.E_total, Rd, cv, 561 tau_strong_conv, tau_strong_conv_conservative); 562 for (CeedInt j=0; j<3; j++) 563 for (CeedInt k=0; k<5; k++) 564 for (CeedInt l=0; l<5; l++) 565 stab[k][j] += jacob_F_conv[j][k][l] * tau_strong_conv_conservative[l]; 566 567 for (CeedInt j=0; j<5; j++) 568 for (CeedInt k=0; k<3; k++) 569 Grad_v[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 570 stab[j][1] * dXdx[k][1] + 571 stab[j][2] * dXdx[k][2]); 572 573 break; 574 case STAB_SUPG: // SUPG 575 Tau_diagPrim(Tau_d, dXdx, s.Y.velocity, cv, context, mu, dt, s.U.density); 576 tau_strong_res[0] = Tau_d[0] * strong_res[0]; 577 tau_strong_res[1] = Tau_d[1] * strong_res[1]; 578 tau_strong_res[2] = Tau_d[1] * strong_res[2]; 579 tau_strong_res[3] = Tau_d[1] * strong_res[3]; 580 tau_strong_res[4] = Tau_d[2] * strong_res[4]; 581 // Alternate route (useful later with primitive variable code) 582 // this function was verified against PHASTA for as IC that was as close as possible 583 // computeFluxJacobian_NSp(jacob_F_conv_p, rho, u, E, Rd, cv); 584 // it has also been verified to compute a correct through the following 585 // stab[k][j] += jacob_F_conv_p[j][k][l] * tau_strong_res[l] // flux Jacobian wrt primitive 586 // applied in the triple loop below 587 // However, it is more flops than using the existing Jacobian wrt q after q_{,Y} viz 588 PrimitiveToConservative_fwd(s.U.density, s.Y.velocity, s.U.E_total, Rd, cv, 589 tau_strong_res, tau_strong_res_conservative); 590 for (CeedInt j=0; j<3; j++) 591 for (CeedInt k=0; k<5; k++) 592 for (CeedInt l=0; l<5; l++) 593 stab[k][j] += jacob_F_conv[j][k][l] * tau_strong_res_conservative[l]; 594 595 for (CeedInt j=0; j<5; j++) 596 for (CeedInt k=0; k<3; k++) 597 Grad_v[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 598 stab[j][1] * dXdx[k][1] + 599 stab[j][2] * dXdx[k][2]); 600 break; 601 } 602 for (CeedInt j=0; j<5; j++) jac_data[j][i] = U[j]; 603 for (CeedInt j=0; j<6; j++) jac_data[5+j][i] = kmstress[j]; 604 for (CeedInt j=0; j<3; j++) jac_data[5+6+j][i] = Tau_d[j]; 605 606 } // End Quadrature Point Loop 607 608 // Return 609 return 0; 610 } 611 612 CEED_QFUNCTION(IJacobian_Newtonian)(void *ctx, CeedInt Q, 613 const CeedScalar *const *in, 614 CeedScalar *const *out) { 615 // *INDENT-OFF* 616 // Inputs 617 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 618 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 619 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 620 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 621 (*jac_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 622 // Outputs 623 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 624 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 625 // *INDENT-ON* 626 // Context 627 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 628 const CeedScalar *g = context->g; 629 const CeedScalar cp = context->cp; 630 const CeedScalar cv = context->cv; 631 const CeedScalar Rd = cp - cv; 632 const CeedScalar gamma = cp / cv; 633 634 CeedPragmaSIMD 635 // Quadrature Point Loop 636 for (CeedInt i=0; i<Q; i++) { 637 // -- Interp-to-Interp q_data 638 const CeedScalar wdetJ = q_data[0][i]; 639 // -- Interp-to-Grad q_data 640 // ---- Inverse of change of coordinate matrix: X_i,j 641 // *INDENT-OFF* 642 const CeedScalar dXdx[3][3] = {{q_data[1][i], 643 q_data[2][i], 644 q_data[3][i]}, 645 {q_data[4][i], 646 q_data[5][i], 647 q_data[6][i]}, 648 {q_data[7][i], 649 q_data[8][i], 650 q_data[9][i]} 651 }; 652 // *INDENT-ON* 653 654 CeedScalar U[5], kmstress[6], Tau_d[3] __attribute((unused)); 655 for (int j=0; j<5; j++) U[j] = jac_data[j][i]; 656 for (int j=0; j<6; j++) kmstress[j] = jac_data[5+j][i]; 657 for (int j=0; j<3; j++) Tau_d[j] = jac_data[5+6+j][i]; 658 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 659 State s = StateFromU(context, U, x_i); 660 661 CeedScalar dU[5], dx0[3] = {0}; 662 for (int j=0; j<5; j++) dU[j] = dq[j][i]; 663 State ds = StateFromU_fwd(context, s, dU, x_i, dx0); 664 665 State grad_ds[3]; 666 for (int j=0; j<3; j++) { 667 CeedScalar dUj[5]; 668 for (int k=0; k<5; k++) dUj[k] = Grad_dq[0][k][i] * dXdx[0][j] 669 + Grad_dq[1][k][i] * dXdx[1][j] 670 + Grad_dq[2][k][i] * dXdx[2][j]; 671 grad_ds[j] = StateFromU_fwd(context, s, dUj, x_i, dx0); 672 } 673 674 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 675 KMStrainRate(grad_ds, dstrain_rate); 676 NewtonianStress(context, dstrain_rate, dkmstress); 677 KMUnpack(dkmstress, dstress); 678 KMUnpack(kmstress, stress); 679 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 680 681 StateConservative dF_inviscid[3]; 682 FluxInviscid_fwd(context, s, ds, dF_inviscid); 683 684 // Total flux 685 CeedScalar dFlux[5][3]; 686 for (int j=0; j<3; j++) { 687 dFlux[0][j] = dF_inviscid[j].density; 688 for (int k=0; k<3; k++) 689 dFlux[k+1][j] = dF_inviscid[j].momentum[k] - dstress[k][j]; 690 dFlux[4][j] = dF_inviscid[j].E_total + dFe[j]; 691 } 692 693 for (int j=0; j<3; j++) { 694 for (int k=0; k<5; k++) { 695 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * dFlux[k][0] + 696 dXdx[j][1] * dFlux[k][1] + 697 dXdx[j][2] * dFlux[k][2]); 698 } 699 } 700 701 const CeedScalar dbody_force[5] = {0, ds.U.density *g[0], ds.U.density *g[1], ds.U.density *g[2], 0}; 702 for (int j=0; j<5; j++) 703 v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]); 704 705 if (1) { 706 CeedScalar jacob_F_conv[3][5][5] = {0}; 707 computeFluxJacobian_NS(jacob_F_conv, s.U.density, s.Y.velocity, s.U.E_total, 708 gamma, g, x_i); 709 CeedScalar grad_dU[5][3]; 710 for (int j=0; j<3; j++) { 711 grad_dU[0][j] = grad_ds[j].U.density; 712 for (int k=0; k<3; k++) grad_dU[k+1][j] = grad_ds[j].U.momentum[k]; 713 grad_dU[4][j] = grad_ds[j].U.E_total; 714 } 715 CeedScalar dstrong_conv[5] = {0}; 716 for (int j=0; j<3; j++) 717 for (int k=0; k<5; k++) 718 for (int l=0; l<5; l++) 719 dstrong_conv[k] += jacob_F_conv[j][k][l] * grad_dU[l][j]; 720 CeedScalar dstrong_res[5]; 721 for (int j=0; j<5; j++) 722 dstrong_res[j] = context->ijacobian_time_shift * dU[j] + dstrong_conv[j] - 723 dbody_force[j]; 724 CeedScalar dtau_strong_res[5] = {0.}, dtau_strong_res_conservative[5] = {0}; 725 dtau_strong_res[0] = Tau_d[0] * dstrong_res[0]; 726 dtau_strong_res[1] = Tau_d[1] * dstrong_res[1]; 727 dtau_strong_res[2] = Tau_d[1] * dstrong_res[2]; 728 dtau_strong_res[3] = Tau_d[1] * dstrong_res[3]; 729 dtau_strong_res[4] = Tau_d[2] * dstrong_res[4]; 730 PrimitiveToConservative_fwd(s.U.density, s.Y.velocity, s.U.E_total, Rd, cv, 731 dtau_strong_res, dtau_strong_res_conservative); 732 CeedScalar dstab[5][3] = {0}; 733 for (int j=0; j<3; j++) 734 for (int k=0; k<5; k++) 735 for (int l=0; l<5; l++) 736 dstab[k][j] += jacob_F_conv[j][k][l] * dtau_strong_res_conservative[l]; 737 for (int j=0; j<5; j++) 738 for (int k=0; k<3; k++) 739 Grad_v[k][j][i] += wdetJ*(dstab[j][0] * dXdx[k][0] + 740 dstab[j][1] * dXdx[k][1] + 741 dstab[j][2] * dXdx[k][2]); 742 743 } 744 } // End Quadrature Point Loop 745 return 0; 746 } 747 748 // Compute boundary integral (ie. for strongly set inflows) 749 CEED_QFUNCTION(BoundaryIntegral)(void *ctx, CeedInt Q, 750 const CeedScalar *const *in, 751 CeedScalar *const *out) { 752 753 //*INDENT-OFF* 754 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 755 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 756 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 757 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 758 759 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA]) out[0], 760 (*jac_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA]) out[1]; 761 762 //*INDENT-ON* 763 764 const NewtonianIdealGasContext context = (NewtonianIdealGasContext) ctx; 765 const bool is_implicit = context->is_implicit; 766 767 CeedPragmaSIMD 768 for(CeedInt i=0; i<Q; i++) { 769 const CeedScalar U[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 770 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 771 const State s = StateFromU(context, U, x_i); 772 773 const CeedScalar wdetJb = (is_implicit ? -1. : 1.) * q_data_sur[0][i]; 774 // ---- Normal vect 775 const CeedScalar norm[3] = {q_data_sur[1][i], 776 q_data_sur[2][i], 777 q_data_sur[3][i] 778 }; 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 State grad_s[3]; 786 for (CeedInt j=0; j<3; j++) { 787 CeedScalar dx_i[3] = {0}, dU[5]; 788 for (CeedInt k=0; k<5; k++) 789 dU[k] = Grad_q[0][k][i] * dXdx[0][j] + 790 Grad_q[1][k][i] * dXdx[1][j]; 791 dx_i[j] = 1.; 792 grad_s[j] = StateFromU_fwd(context, s, dU, x_i, dx_i); 793 } 794 795 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 796 KMStrainRate(grad_s, strain_rate); 797 NewtonianStress(context, strain_rate, kmstress); 798 KMUnpack(kmstress, stress); 799 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 800 801 StateConservative F_inviscid[3]; 802 FluxInviscid(context, s, F_inviscid); 803 804 CeedScalar Flux[5] = {0.}; 805 for (int j=0; j<3; j++) { 806 Flux[0] += F_inviscid[j].density * norm[j]; 807 for (int k=0; k<3; k++) 808 Flux[k+1] += (F_inviscid[j].momentum[k] - stress[k][j]) * norm[j]; 809 Flux[4] += (F_inviscid[j].E_total + Fe[j])*norm[j]; 810 } 811 812 // -- Density 813 v[0][i] = -wdetJb * Flux[0]; 814 815 // -- Momentum 816 for (CeedInt j=0; j<3; j++) 817 v[j+1][i] = -wdetJb * Flux[j+1]; 818 819 // -- Total Energy Density 820 v[4][i] = -wdetJb * Flux[4]; 821 822 jac_data_sur[0][i] = s.U.density; 823 jac_data_sur[1][i] = s.Y.velocity[0]; 824 jac_data_sur[2][i] = s.Y.velocity[1]; 825 jac_data_sur[3][i] = s.Y.velocity[2]; 826 jac_data_sur[4][i] = s.U.E_total; 827 for (int j=0; j<6; j++) jac_data_sur[5+j][i] = kmstress[j]; 828 } 829 return 0; 830 } 831 832 // Jacobian for "set nothing" boundary integral 833 CEED_QFUNCTION(BoundaryIntegral_Jacobian)(void *ctx, CeedInt Q, 834 const CeedScalar *const *in, 835 CeedScalar *const *out) { 836 // *INDENT-OFF* 837 // Inputs 838 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 839 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 840 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 841 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 842 (*jac_data_sur)[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 // *INDENT-ON* 846 847 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 848 const bool implicit = context->is_implicit; 849 850 CeedPragmaSIMD 851 // Quadrature Point Loop 852 for (CeedInt i=0; i<Q; i++) { 853 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 854 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 855 const CeedScalar norm[3] = {q_data_sur[1][i], 856 q_data_sur[2][i], 857 q_data_sur[3][i] 858 }; 859 const CeedScalar dXdx[2][3] = { 860 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 861 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 862 }; 863 864 CeedScalar U[5], kmstress[6], dU[5], dx_i[3] = {0.}; 865 for (int j=0; j<5; j++) U[j] = jac_data_sur[j][i]; 866 for (int j=0; j<6; j++) kmstress[j] = jac_data_sur[5+j][i]; 867 for (int j=0; j<3; j++) U[j+1] *= U[0]; 868 for (int j=0; j<5; j++) dU[j] = dq[j][i]; 869 State s = StateFromU(context, U, x_i); 870 State ds = StateFromU_fwd(context, s, dU, x_i, dx_i); 871 872 State grad_ds[3]; 873 for (CeedInt j=0; j<3; j++) { 874 CeedScalar dx_i[3] = {0}, dUj[5]; 875 for (CeedInt k=0; k<5; k++) 876 dUj[k] = Grad_dq[0][k][i] * dXdx[0][j] + 877 Grad_dq[1][k][i] * dXdx[1][j]; 878 dx_i[j] = 1.; 879 grad_ds[j] = StateFromU_fwd(context, s, dUj, x_i, dx_i); 880 } 881 882 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 883 KMStrainRate(grad_ds, dstrain_rate); 884 NewtonianStress(context, dstrain_rate, dkmstress); 885 KMUnpack(dkmstress, dstress); 886 KMUnpack(kmstress, stress); 887 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 888 889 StateConservative dF_inviscid[3]; 890 FluxInviscid_fwd(context, s, ds, dF_inviscid); 891 892 CeedScalar dFlux[5] = {0.}; 893 for (int j=0; j<3; j++) { 894 dFlux[0] += dF_inviscid[j].density * norm[j]; 895 for (int k=0; k<3; k++) 896 dFlux[k+1] += (dF_inviscid[j].momentum[k] - dstress[k][j]) * norm[j]; 897 dFlux[4] += (dF_inviscid[j].E_total + dFe[j]) * norm[j]; 898 } 899 900 for (int j=0; j<5; j++) 901 v[j][i] = -wdetJb * dFlux[j]; 902 } // End Quadrature Point Loop 903 return 0; 904 } 905 906 // Outflow boundary condition, weakly setting a constant pressure 907 CEED_QFUNCTION(PressureOutflow)(void *ctx, CeedInt Q, 908 const CeedScalar *const *in, 909 CeedScalar *const *out) { 910 // *INDENT-OFF* 911 // Inputs 912 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 913 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 914 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 915 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 916 // Outputs 917 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 918 (*jac_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[1]; 919 // *INDENT-ON* 920 921 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 922 const bool implicit = context->is_implicit; 923 const CeedScalar P0 = context->P0; 924 925 CeedPragmaSIMD 926 // Quadrature Point Loop 927 for (CeedInt i=0; i<Q; i++) { 928 // Setup 929 // -- Interp in 930 const CeedScalar U[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 931 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 932 State s = StateFromU(context, U, x_i); 933 s.Y.pressure = P0; 934 935 // -- Interp-to-Interp q_data 936 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 937 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 938 // We can effect this by swapping the sign on this weight 939 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 940 941 // ---- Normal vect 942 const CeedScalar norm[3] = {q_data_sur[1][i], 943 q_data_sur[2][i], 944 q_data_sur[3][i] 945 }; 946 947 const CeedScalar dXdx[2][3] = { 948 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 949 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 950 }; 951 952 State grad_s[3]; 953 for (CeedInt j=0; j<3; j++) { 954 CeedScalar dx_i[3] = {0}, dU[5]; 955 for (CeedInt k=0; k<5; k++) 956 dU[k] = Grad_q[0][k][i] * dXdx[0][j] + 957 Grad_q[1][k][i] * dXdx[1][j]; 958 dx_i[j] = 1.; 959 grad_s[j] = StateFromU_fwd(context, s, dU, x_i, dx_i); 960 } 961 962 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 963 KMStrainRate(grad_s, strain_rate); 964 NewtonianStress(context, strain_rate, kmstress); 965 KMUnpack(kmstress, stress); 966 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 967 968 StateConservative F_inviscid[3]; 969 FluxInviscid(context, s, F_inviscid); 970 971 CeedScalar Flux[5] = {0.}; 972 for (int j=0; j<3; j++) { 973 Flux[0] += F_inviscid[j].density * norm[j]; 974 for (int k=0; k<3; k++) 975 Flux[k+1] += (F_inviscid[j].momentum[k] - stress[k][j]) * norm[j]; 976 Flux[4] += (F_inviscid[j].E_total + Fe[j])*norm[j]; 977 } 978 979 // -- Density 980 v[0][i] = -wdetJb * Flux[0]; 981 982 // -- Momentum 983 for (CeedInt j=0; j<3; j++) 984 v[j+1][i] = -wdetJb * Flux[j+1]; 985 986 // -- Total Energy Density 987 v[4][i] = -wdetJb * Flux[4]; 988 989 // Save values for Jacobian 990 jac_data_sur[0][i] = s.U.density; 991 jac_data_sur[1][i] = s.Y.velocity[0]; 992 jac_data_sur[2][i] = s.Y.velocity[1]; 993 jac_data_sur[3][i] = s.Y.velocity[2]; 994 jac_data_sur[4][i] = s.U.E_total; 995 for (int j=0; j<6; j++) jac_data_sur[5+j][i] = kmstress[j]; 996 } // End Quadrature Point Loop 997 return 0; 998 } 999 1000 // Jacobian for weak-pressure outflow boundary condition 1001 CEED_QFUNCTION(PressureOutflow_Jacobian)(void *ctx, CeedInt Q, 1002 const CeedScalar *const *in, 1003 CeedScalar *const *out) { 1004 // *INDENT-OFF* 1005 // Inputs 1006 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 1007 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 1008 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 1009 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 1010 (*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 1011 // Outputs 1012 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 1013 // *INDENT-ON* 1014 1015 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 1016 const bool implicit = context->is_implicit; 1017 1018 CeedPragmaSIMD 1019 // Quadrature Point Loop 1020 for (CeedInt i=0; i<Q; i++) { 1021 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 1022 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 1023 const CeedScalar norm[3] = {q_data_sur[1][i], 1024 q_data_sur[2][i], 1025 q_data_sur[3][i] 1026 }; 1027 const CeedScalar dXdx[2][3] = { 1028 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 1029 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 1030 }; 1031 1032 CeedScalar U[5], kmstress[6], dU[5], dx_i[3] = {0.}; 1033 for (int j=0; j<5; j++) U[j] = jac_data_sur[j][i]; 1034 for (int j=0; j<6; j++) kmstress[j] = jac_data_sur[5+j][i]; 1035 for (int j=0; j<3; j++) U[j+1] *= U[0]; 1036 for (int j=0; j<5; j++) dU[j] = dq[j][i]; 1037 State s = StateFromU(context, U, x_i); 1038 State ds = StateFromU_fwd(context, s, dU, x_i, dx_i); 1039 s.Y.pressure = context->P0; 1040 ds.Y.pressure = 0.; 1041 1042 State grad_ds[3]; 1043 for (CeedInt j=0; j<3; j++) { 1044 CeedScalar dx_i[3] = {0}, dUj[5]; 1045 for (CeedInt k=0; k<5; k++) 1046 dUj[k] = Grad_dq[0][k][i] * dXdx[0][j] + 1047 Grad_dq[1][k][i] * dXdx[1][j]; 1048 dx_i[j] = 1.; 1049 grad_ds[j] = StateFromU_fwd(context, s, dUj, x_i, dx_i); 1050 } 1051 1052 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 1053 KMStrainRate(grad_ds, dstrain_rate); 1054 NewtonianStress(context, dstrain_rate, dkmstress); 1055 KMUnpack(dkmstress, dstress); 1056 KMUnpack(kmstress, stress); 1057 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 1058 1059 StateConservative dF_inviscid[3]; 1060 FluxInviscid_fwd(context, s, ds, dF_inviscid); 1061 1062 CeedScalar dFlux[5] = {0.}; 1063 for (int j=0; j<3; j++) { 1064 dFlux[0] += dF_inviscid[j].density * norm[j]; 1065 for (int k=0; k<3; k++) 1066 dFlux[k+1] += (dF_inviscid[j].momentum[k] - dstress[k][j]) * norm[j]; 1067 dFlux[4] += (dF_inviscid[j].E_total + dFe[j]) * norm[j]; 1068 } 1069 1070 for (int j=0; j<5; j++) 1071 v[j][i] = -wdetJb * dFlux[j]; 1072 } // End Quadrature Point Loop 1073 return 0; 1074 } 1075 1076 // ***************************************************************************** 1077 #endif // newtonian_h 1078