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