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, StateVariable state_var) { 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 StateToQ(&context->gas, s, q, state_var); 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, STATEVAR_PRIMITIVE); 61 } 62 CEED_QFUNCTION(ICsNewtonianIG_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 63 return ICsNewtonianIG(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 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 k = 0; k < 3; k++) { 145 CeedScalar dx_i[3] = {0}, dU[5]; 146 for (CeedInt j = 0; j < 5; j++) dU[j] = Grad_q[0][j][i] * dXdx[0][k] + Grad_q[1][j][i] * dXdx[1][k] + Grad_q[2][j][i] * dXdx[2][k]; 147 dx_i[k] = 1.; 148 grad_s[k] = 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_State(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 < 5; j++) { 165 for (CeedInt k = 0; k < 3; k++) Grad_v[k][j][i] = wdetJ * (dXdx[k][0] * Flux[j][0] + dXdx[k][1] * Flux[j][1] + dXdx[k][2] * Flux[j][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, StateVariable state_var) { 193 // Inputs 194 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 195 const CeedScalar(*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 196 const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 197 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 198 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 199 200 // Outputs 201 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 202 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 203 CeedScalar(*jac_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[2]; 204 205 // Context 206 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 207 const CeedScalar *g = context->g; 208 const CeedScalar dt = context->dt; 209 const CeedScalar P0 = context->P0; 210 211 // Quadrature Point Loop 212 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 213 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 214 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 215 const State s = StateFromQ(context, qi, x_i, state_var); 216 217 // -- Interp-to-Interp q_data 218 const CeedScalar wdetJ = q_data[0][i]; 219 // -- Interp-to-Grad q_data 220 // ---- Inverse of change of coordinate matrix: X_i,j 221 const CeedScalar dXdx[3][3] = { 222 {q_data[1][i], q_data[2][i], q_data[3][i]}, 223 {q_data[4][i], q_data[5][i], q_data[6][i]}, 224 {q_data[7][i], q_data[8][i], q_data[9][i]} 225 }; 226 State grad_s[3]; 227 for (CeedInt k = 0; k < 3; k++) { 228 CeedScalar dx_i[3] = {0}, dqi[5]; 229 for (CeedInt j = 0; j < 5; j++) { 230 dqi[j] = Grad_q[0][j][i] * dXdx[0][k] + Grad_q[1][j][i] * dXdx[1][k] + Grad_q[2][j][i] * dXdx[2][k]; 231 } 232 dx_i[k] = 1.; 233 grad_s[k] = StateFromQ_fwd(context, s, dqi, x_i, dx_i, state_var); 234 } 235 236 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 237 KMStrainRate_State(grad_s, strain_rate); 238 NewtonianStress(context, strain_rate, kmstress); 239 KMUnpack(kmstress, stress); 240 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 241 242 StateConservative F_inviscid[3]; 243 FluxInviscid(context, s, F_inviscid); 244 245 // Total flux 246 CeedScalar Flux[5][3]; 247 FluxTotal(F_inviscid, stress, Fe, Flux); 248 249 for (CeedInt j = 0; j < 5; j++) { 250 for (CeedInt k = 0; k < 3; k++) { 251 Grad_v[k][j][i] = -wdetJ * (dXdx[k][0] * Flux[j][0] + dXdx[k][1] * Flux[j][1] + dXdx[k][2] * Flux[j][2]); 252 } 253 } 254 255 const CeedScalar body_force[5] = {0, s.U.density * g[0], s.U.density * g[1], s.U.density * g[2], 0}; 256 257 // -- Stabilization method: none (Galerkin), SU, or SUPG 258 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}, qi_dot[5], dx0[3] = {0}; 259 for (int j = 0; j < 5; j++) qi_dot[j] = q_dot[j][i]; 260 State s_dot = StateFromQ_fwd(context, s, qi_dot, x_i, dx0, state_var); 261 UnpackState_U(s_dot.U, U_dot); 262 263 for (CeedInt j = 0; j < 5; j++) v[j][i] = wdetJ * (U_dot[j] - body_force[j]); 264 if (context->idl_enable) { 265 CeedScalar damp_state[5] = {s.Y.pressure - P0, 0, 0, 0, 0}, idl_residual[5] = {0.}; 266 InternalDampingLayer(context, s, x_i, damp_state, idl_residual); 267 for (int j = 0; j < 5; j++) v[j][i] += wdetJ * idl_residual[j]; 268 } 269 270 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 271 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, x_i, stab); 272 273 for (CeedInt j = 0; j < 5; j++) { 274 for (CeedInt k = 0; k < 3; k++) { 275 Grad_v[k][j][i] += wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXdx[k][1] + stab[j][2] * dXdx[k][2]); 276 } 277 } 278 for (CeedInt j = 0; j < 5; j++) jac_data[j][i] = qi[j]; 279 for (CeedInt j = 0; j < 6; j++) jac_data[5 + j][i] = kmstress[j]; 280 for (CeedInt j = 0; j < 3; j++) jac_data[5 + 6 + j][i] = Tau_d[j]; 281 282 } // End Quadrature Point Loop 283 284 // Return 285 return 0; 286 } 287 288 CEED_QFUNCTION(IFunction_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 289 return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 290 } 291 292 CEED_QFUNCTION(IFunction_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 293 return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 294 } 295 296 // ***************************************************************************** 297 // This QFunction implements the jacobian of the Navier-Stokes equations for implicit time stepping method. 298 // ***************************************************************************** 299 CEED_QFUNCTION_HELPER int IJacobian_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 300 // Inputs 301 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 302 const CeedScalar(*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 303 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 304 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 305 const CeedScalar(*jac_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 306 307 // Outputs 308 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 309 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 310 311 // Context 312 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 313 const CeedScalar *g = context->g; 314 315 // Quadrature Point Loop 316 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 317 // -- Interp-to-Interp q_data 318 const CeedScalar wdetJ = q_data[0][i]; 319 // -- Interp-to-Grad q_data 320 // ---- Inverse of change of coordinate matrix: X_i,j 321 const CeedScalar dXdx[3][3] = { 322 {q_data[1][i], q_data[2][i], q_data[3][i]}, 323 {q_data[4][i], q_data[5][i], q_data[6][i]}, 324 {q_data[7][i], q_data[8][i], q_data[9][i]} 325 }; 326 327 CeedScalar qi[5], kmstress[6], Tau_d[3]; 328 for (int j = 0; j < 5; j++) qi[j] = jac_data[j][i]; 329 for (int j = 0; j < 6; j++) kmstress[j] = jac_data[5 + j][i]; 330 for (int j = 0; j < 3; j++) Tau_d[j] = jac_data[5 + 6 + j][i]; 331 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 332 State s = StateFromQ(context, qi, x_i, state_var); 333 334 CeedScalar dqi[5], dx0[3] = {0}; 335 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 336 State ds = StateFromQ_fwd(context, s, dqi, x_i, dx0, state_var); 337 338 State grad_ds[3]; 339 for (int k = 0; k < 3; k++) { 340 CeedScalar dqi_j[5]; 341 for (int j = 0; j < 5; j++) dqi_j[j] = Grad_dq[0][j][i] * dXdx[0][k] + Grad_dq[1][j][i] * dXdx[1][k] + Grad_dq[2][j][i] * dXdx[2][k]; 342 grad_ds[k] = StateFromQ_fwd(context, s, dqi_j, x_i, dx0, state_var); 343 } 344 345 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 346 KMStrainRate_State(grad_ds, dstrain_rate); 347 NewtonianStress(context, dstrain_rate, dkmstress); 348 KMUnpack(dkmstress, dstress); 349 KMUnpack(kmstress, stress); 350 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 351 352 StateConservative dF_inviscid[3]; 353 FluxInviscid_fwd(context, s, ds, dF_inviscid); 354 355 // Total flux 356 CeedScalar dFlux[5][3]; 357 FluxTotal(dF_inviscid, dstress, dFe, dFlux); 358 359 for (int j = 0; j < 5; j++) { 360 for (int k = 0; k < 3; k++) Grad_v[k][j][i] = -wdetJ * (dXdx[k][0] * dFlux[j][0] + dXdx[k][1] * dFlux[j][1] + dXdx[k][2] * dFlux[j][2]); 361 } 362 363 const CeedScalar dbody_force[5] = {0, ds.U.density * g[0], ds.U.density * g[1], ds.U.density * g[2], 0}; 364 CeedScalar dU[5] = {0.}; 365 UnpackState_U(ds.U, dU); 366 for (int j = 0; j < 5; j++) v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]); 367 368 if (context->idl_enable) { 369 CeedScalar damp_state[5] = {ds.Y.pressure, 0, 0, 0, 0}, idl_residual[5] = {0.}; 370 // This is a Picard-type linearization of the damping and could be replaced by an InternalDampingLayer_fwd that uses s and ds. 371 InternalDampingLayer(context, s, x_i, damp_state, idl_residual); 372 for (int j = 0; j < 5; j++) v[j][i] += wdetJ * idl_residual[j]; 373 } 374 375 // -- Stabilization method: none (Galerkin), SU, or SUPG 376 CeedScalar dstab[5][3], U_dot[5] = {0}; 377 for (CeedInt j = 0; j < 5; j++) U_dot[j] = context->ijacobian_time_shift * dU[j]; 378 Stabilization(context, s, Tau_d, grad_ds, U_dot, dbody_force, x_i, dstab); 379 380 for (int j = 0; j < 5; j++) { 381 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]); 382 } 383 } // End Quadrature Point Loop 384 return 0; 385 } 386 387 CEED_QFUNCTION(IJacobian_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 388 return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 389 } 390 391 CEED_QFUNCTION(IJacobian_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 392 return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 393 } 394 395 // ***************************************************************************** 396 // Compute boundary integral (ie. for strongly set inflows) 397 // ***************************************************************************** 398 CEED_QFUNCTION_HELPER int BoundaryIntegral(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 399 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 400 const CeedScalar(*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 401 const CeedScalar(*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 402 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 403 404 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 405 CeedScalar(*jac_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[1]; 406 407 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 408 const bool is_implicit = context->is_implicit; 409 410 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 411 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 412 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 413 State s = StateFromQ(context, qi, x_i, state_var); 414 415 const CeedScalar wdetJb = (is_implicit ? -1. : 1.) * q_data_sur[0][i]; 416 // ---- Normal vector 417 const CeedScalar norm[3] = {q_data_sur[1][i], q_data_sur[2][i], q_data_sur[3][i]}; 418 419 const CeedScalar dXdx[2][3] = { 420 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 421 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 422 }; 423 424 State grad_s[3]; 425 for (CeedInt k = 0; k < 3; k++) { 426 CeedScalar dx_i[3] = {0}, dqi[5]; 427 for (CeedInt j = 0; j < 5; j++) dqi[j] = Grad_q[0][j][i] * dXdx[0][k] + Grad_q[1][j][i] * dXdx[1][k]; 428 dx_i[k] = 1.; 429 grad_s[k] = StateFromQ_fwd(context, s, dqi, x_i, dx_i, state_var); 430 } 431 432 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 433 KMStrainRate_State(grad_s, strain_rate); 434 NewtonianStress(context, strain_rate, kmstress); 435 KMUnpack(kmstress, stress); 436 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 437 438 StateConservative F_inviscid[3]; 439 FluxInviscid(context, s, F_inviscid); 440 441 CeedScalar Flux[5]; 442 FluxTotal_Boundary(F_inviscid, stress, Fe, norm, Flux); 443 444 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j]; 445 446 for (int j = 0; j < 5; j++) jac_data_sur[j][i] = qi[j]; 447 for (int j = 0; j < 6; j++) jac_data_sur[5 + j][i] = kmstress[j]; 448 } 449 return 0; 450 } 451 452 CEED_QFUNCTION(BoundaryIntegral_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 453 return BoundaryIntegral(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 454 } 455 456 CEED_QFUNCTION(BoundaryIntegral_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 457 return BoundaryIntegral(ctx, Q, in, out, STATEVAR_PRIMITIVE); 458 } 459 460 // ***************************************************************************** 461 // Jacobian for "set nothing" boundary integral 462 // ***************************************************************************** 463 CEED_QFUNCTION_HELPER int BoundaryIntegral_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 464 StateVariable state_var) { 465 // Inputs 466 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 467 const CeedScalar(*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 468 const CeedScalar(*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 469 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 470 const CeedScalar(*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 471 472 // Outputs 473 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 474 475 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 476 const bool implicit = context->is_implicit; 477 478 // Quadrature Point Loop 479 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 480 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 481 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 482 const CeedScalar norm[3] = {q_data_sur[1][i], q_data_sur[2][i], q_data_sur[3][i]}; 483 const CeedScalar dXdx[2][3] = { 484 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 485 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 486 }; 487 488 CeedScalar qi[5], kmstress[6], dqi[5], dx_i[3] = {0.}; 489 for (int j = 0; j < 5; j++) qi[j] = jac_data_sur[j][i]; 490 for (int j = 0; j < 6; j++) kmstress[j] = jac_data_sur[5 + j][i]; 491 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 492 493 State s = StateFromQ(context, qi, x_i, state_var); 494 State ds = StateFromQ_fwd(context, s, dqi, x_i, dx_i, state_var); 495 496 State grad_ds[3]; 497 for (CeedInt k = 0; k < 3; k++) { 498 CeedScalar dx_i[3] = {0}, dqi_j[5]; 499 for (CeedInt j = 0; j < 5; j++) dqi_j[j] = Grad_dq[0][j][i] * dXdx[0][k] + Grad_dq[1][j][i] * dXdx[1][k]; 500 dx_i[k] = 1.; 501 grad_ds[k] = StateFromQ_fwd(context, s, dqi_j, x_i, dx_i, state_var); 502 } 503 504 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 505 KMStrainRate_State(grad_ds, dstrain_rate); 506 NewtonianStress(context, dstrain_rate, dkmstress); 507 KMUnpack(dkmstress, dstress); 508 KMUnpack(kmstress, stress); 509 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 510 511 StateConservative dF_inviscid[3]; 512 FluxInviscid_fwd(context, s, ds, dF_inviscid); 513 514 CeedScalar dFlux[5]; 515 FluxTotal_Boundary(dF_inviscid, dstress, dFe, norm, dFlux); 516 517 for (int j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j]; 518 } // End Quadrature Point Loop 519 return 0; 520 } 521 522 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 523 return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 524 } 525 526 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 527 return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 528 } 529 530 #endif // newtonian_h 531