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 #ifndef newtonian_h 12 #define newtonian_h 13 14 #include <ceed.h> 15 #include <math.h> 16 #include <stdlib.h> 17 18 #include "newtonian_state.h" 19 #include "newtonian_types.h" 20 #include "stabilization.h" 21 #include "utils.h" 22 23 // ***************************************************************************** 24 // This QFunction sets a "still" initial condition for generic Newtonian IG problems 25 // ***************************************************************************** 26 CEED_QFUNCTION(ICsNewtonianIG)(void *ctx, CeedInt Q, 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 for (CeedInt i = 0; i < Q; i++) { 44 CeedScalar q[5] = {0.}; 45 46 // Setup 47 // -- Coordinates 48 const CeedScalar x[3] = {X[0][i], X[1][i], X[2][i]}; 49 const CeedScalar e_potential = -Dot3(g, x); 50 51 // -- Density 52 const CeedScalar rho = P0 / (Rd * theta0); 53 54 // Initial Conditions 55 q[0] = rho; 56 q[1] = 0.0; 57 q[2] = 0.0; 58 q[3] = 0.0; 59 q[4] = rho * (cv * theta0 + e_potential); 60 61 for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j]; 62 63 } // End of Quadrature Point Loop 64 return 0; 65 } 66 67 // ***************************************************************************** 68 // This QFunction sets a "still" initial condition for generic Newtonian IG 69 // problems in primitive variables 70 // ***************************************************************************** 71 CEED_QFUNCTION(ICsNewtonianIG_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 72 // Outputs 73 CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 74 75 // Context 76 const SetupContext context = (SetupContext)ctx; 77 const CeedScalar theta0 = context->theta0; 78 const CeedScalar P0 = context->P0; 79 80 // Quadrature Point Loop 81 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 82 CeedScalar q[5] = {0.}; 83 84 // Initial Conditions 85 q[0] = P0; 86 q[1] = 0.0; 87 q[2] = 0.0; 88 q[3] = 0.0; 89 q[4] = theta0; 90 91 for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j]; 92 93 } // End of Quadrature Point Loop 94 return 0; 95 } 96 97 // ***************************************************************************** 98 // This QFunction implements the following formulation of Navier-Stokes with 99 // explicit time stepping method 100 // 101 // This is 3D compressible Navier-Stokes in conservation form with state 102 // variables of density, momentum density, and total energy density. 103 // 104 // State Variables: q = ( rho, U1, U2, U3, E ) 105 // rho - Mass Density 106 // Ui - Momentum Density, Ui = rho ui 107 // E - Total Energy Density, E = rho (cv T + (u u)/2 + g z) 108 // 109 // Navier-Stokes Equations: 110 // drho/dt + div( U ) = 0 111 // dU/dt + div( rho (u x u) + P I3 ) + rho g khat = div( Fu ) 112 // dE/dt + div( (E + P) u ) = div( Fe ) 113 // 114 // Viscous Stress: 115 // Fu = mu (grad( u ) + grad( u )^T + lambda div ( u ) I3) 116 // 117 // Thermal Stress: 118 // Fe = u Fu + k grad( T ) 119 // Equation of State 120 // P = (gamma - 1) (E - rho (u u) / 2 - rho g z) 121 // 122 // Stabilization: 123 // Tau = diag(TauC, TauM, TauM, TauM, TauE) 124 // f1 = rho sqrt(ui uj gij) 125 // gij = dXi/dX * dXi/dX 126 // TauC = Cc f1 / (8 gii) 127 // TauM = min( 1 , 1 / f1 ) 128 // TauE = TauM / (Ce cv) 129 // 130 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 131 // 132 // Constants: 133 // lambda = - 2 / 3, From Stokes hypothesis 134 // mu , Dynamic viscosity 135 // k , Thermal conductivity 136 // cv , Specific heat, constant volume 137 // cp , Specific heat, constant pressure 138 // g , Gravity 139 // gamma = cp / cv, Specific heat ratio 140 // 141 // We require the product of the inverse of the Jacobian (dXdx_j,k) and 142 // its transpose (dXdx_k,j) to properly compute integrals of the form: 143 // int( gradv gradu ) 144 // 145 // ***************************************************************************** 146 CEED_QFUNCTION(RHSFunction_Newtonian)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 147 // Inputs 148 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 149 const CeedScalar(*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 150 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 151 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 152 153 // Outputs 154 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 155 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 156 157 // Context 158 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 159 const CeedScalar *g = context->g; 160 const CeedScalar dt = context->dt; 161 162 // Quadrature Point Loop 163 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 164 CeedScalar U[5]; 165 for (int j = 0; j < 5; j++) U[j] = q[j][i]; 166 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 167 State s = StateFromU(context, U, x_i); 168 169 // -- Interp-to-Interp q_data 170 const CeedScalar wdetJ = q_data[0][i]; 171 // -- Interp-to-Grad q_data 172 // ---- Inverse of change of coordinate matrix: X_i,j 173 const CeedScalar dXdx[3][3] = { 174 {q_data[1][i], q_data[2][i], q_data[3][i]}, 175 {q_data[4][i], q_data[5][i], q_data[6][i]}, 176 {q_data[7][i], q_data[8][i], q_data[9][i]} 177 }; 178 State grad_s[3]; 179 for (CeedInt j = 0; j < 3; j++) { 180 CeedScalar dx_i[3] = {0}, dU[5]; 181 for (CeedInt k = 0; k < 5; k++) dU[k] = Grad_q[0][k][i] * dXdx[0][j] + Grad_q[1][k][i] * dXdx[1][j] + Grad_q[2][k][i] * dXdx[2][j]; 182 dx_i[j] = 1.; 183 grad_s[j] = StateFromU_fwd(context, s, dU, x_i, dx_i); 184 } 185 186 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 187 KMStrainRate(grad_s, strain_rate); 188 NewtonianStress(context, strain_rate, kmstress); 189 KMUnpack(kmstress, stress); 190 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 191 192 StateConservative F_inviscid[3]; 193 FluxInviscid(context, s, F_inviscid); 194 195 // Total flux 196 CeedScalar Flux[5][3]; 197 FluxTotal(F_inviscid, stress, Fe, Flux); 198 199 for (CeedInt j = 0; j < 3; j++) { 200 for (CeedInt k = 0; k < 5; k++) Grad_v[j][k][i] = wdetJ * (dXdx[j][0] * Flux[k][0] + dXdx[j][1] * Flux[k][1] + dXdx[j][2] * Flux[k][2]); 201 } 202 203 const CeedScalar body_force[5] = {0, s.U.density * g[0], s.U.density * g[1], s.U.density * g[2], 0}; 204 for (int j = 0; j < 5; j++) v[j][i] = wdetJ * body_force[j]; 205 206 // -- Stabilization method: none (Galerkin), SU, or SUPG 207 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}; 208 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 209 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, x_i, stab); 210 211 for (CeedInt j = 0; j < 5; j++) { 212 for (CeedInt k = 0; k < 3; k++) Grad_v[k][j][i] -= wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXdx[k][1] + stab[j][2] * dXdx[k][2]); 213 } 214 } // End Quadrature Point Loop 215 216 // Return 217 return 0; 218 } 219 220 // ***************************************************************************** 221 // This QFunction implements the Navier-Stokes equations (mentioned above) with 222 // implicit time stepping method 223 // 224 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 225 // SUPG = Galerkin + grad(v) . ( Ai^T * Tau * (q_dot + Aj q,j - body force) ) 226 // (diffussive terms will be added later) 227 // 228 // ***************************************************************************** 229 CEED_QFUNCTION_HELPER int IFunction_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateFromQi_t StateFromQi, 230 StateFromQi_fwd_t StateFromQi_fwd) { 231 // Inputs 232 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 233 const CeedScalar(*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 234 const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 235 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 236 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 237 238 // Outputs 239 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 240 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 241 CeedScalar(*jac_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[2]; 242 243 // Context 244 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 245 const CeedScalar *g = context->g; 246 const CeedScalar dt = context->dt; 247 248 // Quadrature Point Loop 249 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 250 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 251 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 252 const State s = StateFromQi(context, qi, x_i); 253 254 // -- Interp-to-Interp q_data 255 const CeedScalar wdetJ = q_data[0][i]; 256 // -- Interp-to-Grad q_data 257 // ---- Inverse of change of coordinate matrix: X_i,j 258 const CeedScalar dXdx[3][3] = { 259 {q_data[1][i], q_data[2][i], q_data[3][i]}, 260 {q_data[4][i], q_data[5][i], q_data[6][i]}, 261 {q_data[7][i], q_data[8][i], q_data[9][i]} 262 }; 263 State grad_s[3]; 264 for (CeedInt j = 0; j < 3; j++) { 265 CeedScalar dx_i[3] = {0}, dqi[5]; 266 for (CeedInt k = 0; k < 5; k++) { 267 dqi[k] = Grad_q[0][k][i] * dXdx[0][j] + Grad_q[1][k][i] * dXdx[1][j] + Grad_q[2][k][i] * dXdx[2][j]; 268 } 269 dx_i[j] = 1.; 270 grad_s[j] = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 271 } 272 273 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 274 KMStrainRate(grad_s, strain_rate); 275 NewtonianStress(context, strain_rate, kmstress); 276 KMUnpack(kmstress, stress); 277 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 278 279 StateConservative F_inviscid[3]; 280 FluxInviscid(context, s, F_inviscid); 281 282 // Total flux 283 CeedScalar Flux[5][3]; 284 FluxTotal(F_inviscid, stress, Fe, Flux); 285 286 for (CeedInt j = 0; j < 3; j++) { 287 for (CeedInt k = 0; k < 5; k++) { 288 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * Flux[k][0] + dXdx[j][1] * Flux[k][1] + dXdx[j][2] * Flux[k][2]); 289 } 290 } 291 292 const CeedScalar body_force[5] = {0, s.U.density * g[0], s.U.density * g[1], s.U.density * g[2], 0}; 293 294 // -- Stabilization method: none (Galerkin), SU, or SUPG 295 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}, qi_dot[5], dx0[3] = {0}; 296 for (int j = 0; j < 5; j++) qi_dot[j] = q_dot[j][i]; 297 State s_dot = StateFromQi_fwd(context, s, qi_dot, x_i, dx0); 298 UnpackState_U(s_dot.U, U_dot); 299 300 for (CeedInt j = 0; j < 5; j++) v[j][i] = wdetJ * (U_dot[j] - body_force[j]); 301 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 302 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, x_i, stab); 303 304 for (CeedInt j = 0; j < 5; j++) { 305 for (CeedInt k = 0; k < 3; k++) { 306 Grad_v[k][j][i] += wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXdx[k][1] + stab[j][2] * dXdx[k][2]); 307 } 308 } 309 for (CeedInt j = 0; j < 5; j++) jac_data[j][i] = qi[j]; 310 for (CeedInt j = 0; j < 6; j++) jac_data[5 + j][i] = kmstress[j]; 311 for (CeedInt j = 0; j < 3; j++) jac_data[5 + 6 + j][i] = Tau_d[j]; 312 313 } // End Quadrature Point Loop 314 315 // Return 316 return 0; 317 } 318 319 CEED_QFUNCTION(IFunction_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 320 return IFunction_Newtonian(ctx, Q, in, out, StateFromU, StateFromU_fwd); 321 } 322 323 CEED_QFUNCTION(IFunction_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 324 return IFunction_Newtonian(ctx, Q, in, out, StateFromY, StateFromY_fwd); 325 } 326 327 // ***************************************************************************** 328 // This QFunction implements the jacobian of the Navier-Stokes equations 329 // for implicit time stepping method. 330 // ***************************************************************************** 331 CEED_QFUNCTION_HELPER int IJacobian_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateFromQi_t StateFromQi, 332 StateFromQi_fwd_t StateFromQi_fwd) { 333 // Inputs 334 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 335 const CeedScalar(*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 336 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 337 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 338 const CeedScalar(*jac_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 339 340 // Outputs 341 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 342 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 343 344 // Context 345 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 346 const CeedScalar *g = context->g; 347 348 // Quadrature Point Loop 349 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 350 // -- Interp-to-Interp q_data 351 const CeedScalar wdetJ = q_data[0][i]; 352 // -- Interp-to-Grad q_data 353 // ---- Inverse of change of coordinate matrix: X_i,j 354 const CeedScalar dXdx[3][3] = { 355 {q_data[1][i], q_data[2][i], q_data[3][i]}, 356 {q_data[4][i], q_data[5][i], q_data[6][i]}, 357 {q_data[7][i], q_data[8][i], q_data[9][i]} 358 }; 359 360 CeedScalar qi[5], kmstress[6], Tau_d[3]; 361 for (int j = 0; j < 5; j++) qi[j] = jac_data[j][i]; 362 for (int j = 0; j < 6; j++) kmstress[j] = jac_data[5 + j][i]; 363 for (int j = 0; j < 3; j++) Tau_d[j] = jac_data[5 + 6 + j][i]; 364 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 365 State s = StateFromQi(context, qi, x_i); 366 367 CeedScalar dqi[5], dx0[3] = {0}; 368 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 369 State ds = StateFromQi_fwd(context, s, dqi, x_i, dx0); 370 371 State grad_ds[3]; 372 for (int j = 0; j < 3; j++) { 373 CeedScalar dqi_j[5]; 374 for (int k = 0; k < 5; k++) dqi_j[k] = Grad_dq[0][k][i] * dXdx[0][j] + Grad_dq[1][k][i] * dXdx[1][j] + Grad_dq[2][k][i] * dXdx[2][j]; 375 grad_ds[j] = StateFromQi_fwd(context, s, dqi_j, x_i, dx0); 376 } 377 378 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 379 KMStrainRate(grad_ds, dstrain_rate); 380 NewtonianStress(context, dstrain_rate, dkmstress); 381 KMUnpack(dkmstress, dstress); 382 KMUnpack(kmstress, stress); 383 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 384 385 StateConservative dF_inviscid[3]; 386 FluxInviscid_fwd(context, s, ds, dF_inviscid); 387 388 // Total flux 389 CeedScalar dFlux[5][3]; 390 FluxTotal(dF_inviscid, dstress, dFe, dFlux); 391 392 for (int j = 0; j < 3; j++) { 393 for (int k = 0; k < 5; k++) Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * dFlux[k][0] + dXdx[j][1] * dFlux[k][1] + dXdx[j][2] * dFlux[k][2]); 394 } 395 396 const CeedScalar dbody_force[5] = {0, ds.U.density * g[0], ds.U.density * g[1], ds.U.density * g[2], 0}; 397 CeedScalar dU[5] = {0.}; 398 UnpackState_U(ds.U, dU); 399 for (int j = 0; j < 5; j++) v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]); 400 401 // -- Stabilization method: none (Galerkin), SU, or SUPG 402 CeedScalar dstab[5][3], U_dot[5] = {0}; 403 for (CeedInt j = 0; j < 5; j++) U_dot[j] = context->ijacobian_time_shift * dU[j]; 404 Stabilization(context, s, Tau_d, grad_ds, U_dot, dbody_force, x_i, dstab); 405 406 for (int j = 0; j < 5; j++) { 407 for (int k = 0; k < 3; k++) Grad_v[k][j][i] += wdetJ * (dstab[j][0] * dXdx[k][0] + dstab[j][1] * dXdx[k][1] + dstab[j][2] * dXdx[k][2]); 408 } 409 } // End Quadrature Point Loop 410 return 0; 411 } 412 413 CEED_QFUNCTION(IJacobian_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 414 return IJacobian_Newtonian(ctx, Q, in, out, StateFromU, StateFromU_fwd); 415 } 416 417 CEED_QFUNCTION(IJacobian_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 418 return IJacobian_Newtonian(ctx, Q, in, out, StateFromY, StateFromY_fwd); 419 } 420 421 // ***************************************************************************** 422 // Compute boundary integral (ie. for strongly set inflows) 423 // ***************************************************************************** 424 CEED_QFUNCTION_HELPER int BoundaryIntegral(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateFromQi_t StateFromQi, 425 StateFromQi_fwd_t StateFromQi_fwd) { 426 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 427 const CeedScalar(*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 428 const CeedScalar(*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 429 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 430 431 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 432 CeedScalar(*jac_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[1]; 433 434 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 435 const bool is_implicit = context->is_implicit; 436 437 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 438 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 439 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 440 State s = StateFromQi(context, qi, x_i); 441 442 const CeedScalar wdetJb = (is_implicit ? -1. : 1.) * q_data_sur[0][i]; 443 // ---- Normal vector 444 const CeedScalar norm[3] = {q_data_sur[1][i], q_data_sur[2][i], q_data_sur[3][i]}; 445 446 const CeedScalar dXdx[2][3] = { 447 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 448 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 449 }; 450 451 State grad_s[3]; 452 for (CeedInt j = 0; j < 3; j++) { 453 CeedScalar dx_i[3] = {0}, dqi[5]; 454 for (CeedInt k = 0; k < 5; k++) dqi[k] = Grad_q[0][k][i] * dXdx[0][j] + Grad_q[1][k][i] * dXdx[1][j]; 455 dx_i[j] = 1.; 456 grad_s[j] = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 457 } 458 459 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 460 KMStrainRate(grad_s, strain_rate); 461 NewtonianStress(context, strain_rate, kmstress); 462 KMUnpack(kmstress, stress); 463 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 464 465 StateConservative F_inviscid[3]; 466 FluxInviscid(context, s, F_inviscid); 467 468 CeedScalar Flux[5]; 469 FluxTotal_Boundary(F_inviscid, stress, Fe, norm, Flux); 470 471 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j]; 472 473 for (int j = 0; j < 5; j++) jac_data_sur[j][i] = qi[j]; 474 for (int j = 0; j < 6; j++) jac_data_sur[5 + j][i] = kmstress[j]; 475 } 476 return 0; 477 } 478 479 CEED_QFUNCTION(BoundaryIntegral_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 480 return BoundaryIntegral(ctx, Q, in, out, StateFromU, StateFromU_fwd); 481 } 482 483 CEED_QFUNCTION(BoundaryIntegral_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 484 return BoundaryIntegral(ctx, Q, in, out, StateFromY, StateFromY_fwd); 485 } 486 487 // ***************************************************************************** 488 // Jacobian for "set nothing" boundary integral 489 // ***************************************************************************** 490 CEED_QFUNCTION_HELPER int BoundaryIntegral_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 491 StateFromQi_t StateFromQi, StateFromQi_fwd_t StateFromQi_fwd) { 492 // Inputs 493 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 494 const CeedScalar(*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 495 const CeedScalar(*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 496 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 497 const CeedScalar(*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 498 499 // Outputs 500 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 501 502 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 503 const bool implicit = context->is_implicit; 504 505 // Quadrature Point Loop 506 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 507 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 508 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 509 const CeedScalar norm[3] = {q_data_sur[1][i], q_data_sur[2][i], q_data_sur[3][i]}; 510 const CeedScalar dXdx[2][3] = { 511 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 512 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 513 }; 514 515 CeedScalar qi[5], kmstress[6], dqi[5], dx_i[3] = {0.}; 516 for (int j = 0; j < 5; j++) qi[j] = jac_data_sur[j][i]; 517 for (int j = 0; j < 6; j++) kmstress[j] = jac_data_sur[5 + j][i]; 518 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 519 520 State s = StateFromQi(context, qi, x_i); 521 State ds = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 522 523 State grad_ds[3]; 524 for (CeedInt j = 0; j < 3; j++) { 525 CeedScalar dx_i[3] = {0}, dqi_j[5]; 526 for (CeedInt k = 0; k < 5; k++) dqi_j[k] = Grad_dq[0][k][i] * dXdx[0][j] + Grad_dq[1][k][i] * dXdx[1][j]; 527 dx_i[j] = 1.; 528 grad_ds[j] = StateFromQi_fwd(context, s, dqi_j, x_i, dx_i); 529 } 530 531 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 532 KMStrainRate(grad_ds, dstrain_rate); 533 NewtonianStress(context, dstrain_rate, dkmstress); 534 KMUnpack(dkmstress, dstress); 535 KMUnpack(kmstress, stress); 536 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 537 538 StateConservative dF_inviscid[3]; 539 FluxInviscid_fwd(context, s, ds, dF_inviscid); 540 541 CeedScalar dFlux[5]; 542 FluxTotal_Boundary(dF_inviscid, dstress, dFe, norm, dFlux); 543 544 for (int j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j]; 545 } // End Quadrature Point Loop 546 return 0; 547 } 548 549 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 550 return BoundaryIntegral_Jacobian(ctx, Q, in, out, StateFromU, StateFromU_fwd); 551 } 552 553 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 554 return BoundaryIntegral_Jacobian(ctx, Q, in, out, StateFromY, StateFromY_fwd); 555 } 556 557 // ***************************************************************************** 558 // Outflow boundary condition, weakly setting a constant pressure 559 // ***************************************************************************** 560 CEED_QFUNCTION_HELPER int PressureOutflow(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateFromQi_t StateFromQi, 561 StateFromQi_fwd_t StateFromQi_fwd) { 562 // Inputs 563 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 564 const CeedScalar(*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 565 const CeedScalar(*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 566 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 567 568 // Outputs 569 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 570 CeedScalar(*jac_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[1]; 571 572 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 573 const bool implicit = context->is_implicit; 574 const CeedScalar P0 = context->P0; 575 576 // Quadrature Point Loop 577 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 578 // Setup 579 // -- Interp in 580 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 581 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 582 State s = StateFromQi(context, qi, x_i); 583 s.Y.pressure = P0; 584 585 // -- Interp-to-Interp q_data 586 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 587 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 588 // We can effect this by swapping the sign on this weight 589 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 590 591 // ---- Normal vector 592 const CeedScalar norm[3] = {q_data_sur[1][i], q_data_sur[2][i], q_data_sur[3][i]}; 593 594 const CeedScalar dXdx[2][3] = { 595 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 596 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 597 }; 598 599 State grad_s[3]; 600 for (CeedInt j = 0; j < 3; j++) { 601 CeedScalar dx_i[3] = {0}, dqi[5]; 602 for (CeedInt k = 0; k < 5; k++) dqi[k] = Grad_q[0][k][i] * dXdx[0][j] + Grad_q[1][k][i] * dXdx[1][j]; 603 dx_i[j] = 1.; 604 grad_s[j] = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 605 } 606 607 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 608 KMStrainRate(grad_s, strain_rate); 609 NewtonianStress(context, strain_rate, kmstress); 610 KMUnpack(kmstress, stress); 611 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 612 613 StateConservative F_inviscid[3]; 614 FluxInviscid(context, s, F_inviscid); 615 616 CeedScalar Flux[5]; 617 FluxTotal_Boundary(F_inviscid, stress, Fe, norm, Flux); 618 619 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j]; 620 621 // Save values for Jacobian 622 for (int j = 0; j < 5; j++) jac_data_sur[j][i] = qi[j]; 623 for (int j = 0; j < 6; j++) jac_data_sur[5 + j][i] = kmstress[j]; 624 } // End Quadrature Point Loop 625 return 0; 626 } 627 628 CEED_QFUNCTION(PressureOutflow_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 629 return PressureOutflow(ctx, Q, in, out, StateFromU, StateFromU_fwd); 630 } 631 632 CEED_QFUNCTION(PressureOutflow_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 633 return PressureOutflow(ctx, Q, in, out, StateFromY, StateFromY_fwd); 634 } 635 636 // ***************************************************************************** 637 // Jacobian for weak-pressure outflow boundary condition 638 // ***************************************************************************** 639 CEED_QFUNCTION_HELPER int PressureOutflow_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 640 StateFromQi_t StateFromQi, StateFromQi_fwd_t StateFromQi_fwd) { 641 // Inputs 642 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 643 const CeedScalar(*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 644 const CeedScalar(*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 645 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 646 const CeedScalar(*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 647 648 // Outputs 649 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 650 651 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 652 const bool implicit = context->is_implicit; 653 654 // Quadrature Point Loop 655 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 656 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 657 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 658 const CeedScalar norm[3] = {q_data_sur[1][i], q_data_sur[2][i], q_data_sur[3][i]}; 659 const CeedScalar dXdx[2][3] = { 660 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 661 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 662 }; 663 664 CeedScalar qi[5], kmstress[6], dqi[5], dx_i[3] = {0.}; 665 for (int j = 0; j < 5; j++) qi[j] = jac_data_sur[j][i]; 666 for (int j = 0; j < 6; j++) kmstress[j] = jac_data_sur[5 + j][i]; 667 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 668 669 State s = StateFromQi(context, qi, x_i); 670 State ds = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 671 s.Y.pressure = context->P0; 672 ds.Y.pressure = 0.; 673 674 State grad_ds[3]; 675 for (CeedInt j = 0; j < 3; j++) { 676 CeedScalar dx_i[3] = {0}, dqi_j[5]; 677 for (CeedInt k = 0; k < 5; k++) dqi_j[k] = Grad_dq[0][k][i] * dXdx[0][j] + Grad_dq[1][k][i] * dXdx[1][j]; 678 dx_i[j] = 1.; 679 grad_ds[j] = StateFromQi_fwd(context, s, dqi_j, x_i, dx_i); 680 } 681 682 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 683 KMStrainRate(grad_ds, dstrain_rate); 684 NewtonianStress(context, dstrain_rate, dkmstress); 685 KMUnpack(dkmstress, dstress); 686 KMUnpack(kmstress, stress); 687 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 688 689 StateConservative dF_inviscid[3]; 690 FluxInviscid_fwd(context, s, ds, dF_inviscid); 691 692 CeedScalar dFlux[5]; 693 FluxTotal_Boundary(dF_inviscid, dstress, dFe, norm, dFlux); 694 695 for (int j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j]; 696 } // End Quadrature Point Loop 697 return 0; 698 } 699 700 CEED_QFUNCTION(PressureOutflow_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 701 return PressureOutflow_Jacobian(ctx, Q, in, out, StateFromU, StateFromU_fwd); 702 } 703 704 CEED_QFUNCTION(PressureOutflow_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 705 return PressureOutflow_Jacobian(ctx, Q, in, out, StateFromY, StateFromY_fwd); 706 } 707 708 #endif // newtonian_h 709