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], (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 149 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 150 // Outputs 151 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 152 153 // Context 154 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 155 const CeedScalar *g = context->g; 156 const CeedScalar dt = context->dt; 157 158 CeedPragmaSIMD 159 // Quadrature Point Loop 160 for (CeedInt i = 0; i < Q; i++) { 161 CeedScalar U[5]; 162 for (int j = 0; j < 5; j++) U[j] = q[j][i]; 163 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 164 State s = StateFromU(context, U, x_i); 165 166 // -- Interp-to-Interp q_data 167 const CeedScalar wdetJ = q_data[0][i]; 168 // -- Interp-to-Grad q_data 169 // ---- Inverse of change of coordinate matrix: X_i,j 170 const CeedScalar dXdx[3][3] = { 171 {q_data[1][i], q_data[2][i], q_data[3][i]}, 172 {q_data[4][i], q_data[5][i], q_data[6][i]}, 173 {q_data[7][i], q_data[8][i], q_data[9][i]} 174 }; 175 State grad_s[3]; 176 for (CeedInt j = 0; j < 3; j++) { 177 CeedScalar dx_i[3] = {0}, dU[5]; 178 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]; 179 dx_i[j] = 1.; 180 grad_s[j] = StateFromU_fwd(context, s, dU, x_i, dx_i); 181 } 182 183 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 184 KMStrainRate(grad_s, strain_rate); 185 NewtonianStress(context, strain_rate, kmstress); 186 KMUnpack(kmstress, stress); 187 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 188 189 StateConservative F_inviscid[3]; 190 FluxInviscid(context, s, F_inviscid); 191 192 // Total flux 193 CeedScalar Flux[5][3]; 194 FluxTotal(F_inviscid, stress, Fe, Flux); 195 196 for (CeedInt j = 0; j < 3; j++) { 197 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]); 198 } 199 200 const CeedScalar body_force[5] = {0, s.U.density * g[0], s.U.density * g[1], s.U.density * g[2], 0}; 201 for (int j = 0; j < 5; j++) v[j][i] = wdetJ * body_force[j]; 202 203 // -- Stabilization method: none (Galerkin), SU, or SUPG 204 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}; 205 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 206 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, x_i, stab); 207 208 for (CeedInt j = 0; j < 5; j++) { 209 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]); 210 } 211 } // End Quadrature Point Loop 212 213 // Return 214 return 0; 215 } 216 217 // ***************************************************************************** 218 // This QFunction implements the Navier-Stokes equations (mentioned above) with 219 // implicit time stepping method 220 // 221 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 222 // SUPG = Galerkin + grad(v) . ( Ai^T * Tau * (q_dot + Aj q,j - body force) ) 223 // (diffussive terms will be added later) 224 // 225 // ***************************************************************************** 226 CEED_QFUNCTION_HELPER int IFunction_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateFromQi_t StateFromQi, 227 StateFromQi_fwd_t StateFromQi_fwd) { 228 // Inputs 229 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 230 (*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 231 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 232 // Outputs 233 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1], 234 (*jac_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[2]; 235 // Context 236 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 237 const CeedScalar *g = context->g; 238 const CeedScalar dt = context->dt; 239 240 CeedPragmaSIMD 241 // Quadrature Point Loop 242 for (CeedInt i = 0; i < Q; i++) { 243 CeedScalar qi[5]; 244 for (CeedInt j = 0; j < 5; j++) qi[j] = q[j][i]; 245 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 246 State s = StateFromQi(context, qi, x_i); 247 248 // -- Interp-to-Interp q_data 249 const CeedScalar wdetJ = q_data[0][i]; 250 // -- Interp-to-Grad q_data 251 // ---- Inverse of change of coordinate matrix: X_i,j 252 const CeedScalar dXdx[3][3] = { 253 {q_data[1][i], q_data[2][i], q_data[3][i]}, 254 {q_data[4][i], q_data[5][i], q_data[6][i]}, 255 {q_data[7][i], q_data[8][i], q_data[9][i]} 256 }; 257 State grad_s[3]; 258 for (CeedInt j = 0; j < 3; j++) { 259 CeedScalar dx_i[3] = {0}, dqi[5]; 260 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] + Grad_q[2][k][i] * dXdx[2][j]; 261 dx_i[j] = 1.; 262 grad_s[j] = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 263 } 264 265 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 266 KMStrainRate(grad_s, strain_rate); 267 NewtonianStress(context, strain_rate, kmstress); 268 KMUnpack(kmstress, stress); 269 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 270 271 StateConservative F_inviscid[3]; 272 FluxInviscid(context, s, F_inviscid); 273 274 // Total flux 275 CeedScalar Flux[5][3]; 276 FluxTotal(F_inviscid, stress, Fe, Flux); 277 278 for (CeedInt j = 0; j < 3; j++) { 279 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]); 280 } 281 282 const CeedScalar body_force[5] = {0, s.U.density * g[0], s.U.density * g[1], s.U.density * g[2], 0}; 283 284 // -- Stabilization method: none (Galerkin), SU, or SUPG 285 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}, qi_dot[5], dx0[3] = {0}; 286 for (int j = 0; j < 5; j++) qi_dot[j] = q_dot[j][i]; 287 State s_dot = StateFromQi_fwd(context, s, qi_dot, x_i, dx0); 288 UnpackState_U(s_dot.U, U_dot); 289 290 for (CeedInt j = 0; j < 5; j++) v[j][i] = wdetJ * (U_dot[j] - body_force[j]); 291 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 292 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, x_i, stab); 293 294 for (CeedInt j = 0; j < 5; j++) { 295 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]); 296 } 297 for (CeedInt j = 0; j < 5; j++) jac_data[j][i] = qi[j]; 298 for (CeedInt j = 0; j < 6; j++) jac_data[5 + j][i] = kmstress[j]; 299 for (CeedInt j = 0; j < 3; j++) jac_data[5 + 6 + j][i] = Tau_d[j]; 300 301 } // End Quadrature Point Loop 302 303 // Return 304 return 0; 305 } 306 307 CEED_QFUNCTION(IFunction_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 308 return IFunction_Newtonian(ctx, Q, in, out, StateFromU, StateFromU_fwd); 309 } 310 311 CEED_QFUNCTION(IFunction_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 312 return IFunction_Newtonian(ctx, Q, in, out, StateFromY, StateFromY_fwd); 313 } 314 315 // ***************************************************************************** 316 // This QFunction implements the jacobian of the Navier-Stokes equations 317 // for implicit time stepping method. 318 // ***************************************************************************** 319 CEED_QFUNCTION_HELPER int IJacobian_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateFromQi_t StateFromQi, 320 StateFromQi_fwd_t StateFromQi_fwd) { 321 // Inputs 322 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 323 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 324 (*jac_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 325 // Outputs 326 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 327 // Context 328 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 329 const CeedScalar *g = context->g; 330 331 CeedPragmaSIMD 332 // Quadrature Point Loop 333 for (CeedInt i = 0; i < Q; i++) { 334 // -- Interp-to-Interp q_data 335 const CeedScalar wdetJ = q_data[0][i]; 336 // -- Interp-to-Grad q_data 337 // ---- Inverse of change of coordinate matrix: X_i,j 338 const CeedScalar dXdx[3][3] = { 339 {q_data[1][i], q_data[2][i], q_data[3][i]}, 340 {q_data[4][i], q_data[5][i], q_data[6][i]}, 341 {q_data[7][i], q_data[8][i], q_data[9][i]} 342 }; 343 344 CeedScalar qi[5], kmstress[6], Tau_d[3]; 345 for (int j = 0; j < 5; j++) qi[j] = jac_data[j][i]; 346 for (int j = 0; j < 6; j++) kmstress[j] = jac_data[5 + j][i]; 347 for (int j = 0; j < 3; j++) Tau_d[j] = jac_data[5 + 6 + j][i]; 348 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 349 State s = StateFromQi(context, qi, x_i); 350 351 CeedScalar dqi[5], dx0[3] = {0}; 352 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 353 State ds = StateFromQi_fwd(context, s, dqi, x_i, dx0); 354 355 State grad_ds[3]; 356 for (int j = 0; j < 3; j++) { 357 CeedScalar dqi_j[5]; 358 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]; 359 grad_ds[j] = StateFromQi_fwd(context, s, dqi_j, x_i, dx0); 360 } 361 362 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 363 KMStrainRate(grad_ds, dstrain_rate); 364 NewtonianStress(context, dstrain_rate, dkmstress); 365 KMUnpack(dkmstress, dstress); 366 KMUnpack(kmstress, stress); 367 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 368 369 StateConservative dF_inviscid[3]; 370 FluxInviscid_fwd(context, s, ds, dF_inviscid); 371 372 // Total flux 373 CeedScalar dFlux[5][3]; 374 FluxTotal(dF_inviscid, dstress, dFe, dFlux); 375 376 for (int j = 0; j < 3; j++) { 377 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]); 378 } 379 380 const CeedScalar dbody_force[5] = {0, ds.U.density * g[0], ds.U.density * g[1], ds.U.density * g[2], 0}; 381 CeedScalar dU[5] = {0.}; 382 UnpackState_U(ds.U, dU); 383 for (int j = 0; j < 5; j++) v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]); 384 385 // -- Stabilization method: none (Galerkin), SU, or SUPG 386 CeedScalar dstab[5][3], U_dot[5] = {0}; 387 for (CeedInt j = 0; j < 5; j++) U_dot[j] = context->ijacobian_time_shift * dU[j]; 388 Stabilization(context, s, Tau_d, grad_ds, U_dot, dbody_force, x_i, dstab); 389 390 for (int j = 0; j < 5; j++) { 391 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]); 392 } 393 } // End Quadrature Point Loop 394 return 0; 395 } 396 397 CEED_QFUNCTION(IJacobian_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 398 return IJacobian_Newtonian(ctx, Q, in, out, StateFromU, StateFromU_fwd); 399 } 400 401 CEED_QFUNCTION(IJacobian_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 402 return IJacobian_Newtonian(ctx, Q, in, out, StateFromY, StateFromY_fwd); 403 } 404 405 // ***************************************************************************** 406 // Compute boundary integral (ie. for strongly set inflows) 407 // ***************************************************************************** 408 CEED_QFUNCTION_HELPER int BoundaryIntegral(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateFromQi_t StateFromQi, 409 StateFromQi_fwd_t StateFromQi_fwd) { 410 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 411 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 412 413 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], (*jac_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[1]; 414 415 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 416 const bool is_implicit = context->is_implicit; 417 418 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 419 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 420 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 421 State s = StateFromQi(context, qi, x_i); 422 423 const CeedScalar wdetJb = (is_implicit ? -1. : 1.) * q_data_sur[0][i]; 424 // ---- Normal vector 425 const CeedScalar norm[3] = {q_data_sur[1][i], q_data_sur[2][i], q_data_sur[3][i]}; 426 427 const CeedScalar dXdx[2][3] = { 428 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 429 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 430 }; 431 432 State grad_s[3]; 433 for (CeedInt j = 0; j < 3; j++) { 434 CeedScalar dx_i[3] = {0}, dqi[5]; 435 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]; 436 dx_i[j] = 1.; 437 grad_s[j] = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 438 } 439 440 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 441 KMStrainRate(grad_s, strain_rate); 442 NewtonianStress(context, strain_rate, kmstress); 443 KMUnpack(kmstress, stress); 444 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 445 446 StateConservative F_inviscid[3]; 447 FluxInviscid(context, s, F_inviscid); 448 449 CeedScalar Flux[5]; 450 FluxTotal_Boundary(F_inviscid, stress, Fe, norm, Flux); 451 452 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j]; 453 454 for (int j = 0; j < 5; j++) jac_data_sur[j][i] = qi[j]; 455 for (int j = 0; j < 6; j++) jac_data_sur[5 + j][i] = kmstress[j]; 456 } 457 return 0; 458 } 459 460 CEED_QFUNCTION(BoundaryIntegral_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 461 return BoundaryIntegral(ctx, Q, in, out, StateFromU, StateFromU_fwd); 462 } 463 464 CEED_QFUNCTION(BoundaryIntegral_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 465 return BoundaryIntegral(ctx, Q, in, out, StateFromY, StateFromY_fwd); 466 } 467 468 // ***************************************************************************** 469 // Jacobian for "set nothing" boundary integral 470 // ***************************************************************************** 471 CEED_QFUNCTION_HELPER int BoundaryIntegral_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 472 StateFromQi_t StateFromQi, StateFromQi_fwd_t StateFromQi_fwd) { 473 // Inputs 474 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 475 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 476 (*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 477 // Outputs 478 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 479 480 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 481 const bool implicit = context->is_implicit; 482 483 CeedPragmaSIMD 484 // Quadrature Point Loop 485 for (CeedInt i = 0; i < Q; i++) { 486 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 487 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 488 const CeedScalar norm[3] = {q_data_sur[1][i], q_data_sur[2][i], q_data_sur[3][i]}; 489 const CeedScalar dXdx[2][3] = { 490 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 491 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 492 }; 493 494 CeedScalar qi[5], kmstress[6], dqi[5], dx_i[3] = {0.}; 495 for (int j = 0; j < 5; j++) qi[j] = jac_data_sur[j][i]; 496 for (int j = 0; j < 6; j++) kmstress[j] = jac_data_sur[5 + j][i]; 497 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 498 499 State s = StateFromQi(context, qi, x_i); 500 State ds = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 501 502 State grad_ds[3]; 503 for (CeedInt j = 0; j < 3; j++) { 504 CeedScalar dx_i[3] = {0}, dqi_j[5]; 505 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]; 506 dx_i[j] = 1.; 507 grad_ds[j] = StateFromQi_fwd(context, s, dqi_j, x_i, dx_i); 508 } 509 510 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 511 KMStrainRate(grad_ds, dstrain_rate); 512 NewtonianStress(context, dstrain_rate, dkmstress); 513 KMUnpack(dkmstress, dstress); 514 KMUnpack(kmstress, stress); 515 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 516 517 StateConservative dF_inviscid[3]; 518 FluxInviscid_fwd(context, s, ds, dF_inviscid); 519 520 CeedScalar dFlux[5]; 521 FluxTotal_Boundary(dF_inviscid, dstress, dFe, norm, dFlux); 522 523 for (int j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j]; 524 } // End Quadrature Point Loop 525 return 0; 526 } 527 528 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 529 return BoundaryIntegral_Jacobian(ctx, Q, in, out, StateFromU, StateFromU_fwd); 530 } 531 532 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 533 return BoundaryIntegral_Jacobian(ctx, Q, in, out, StateFromY, StateFromY_fwd); 534 } 535 536 // ***************************************************************************** 537 // Outflow boundary condition, weakly setting a constant pressure 538 // ***************************************************************************** 539 CEED_QFUNCTION_HELPER int PressureOutflow(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateFromQi_t StateFromQi, 540 StateFromQi_fwd_t StateFromQi_fwd) { 541 // Inputs 542 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 543 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 544 // Outputs 545 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], (*jac_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[1]; 546 547 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 548 const bool implicit = context->is_implicit; 549 const CeedScalar P0 = context->P0; 550 551 CeedPragmaSIMD 552 // Quadrature Point Loop 553 for (CeedInt i = 0; i < Q; i++) { 554 // Setup 555 // -- Interp in 556 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 557 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 558 State s = StateFromQi(context, qi, x_i); 559 s.Y.pressure = P0; 560 561 // -- Interp-to-Interp q_data 562 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 563 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 564 // We can effect this by swapping the sign on this weight 565 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 566 567 // ---- Normal vector 568 const CeedScalar norm[3] = {q_data_sur[1][i], q_data_sur[2][i], q_data_sur[3][i]}; 569 570 const CeedScalar dXdx[2][3] = { 571 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 572 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 573 }; 574 575 State grad_s[3]; 576 for (CeedInt j = 0; j < 3; j++) { 577 CeedScalar dx_i[3] = {0}, dqi[5]; 578 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]; 579 dx_i[j] = 1.; 580 grad_s[j] = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 581 } 582 583 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 584 KMStrainRate(grad_s, strain_rate); 585 NewtonianStress(context, strain_rate, kmstress); 586 KMUnpack(kmstress, stress); 587 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 588 589 StateConservative F_inviscid[3]; 590 FluxInviscid(context, s, F_inviscid); 591 592 CeedScalar Flux[5]; 593 FluxTotal_Boundary(F_inviscid, stress, Fe, norm, Flux); 594 595 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j]; 596 597 // Save values for Jacobian 598 for (int j = 0; j < 5; j++) jac_data_sur[j][i] = qi[j]; 599 for (int j = 0; j < 6; j++) jac_data_sur[5 + j][i] = kmstress[j]; 600 } // End Quadrature Point Loop 601 return 0; 602 } 603 604 CEED_QFUNCTION(PressureOutflow_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 605 return PressureOutflow(ctx, Q, in, out, StateFromU, StateFromU_fwd); 606 } 607 608 CEED_QFUNCTION(PressureOutflow_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 609 return PressureOutflow(ctx, Q, in, out, StateFromY, StateFromY_fwd); 610 } 611 612 // ***************************************************************************** 613 // Jacobian for weak-pressure outflow boundary condition 614 // ***************************************************************************** 615 CEED_QFUNCTION_HELPER int PressureOutflow_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 616 StateFromQi_t StateFromQi, StateFromQi_fwd_t StateFromQi_fwd) { 617 // Inputs 618 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 619 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 620 (*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 621 // Outputs 622 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 623 624 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 625 const bool implicit = context->is_implicit; 626 627 CeedPragmaSIMD 628 // Quadrature Point Loop 629 for (CeedInt i = 0; i < Q; i++) { 630 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 631 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 632 const CeedScalar norm[3] = {q_data_sur[1][i], q_data_sur[2][i], q_data_sur[3][i]}; 633 const CeedScalar dXdx[2][3] = { 634 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 635 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 636 }; 637 638 CeedScalar qi[5], kmstress[6], dqi[5], dx_i[3] = {0.}; 639 for (int j = 0; j < 5; j++) qi[j] = jac_data_sur[j][i]; 640 for (int j = 0; j < 6; j++) kmstress[j] = jac_data_sur[5 + j][i]; 641 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 642 643 State s = StateFromQi(context, qi, x_i); 644 State ds = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 645 s.Y.pressure = context->P0; 646 ds.Y.pressure = 0.; 647 648 State grad_ds[3]; 649 for (CeedInt j = 0; j < 3; j++) { 650 CeedScalar dx_i[3] = {0}, dqi_j[5]; 651 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]; 652 dx_i[j] = 1.; 653 grad_ds[j] = StateFromQi_fwd(context, s, dqi_j, x_i, dx_i); 654 } 655 656 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 657 KMStrainRate(grad_ds, dstrain_rate); 658 NewtonianStress(context, dstrain_rate, dkmstress); 659 KMUnpack(dkmstress, dstress); 660 KMUnpack(kmstress, stress); 661 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 662 663 StateConservative dF_inviscid[3]; 664 FluxInviscid_fwd(context, s, ds, dF_inviscid); 665 666 CeedScalar dFlux[5]; 667 FluxTotal_Boundary(dF_inviscid, dstress, dFe, norm, dFlux); 668 669 for (int j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j]; 670 } // End Quadrature Point Loop 671 return 0; 672 } 673 674 CEED_QFUNCTION(PressureOutflow_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 675 return PressureOutflow_Jacobian(ctx, Q, in, out, StateFromU, StateFromU_fwd); 676 } 677 678 CEED_QFUNCTION(PressureOutflow_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 679 return PressureOutflow_Jacobian(ctx, Q, in, out, StateFromY, StateFromY_fwd); 680 } 681 682 #endif // newtonian_h 683