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