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