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