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 CEED_QFUNCTION_HELPER void InternalDampingLayer(const NewtonianIdealGasContext context, const State s, const CeedScalar x_i[3], CeedScalar damp_Y[5], 24 CeedScalar damp_residual[5]) { 25 const CeedScalar sigma = LinearRampCoefficient(context->idl_amplitude, context->idl_length, context->idl_start, x_i[0]); 26 ScaleN(damp_Y, sigma, 5); 27 CeedScalar dx_i[3] = {0}; 28 State damp_s = StateFromY_fwd(context, s, damp_Y, x_i, dx_i); 29 30 CeedScalar U[5]; 31 UnpackState_U(damp_s.U, U); 32 for (int i = 0; i < 5; i++) damp_residual[i] += U[i]; 33 } 34 35 // ***************************************************************************** 36 // This QFunction sets a "still" initial condition for generic Newtonian IG problems 37 // ***************************************************************************** 38 CEED_QFUNCTION_HELPER int ICsNewtonianIG(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateToQi_t StateToQi) { 39 // Inputs 40 const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 41 42 // Outputs 43 CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 44 45 // Context 46 const SetupContext context = (SetupContext)ctx; 47 48 // Quadrature Point Loop 49 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 50 CeedScalar x[3] = {X[0][i], X[1][i], X[2][i]}; 51 CeedScalar q[5] = {0.}; 52 State s = StateFromPrimitive(&context->gas, context->reference, x); 53 StateToQi(&context->gas, s, q); 54 for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j]; 55 } // End of Quadrature Point Loop 56 return 0; 57 } 58 59 CEED_QFUNCTION(ICsNewtonianIG_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 60 return ICsNewtonianIG(ctx, Q, in, out, StateToY); 61 } 62 CEED_QFUNCTION(ICsNewtonianIG_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 63 return ICsNewtonianIG(ctx, Q, in, out, StateToU); 64 } 65 66 // ***************************************************************************** 67 // This QFunction implements the following formulation of Navier-Stokes with explicit time stepping method 68 // 69 // This is 3D compressible Navier-Stokes in conservation form with state variables of density, momentum density, and total energy density. 70 // 71 // State Variables: q = ( rho, U1, U2, U3, E ) 72 // rho - Mass Density 73 // Ui - Momentum Density, Ui = rho ui 74 // E - Total Energy Density, E = rho (cv T + (u u)/2 + g z) 75 // 76 // Navier-Stokes Equations: 77 // drho/dt + div( U ) = 0 78 // dU/dt + div( rho (u x u) + P I3 ) + rho g khat = div( Fu ) 79 // dE/dt + div( (E + P) u ) = div( Fe ) 80 // 81 // Viscous Stress: 82 // Fu = mu (grad( u ) + grad( u )^T + lambda div ( u ) I3) 83 // 84 // Thermal Stress: 85 // Fe = u Fu + k grad( T ) 86 // Equation of State 87 // P = (gamma - 1) (E - rho (u u) / 2 - rho g z) 88 // 89 // Stabilization: 90 // Tau = diag(TauC, TauM, TauM, TauM, TauE) 91 // f1 = rho sqrt(ui uj gij) 92 // gij = dXi/dX * dXi/dX 93 // TauC = Cc f1 / (8 gii) 94 // TauM = min( 1 , 1 / f1 ) 95 // TauE = TauM / (Ce cv) 96 // 97 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 98 // 99 // Constants: 100 // lambda = - 2 / 3, From Stokes hypothesis 101 // mu , Dynamic viscosity 102 // k , Thermal conductivity 103 // cv , Specific heat, constant volume 104 // cp , Specific heat, constant pressure 105 // g , Gravity 106 // gamma = cp / cv, Specific heat ratio 107 // 108 // 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 109 // gradu ) 110 // ***************************************************************************** 111 CEED_QFUNCTION(RHSFunction_Newtonian)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 112 // Inputs 113 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 114 const CeedScalar(*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 115 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 116 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 117 118 // Outputs 119 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 120 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 121 122 // Context 123 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 124 const CeedScalar *g = context->g; 125 const CeedScalar dt = context->dt; 126 127 // Quadrature Point Loop 128 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 129 CeedScalar U[5]; 130 for (int j = 0; j < 5; j++) U[j] = q[j][i]; 131 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 132 State s = StateFromU(context, U, x_i); 133 134 // -- Interp-to-Interp q_data 135 const CeedScalar wdetJ = q_data[0][i]; 136 // -- Interp-to-Grad q_data 137 // ---- Inverse of change of coordinate matrix: X_i,j 138 const CeedScalar dXdx[3][3] = { 139 {q_data[1][i], q_data[2][i], q_data[3][i]}, 140 {q_data[4][i], q_data[5][i], q_data[6][i]}, 141 {q_data[7][i], q_data[8][i], q_data[9][i]} 142 }; 143 State grad_s[3]; 144 for (CeedInt j = 0; j < 3; j++) { 145 CeedScalar dx_i[3] = {0}, dU[5]; 146 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]; 147 dx_i[j] = 1.; 148 grad_s[j] = StateFromU_fwd(context, s, dU, x_i, dx_i); 149 } 150 151 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 152 KMStrainRate(grad_s, strain_rate); 153 NewtonianStress(context, strain_rate, kmstress); 154 KMUnpack(kmstress, stress); 155 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 156 157 StateConservative F_inviscid[3]; 158 FluxInviscid(context, s, F_inviscid); 159 160 // Total flux 161 CeedScalar Flux[5][3]; 162 FluxTotal(F_inviscid, stress, Fe, Flux); 163 164 for (CeedInt j = 0; j < 3; j++) { 165 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]); 166 } 167 168 const CeedScalar body_force[5] = {0, s.U.density * g[0], s.U.density * g[1], s.U.density * g[2], 0}; 169 for (int j = 0; j < 5; j++) v[j][i] = wdetJ * body_force[j]; 170 171 // -- Stabilization method: none (Galerkin), SU, or SUPG 172 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}; 173 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 174 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, x_i, stab); 175 176 for (CeedInt j = 0; j < 5; j++) { 177 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]); 178 } 179 } // End Quadrature Point Loop 180 181 // Return 182 return 0; 183 } 184 185 // ***************************************************************************** 186 // This QFunction implements the Navier-Stokes equations (mentioned above) with implicit time stepping method 187 // 188 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 189 // SUPG = Galerkin + grad(v) . ( Ai^T * Tau * (q_dot + Aj q,j - body force) ) 190 // (diffusive terms will be added later) 191 // ***************************************************************************** 192 CEED_QFUNCTION_HELPER int IFunction_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateFromQi_t StateFromQi, 193 StateFromQi_fwd_t StateFromQi_fwd) { 194 // Inputs 195 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 196 const CeedScalar(*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 197 const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 198 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 199 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 200 201 // Outputs 202 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 203 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 204 CeedScalar(*jac_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[2]; 205 206 // Context 207 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 208 const CeedScalar *g = context->g; 209 const CeedScalar dt = context->dt; 210 const CeedScalar P0 = context->P0; 211 212 // Quadrature Point Loop 213 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 214 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 215 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 216 const State s = StateFromQi(context, qi, x_i); 217 218 // -- Interp-to-Interp q_data 219 const CeedScalar wdetJ = q_data[0][i]; 220 // -- Interp-to-Grad q_data 221 // ---- Inverse of change of coordinate matrix: X_i,j 222 const CeedScalar dXdx[3][3] = { 223 {q_data[1][i], q_data[2][i], q_data[3][i]}, 224 {q_data[4][i], q_data[5][i], q_data[6][i]}, 225 {q_data[7][i], q_data[8][i], q_data[9][i]} 226 }; 227 State grad_s[3]; 228 for (CeedInt j = 0; j < 3; j++) { 229 CeedScalar dx_i[3] = {0}, dqi[5]; 230 for (CeedInt k = 0; k < 5; k++) { 231 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]; 232 } 233 dx_i[j] = 1.; 234 grad_s[j] = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 235 } 236 237 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 238 KMStrainRate(grad_s, strain_rate); 239 NewtonianStress(context, strain_rate, kmstress); 240 KMUnpack(kmstress, stress); 241 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 242 243 StateConservative F_inviscid[3]; 244 FluxInviscid(context, s, F_inviscid); 245 246 // Total flux 247 CeedScalar Flux[5][3]; 248 FluxTotal(F_inviscid, stress, Fe, Flux); 249 250 for (CeedInt j = 0; j < 3; j++) { 251 for (CeedInt k = 0; k < 5; k++) { 252 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * Flux[k][0] + dXdx[j][1] * Flux[k][1] + dXdx[j][2] * Flux[k][2]); 253 } 254 } 255 256 const CeedScalar body_force[5] = {0, s.U.density * g[0], s.U.density * g[1], s.U.density * g[2], 0}; 257 258 // -- Stabilization method: none (Galerkin), SU, or SUPG 259 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}, qi_dot[5], dx0[3] = {0}; 260 for (int j = 0; j < 5; j++) qi_dot[j] = q_dot[j][i]; 261 State s_dot = StateFromQi_fwd(context, s, qi_dot, x_i, dx0); 262 UnpackState_U(s_dot.U, U_dot); 263 264 for (CeedInt j = 0; j < 5; j++) v[j][i] = wdetJ * (U_dot[j] - body_force[j]); 265 if (context->idl_enable) { 266 CeedScalar damp_state[5] = {s.Y.pressure - P0, 0, 0, 0, 0}, idl_residual[5] = {0.}; 267 InternalDampingLayer(context, s, x_i, damp_state, idl_residual); 268 for (int j = 0; j < 5; j++) v[j][i] += wdetJ * idl_residual[j]; 269 } 270 271 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 272 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, x_i, stab); 273 274 for (CeedInt j = 0; j < 5; j++) { 275 for (CeedInt k = 0; k < 3; k++) { 276 Grad_v[k][j][i] += wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXdx[k][1] + stab[j][2] * dXdx[k][2]); 277 } 278 } 279 for (CeedInt j = 0; j < 5; j++) jac_data[j][i] = qi[j]; 280 for (CeedInt j = 0; j < 6; j++) jac_data[5 + j][i] = kmstress[j]; 281 for (CeedInt j = 0; j < 3; j++) jac_data[5 + 6 + j][i] = Tau_d[j]; 282 283 } // End Quadrature Point Loop 284 285 // Return 286 return 0; 287 } 288 289 CEED_QFUNCTION(IFunction_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 290 return IFunction_Newtonian(ctx, Q, in, out, StateFromU, StateFromU_fwd); 291 } 292 293 CEED_QFUNCTION(IFunction_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 294 return IFunction_Newtonian(ctx, Q, in, out, StateFromY, StateFromY_fwd); 295 } 296 297 // ***************************************************************************** 298 // This QFunction implements the jacobian of the Navier-Stokes equations for implicit time stepping method. 299 // ***************************************************************************** 300 CEED_QFUNCTION_HELPER int IJacobian_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateFromQi_t StateFromQi, 301 StateFromQi_fwd_t StateFromQi_fwd) { 302 // Inputs 303 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 304 const CeedScalar(*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 305 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 306 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 307 const CeedScalar(*jac_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 308 309 // Outputs 310 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 311 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 312 313 // Context 314 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 315 const CeedScalar *g = context->g; 316 317 // Quadrature Point Loop 318 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 319 // -- Interp-to-Interp q_data 320 const CeedScalar wdetJ = q_data[0][i]; 321 // -- Interp-to-Grad q_data 322 // ---- Inverse of change of coordinate matrix: X_i,j 323 const CeedScalar dXdx[3][3] = { 324 {q_data[1][i], q_data[2][i], q_data[3][i]}, 325 {q_data[4][i], q_data[5][i], q_data[6][i]}, 326 {q_data[7][i], q_data[8][i], q_data[9][i]} 327 }; 328 329 CeedScalar qi[5], kmstress[6], Tau_d[3]; 330 for (int j = 0; j < 5; j++) qi[j] = jac_data[j][i]; 331 for (int j = 0; j < 6; j++) kmstress[j] = jac_data[5 + j][i]; 332 for (int j = 0; j < 3; j++) Tau_d[j] = jac_data[5 + 6 + j][i]; 333 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 334 State s = StateFromQi(context, qi, x_i); 335 336 CeedScalar dqi[5], dx0[3] = {0}; 337 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 338 State ds = StateFromQi_fwd(context, s, dqi, x_i, dx0); 339 340 State grad_ds[3]; 341 for (int j = 0; j < 3; j++) { 342 CeedScalar dqi_j[5]; 343 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]; 344 grad_ds[j] = StateFromQi_fwd(context, s, dqi_j, x_i, dx0); 345 } 346 347 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 348 KMStrainRate(grad_ds, dstrain_rate); 349 NewtonianStress(context, dstrain_rate, dkmstress); 350 KMUnpack(dkmstress, dstress); 351 KMUnpack(kmstress, stress); 352 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 353 354 StateConservative dF_inviscid[3]; 355 FluxInviscid_fwd(context, s, ds, dF_inviscid); 356 357 // Total flux 358 CeedScalar dFlux[5][3]; 359 FluxTotal(dF_inviscid, dstress, dFe, dFlux); 360 361 for (int j = 0; j < 3; j++) { 362 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]); 363 } 364 365 const CeedScalar dbody_force[5] = {0, ds.U.density * g[0], ds.U.density * g[1], ds.U.density * g[2], 0}; 366 CeedScalar dU[5] = {0.}; 367 UnpackState_U(ds.U, dU); 368 for (int j = 0; j < 5; j++) v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]); 369 370 if (context->idl_enable) { 371 CeedScalar damp_state[5] = {ds.Y.pressure, 0, 0, 0, 0}, idl_residual[5] = {0.}; 372 // This is a Picard-type linearization of the damping and could be replaced by an InternalDampingLayer_fwd that uses s and ds. 373 InternalDampingLayer(context, s, x_i, damp_state, idl_residual); 374 for (int j = 0; j < 5; j++) v[j][i] += wdetJ * idl_residual[j]; 375 } 376 377 // -- Stabilization method: none (Galerkin), SU, or SUPG 378 CeedScalar dstab[5][3], U_dot[5] = {0}; 379 for (CeedInt j = 0; j < 5; j++) U_dot[j] = context->ijacobian_time_shift * dU[j]; 380 Stabilization(context, s, Tau_d, grad_ds, U_dot, dbody_force, x_i, dstab); 381 382 for (int j = 0; j < 5; j++) { 383 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]); 384 } 385 } // End Quadrature Point Loop 386 return 0; 387 } 388 389 CEED_QFUNCTION(IJacobian_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 390 return IJacobian_Newtonian(ctx, Q, in, out, StateFromU, StateFromU_fwd); 391 } 392 393 CEED_QFUNCTION(IJacobian_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 394 return IJacobian_Newtonian(ctx, Q, in, out, StateFromY, StateFromY_fwd); 395 } 396 397 // ***************************************************************************** 398 // Compute boundary integral (ie. for strongly set inflows) 399 // ***************************************************************************** 400 CEED_QFUNCTION_HELPER int BoundaryIntegral(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateFromQi_t StateFromQi, 401 StateFromQi_fwd_t StateFromQi_fwd) { 402 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 403 const CeedScalar(*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 404 const CeedScalar(*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 405 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 406 407 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 408 CeedScalar(*jac_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[1]; 409 410 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 411 const bool is_implicit = context->is_implicit; 412 413 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 414 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 415 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 416 State s = StateFromQi(context, qi, x_i); 417 418 const CeedScalar wdetJb = (is_implicit ? -1. : 1.) * q_data_sur[0][i]; 419 // ---- Normal vector 420 const CeedScalar norm[3] = {q_data_sur[1][i], q_data_sur[2][i], q_data_sur[3][i]}; 421 422 const CeedScalar dXdx[2][3] = { 423 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 424 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 425 }; 426 427 State grad_s[3]; 428 for (CeedInt j = 0; j < 3; j++) { 429 CeedScalar dx_i[3] = {0}, dqi[5]; 430 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]; 431 dx_i[j] = 1.; 432 grad_s[j] = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 433 } 434 435 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 436 KMStrainRate(grad_s, strain_rate); 437 NewtonianStress(context, strain_rate, kmstress); 438 KMUnpack(kmstress, stress); 439 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 440 441 StateConservative F_inviscid[3]; 442 FluxInviscid(context, s, F_inviscid); 443 444 CeedScalar Flux[5]; 445 FluxTotal_Boundary(F_inviscid, stress, Fe, norm, Flux); 446 447 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j]; 448 449 for (int j = 0; j < 5; j++) jac_data_sur[j][i] = qi[j]; 450 for (int j = 0; j < 6; j++) jac_data_sur[5 + j][i] = kmstress[j]; 451 } 452 return 0; 453 } 454 455 CEED_QFUNCTION(BoundaryIntegral_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 456 return BoundaryIntegral(ctx, Q, in, out, StateFromU, StateFromU_fwd); 457 } 458 459 CEED_QFUNCTION(BoundaryIntegral_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 460 return BoundaryIntegral(ctx, Q, in, out, StateFromY, StateFromY_fwd); 461 } 462 463 // ***************************************************************************** 464 // Jacobian for "set nothing" boundary integral 465 // ***************************************************************************** 466 CEED_QFUNCTION_HELPER int BoundaryIntegral_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 467 StateFromQi_t StateFromQi, StateFromQi_fwd_t StateFromQi_fwd) { 468 // Inputs 469 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 470 const CeedScalar(*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 471 const CeedScalar(*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 472 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 473 const CeedScalar(*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 474 475 // Outputs 476 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 477 478 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 479 const bool implicit = context->is_implicit; 480 481 // Quadrature Point Loop 482 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 483 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 484 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 485 const CeedScalar norm[3] = {q_data_sur[1][i], q_data_sur[2][i], q_data_sur[3][i]}; 486 const CeedScalar dXdx[2][3] = { 487 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 488 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 489 }; 490 491 CeedScalar qi[5], kmstress[6], dqi[5], dx_i[3] = {0.}; 492 for (int j = 0; j < 5; j++) qi[j] = jac_data_sur[j][i]; 493 for (int j = 0; j < 6; j++) kmstress[j] = jac_data_sur[5 + j][i]; 494 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 495 496 State s = StateFromQi(context, qi, x_i); 497 State ds = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 498 499 State grad_ds[3]; 500 for (CeedInt j = 0; j < 3; j++) { 501 CeedScalar dx_i[3] = {0}, dqi_j[5]; 502 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]; 503 dx_i[j] = 1.; 504 grad_ds[j] = StateFromQi_fwd(context, s, dqi_j, x_i, dx_i); 505 } 506 507 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 508 KMStrainRate(grad_ds, dstrain_rate); 509 NewtonianStress(context, dstrain_rate, dkmstress); 510 KMUnpack(dkmstress, dstress); 511 KMUnpack(kmstress, stress); 512 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 513 514 StateConservative dF_inviscid[3]; 515 FluxInviscid_fwd(context, s, ds, dF_inviscid); 516 517 CeedScalar dFlux[5]; 518 FluxTotal_Boundary(dF_inviscid, dstress, dFe, norm, dFlux); 519 520 for (int j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j]; 521 } // End Quadrature Point Loop 522 return 0; 523 } 524 525 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 526 return BoundaryIntegral_Jacobian(ctx, Q, in, out, StateFromU, StateFromU_fwd); 527 } 528 529 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 530 return BoundaryIntegral_Jacobian(ctx, Q, in, out, StateFromY, StateFromY_fwd); 531 } 532 533 #endif // newtonian_h 534