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