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.}, divFdiff[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, divFdiff, 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}, zeroFlux[5] = {0.}; 195 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 196 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, zeroFlux, 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 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 214 215 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 216 const CeedScalar(*Grad_q) = in[1]; 217 const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 218 const CeedScalar(*q_data) = in[3]; 219 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 220 const CeedScalar(*divFdiff)[CEED_Q_VLA] = context->divFdiff_method != DIV_DIFF_FLUX_PROJ_NONE ? (const CeedScalar(*)[CEED_Q_VLA])in[5] : NULL; 221 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 222 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 223 CeedScalar(*jac_data) = out[2]; 224 225 const CeedScalar *g = context->g; 226 const CeedScalar dt = context->dt; 227 const CeedScalar P0 = context->idl_pressure; 228 229 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 230 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 231 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 232 const State s = StateFromQ(context, qi, state_var); 233 234 CeedScalar wdetJ, dXdx[3][3]; 235 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 236 State grad_s[3]; 237 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 238 239 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 240 KMStrainRate_State(grad_s, strain_rate); 241 NewtonianStress(context, strain_rate, kmstress); 242 KMUnpack(kmstress, stress); 243 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 244 245 StateConservative F_inviscid[3]; 246 FluxInviscid(context, s, F_inviscid); 247 248 // Total flux 249 CeedScalar Flux[5][3]; 250 FluxTotal(F_inviscid, stress, Fe, Flux); 251 252 for (CeedInt j = 0; j < 5; j++) { 253 for (CeedInt k = 0; k < 3; k++) { 254 Grad_v[k][j][i] = -wdetJ * (dXdx[k][0] * Flux[j][0] + dXdx[k][1] * Flux[j][1] + dXdx[k][2] * Flux[j][2]); 255 } 256 } 257 258 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)}; 259 260 // -- Stabilization method: none (Galerkin), SU, or SUPG 261 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}, qi_dot[5]; 262 for (int j = 0; j < 5; j++) qi_dot[j] = q_dot[j][i]; 263 State s_dot = StateFromQ_fwd(context, s, qi_dot, state_var); 264 UnpackState_U(s_dot.U, U_dot); 265 266 for (CeedInt j = 0; j < 5; j++) v[j][i] = wdetJ * (U_dot[j] - body_force[j]); 267 if (context->idl_enable) { 268 const CeedScalar sigma = LinearRampCoefficient(context->idl_amplitude, context->idl_length, context->idl_start, x_i[0]); 269 StoredValuesPack(Q, i, 14, 1, &sigma, jac_data); 270 CeedScalar damp_state[5] = {s.Y.pressure - P0, 0, 0, 0, 0}, idl_residual[5] = {0.}; 271 InternalDampingLayer(context, s, sigma, damp_state, idl_residual); 272 for (int j = 0; j < 5; j++) v[j][i] += wdetJ * idl_residual[j]; 273 } 274 275 CeedScalar divFdiff_i[5] = {0.}; 276 if (context->divFdiff_method != DIV_DIFF_FLUX_PROJ_NONE) { 277 for (int j = 1; j < 5; j++) divFdiff_i[j] = divFdiff[j - 1][i]; 278 } 279 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 280 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, divFdiff_i, stab); 281 282 for (CeedInt j = 0; j < 5; j++) { 283 for (CeedInt k = 0; k < 3; k++) { 284 Grad_v[k][j][i] += wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXdx[k][1] + stab[j][2] * dXdx[k][2]); 285 } 286 } 287 StoredValuesPack(Q, i, 0, 5, qi, jac_data); 288 StoredValuesPack(Q, i, 5, 6, kmstress, jac_data); 289 StoredValuesPack(Q, i, 11, 3, Tau_d, jac_data); 290 } 291 return 0; 292 } 293 294 CEED_QFUNCTION(IFunction_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 295 return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 296 } 297 298 CEED_QFUNCTION(IFunction_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 299 return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 300 } 301 302 CEED_QFUNCTION(IFunction_Newtonian_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 303 return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_ENTROPY); 304 } 305 306 // ***************************************************************************** 307 // This QFunction implements the jacobian of the Navier-Stokes equations for implicit time stepping method. 308 // ***************************************************************************** 309 CEED_QFUNCTION_HELPER int IJacobian_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 310 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 311 const CeedScalar(*Grad_dq) = in[1]; 312 const CeedScalar(*q_data) = in[2]; 313 const CeedScalar(*jac_data) = in[3]; 314 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 315 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 316 317 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 318 const CeedScalar *g = context->g; 319 320 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 321 CeedScalar wdetJ, dXdx[3][3]; 322 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 323 324 CeedScalar qi[5], kmstress[6], Tau_d[3]; 325 StoredValuesUnpack(Q, i, 0, 5, jac_data, qi); 326 StoredValuesUnpack(Q, i, 5, 6, jac_data, kmstress); 327 StoredValuesUnpack(Q, i, 11, 3, jac_data, Tau_d); 328 State s = StateFromQ(context, qi, state_var); 329 330 CeedScalar dqi[5]; 331 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 332 State ds = StateFromQ_fwd(context, s, dqi, state_var); 333 334 State grad_ds[3]; 335 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_dq, dXdx, grad_ds); 336 337 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 338 KMStrainRate_State(grad_ds, dstrain_rate); 339 NewtonianStress(context, dstrain_rate, dkmstress); 340 KMUnpack(dkmstress, dstress); 341 KMUnpack(kmstress, stress); 342 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 343 344 StateConservative dF_inviscid[3]; 345 FluxInviscid_fwd(context, s, ds, dF_inviscid); 346 347 // Total flux 348 CeedScalar dFlux[5][3]; 349 FluxTotal(dF_inviscid, dstress, dFe, dFlux); 350 351 for (int j = 0; j < 5; j++) { 352 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]); 353 } 354 355 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)}; 356 CeedScalar dU[5] = {0.}; 357 UnpackState_U(ds.U, dU); 358 for (int j = 0; j < 5; j++) v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]); 359 360 if (context->idl_enable) { 361 const CeedScalar sigma = jac_data[14 * Q + i]; 362 CeedScalar damp_state[5] = {ds.Y.pressure, 0, 0, 0, 0}, idl_residual[5] = {0.}; 363 // This is a Picard-type linearization of the damping and could be replaced by an InternalDampingLayer_fwd that uses s and ds. 364 InternalDampingLayer(context, s, sigma, damp_state, idl_residual); 365 for (int j = 0; j < 5; j++) v[j][i] += wdetJ * idl_residual[j]; 366 } 367 368 // -- Stabilization method: none (Galerkin), SU, or SUPG 369 CeedScalar dstab[5][3], U_dot[5] = {0}; 370 for (CeedInt j = 0; j < 5; j++) U_dot[j] = context->ijacobian_time_shift * dU[j]; 371 const CeedScalar zeroFlux[5] = {0.}; 372 Stabilization(context, s, Tau_d, grad_ds, U_dot, dbody_force, zeroFlux, dstab); 373 374 for (int j = 0; j < 5; j++) { 375 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]); 376 } 377 } 378 return 0; 379 } 380 381 CEED_QFUNCTION(IJacobian_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 382 return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 383 } 384 385 CEED_QFUNCTION(IJacobian_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 386 return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 387 } 388 389 CEED_QFUNCTION(IJacobian_Newtonian_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 390 return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_ENTROPY); 391 } 392 393 // ***************************************************************************** 394 // Compute boundary integral (ie. for strongly set inflows) 395 // ***************************************************************************** 396 CEED_QFUNCTION_HELPER int BoundaryIntegral(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 397 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 398 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 399 const CeedScalar(*Grad_q) = in[1]; 400 const CeedScalar(*q_data_sur) = in[2]; 401 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 402 CeedScalar(*jac_data_sur) = context->is_implicit ? out[1] : NULL; 403 404 const bool is_implicit = context->is_implicit; 405 406 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 407 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 408 State s = StateFromQ(context, qi, state_var); 409 410 CeedScalar wdetJb, dXdx[2][3], normal[3]; 411 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, normal); 412 wdetJb *= is_implicit ? -1. : 1.; 413 414 State grad_s[3]; 415 StatePhysicalGradientFromReference_Boundary(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 416 417 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 418 KMStrainRate_State(grad_s, strain_rate); 419 NewtonianStress(context, strain_rate, kmstress); 420 KMUnpack(kmstress, stress); 421 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 422 423 StateConservative F_inviscid[3]; 424 FluxInviscid(context, s, F_inviscid); 425 426 CeedScalar Flux[5]; 427 FluxTotal_Boundary(F_inviscid, stress, Fe, normal, Flux); 428 429 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j]; 430 431 if (is_implicit) { 432 StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur); 433 StoredValuesPack(Q, i, 5, 6, kmstress, jac_data_sur); 434 } 435 } 436 return 0; 437 } 438 439 CEED_QFUNCTION(BoundaryIntegral_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 440 return BoundaryIntegral(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 441 } 442 443 CEED_QFUNCTION(BoundaryIntegral_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 444 return BoundaryIntegral(ctx, Q, in, out, STATEVAR_PRIMITIVE); 445 } 446 447 CEED_QFUNCTION(BoundaryIntegral_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 448 return BoundaryIntegral(ctx, Q, in, out, STATEVAR_ENTROPY); 449 } 450 451 // ***************************************************************************** 452 // Jacobian for "set nothing" boundary integral 453 // ***************************************************************************** 454 CEED_QFUNCTION_HELPER int BoundaryIntegral_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 455 StateVariable state_var) { 456 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 457 const CeedScalar(*Grad_dq) = in[1]; 458 const CeedScalar(*q_data_sur) = in[2]; 459 const CeedScalar(*jac_data_sur) = in[4]; 460 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 461 462 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 463 const bool is_implicit = context->is_implicit; 464 465 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 466 CeedScalar wdetJb, dXdx[2][3], normal[3]; 467 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, normal); 468 wdetJb *= is_implicit ? -1. : 1.; 469 470 CeedScalar qi[5], kmstress[6], dqi[5]; 471 StoredValuesUnpack(Q, i, 0, 5, jac_data_sur, qi); 472 StoredValuesUnpack(Q, i, 5, 6, jac_data_sur, kmstress); 473 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 474 475 State s = StateFromQ(context, qi, state_var); 476 State ds = StateFromQ_fwd(context, s, dqi, state_var); 477 478 State grad_ds[3]; 479 StatePhysicalGradientFromReference_Boundary(Q, i, context, s, state_var, Grad_dq, dXdx, grad_ds); 480 481 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 482 KMStrainRate_State(grad_ds, dstrain_rate); 483 NewtonianStress(context, dstrain_rate, dkmstress); 484 KMUnpack(dkmstress, dstress); 485 KMUnpack(kmstress, stress); 486 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 487 488 StateConservative dF_inviscid[3]; 489 FluxInviscid_fwd(context, s, ds, dF_inviscid); 490 491 CeedScalar dFlux[5]; 492 FluxTotal_Boundary(dF_inviscid, dstress, dFe, normal, dFlux); 493 494 for (int j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j]; 495 } 496 return 0; 497 } 498 499 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 500 return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 501 } 502 503 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 504 return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 505 } 506 507 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 508 return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_ENTROPY); 509 } 510 511 // @brief Volume integral for RHS of divergence of diffusive flux direct projection 512 CEED_QFUNCTION_HELPER int DivDiffusiveFluxVolumeRHS_NS(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 513 StateVariable state_var) { 514 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 515 const CeedScalar(*Grad_q) = in[1]; 516 const CeedScalar(*q_data) = in[2]; 517 CeedScalar(*Grad_v)[4][CEED_Q_VLA] = (CeedScalar(*)[4][CEED_Q_VLA])out[0]; 518 519 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 520 const StateConservative ZeroInviscidFluxes[3] = {{0}}; 521 522 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 523 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 524 const State s = StateFromQ(context, qi, state_var); 525 CeedScalar wdetJ, dXdx[3][3]; 526 CeedScalar stress[3][3], Fe[3], Fdiff[5][3]; 527 528 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 529 { // Get stress and Fe 530 State grad_s[3]; 531 CeedScalar strain_rate[6], kmstress[6]; 532 533 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 534 KMStrainRate_State(grad_s, strain_rate); 535 NewtonianStress(context, strain_rate, kmstress); 536 KMUnpack(kmstress, stress); 537 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 538 } 539 540 FluxTotal(ZeroInviscidFluxes, stress, Fe, Fdiff); 541 542 for (CeedInt j = 1; j < 5; j++) { // Continuity has no diffusive flux, therefore skip 543 for (CeedInt k = 0; k < 3; k++) { 544 Grad_v[k][j - 1][i] = -wdetJ * Dot3(dXdx[k], Fdiff[j]); 545 } 546 } 547 } 548 return 0; 549 } 550 551 CEED_QFUNCTION(DivDiffusiveFluxVolumeRHS_NS_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 552 return DivDiffusiveFluxVolumeRHS_NS(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 553 } 554 555 CEED_QFUNCTION(DivDiffusiveFluxVolumeRHS_NS_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 556 return DivDiffusiveFluxVolumeRHS_NS(ctx, Q, in, out, STATEVAR_PRIMITIVE); 557 } 558 559 CEED_QFUNCTION(DivDiffusiveFluxVolumeRHS_NS_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 560 return DivDiffusiveFluxVolumeRHS_NS(ctx, Q, in, out, STATEVAR_ENTROPY); 561 } 562 563 // @brief Boundary integral for RHS of divergence of diffusive flux direct projection 564 CEED_QFUNCTION_HELPER int DivDiffusiveFluxBoundaryRHS_NS(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 565 StateVariable state_var) { 566 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 567 const CeedScalar(*Grad_q) = in[1]; 568 const CeedScalar(*q_data) = in[2]; 569 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 570 571 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 572 const StateConservative ZeroInviscidFluxes[3] = {{0}}; 573 574 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 575 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 576 const State s = StateFromQ(context, qi, state_var); 577 CeedScalar wdetJ, dXdx[3][3], normal[3]; 578 CeedScalar stress[3][3], Fe[3], Fdiff[5]; 579 580 QdataBoundaryGradientUnpack_3D(Q, i, q_data, &wdetJ, dXdx, normal); 581 { // Get stress and Fe 582 State grad_s[3]; 583 CeedScalar strain_rate[6], kmstress[6]; 584 585 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 586 KMStrainRate_State(grad_s, strain_rate); 587 NewtonianStress(context, strain_rate, kmstress); 588 KMUnpack(kmstress, stress); 589 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 590 } 591 592 FluxTotal_Boundary(ZeroInviscidFluxes, stress, Fe, normal, Fdiff); 593 594 // Continuity has no diffusive flux, therefore skip 595 for (CeedInt j = 1; j < 5; j++) v[j - 1][i] = wdetJ * Fdiff[j]; 596 } 597 return 0; 598 } 599 600 CEED_QFUNCTION(DivDiffusiveFluxBoundaryRHS_NS_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 601 return DivDiffusiveFluxBoundaryRHS_NS(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 602 } 603 604 CEED_QFUNCTION(DivDiffusiveFluxBoundaryRHS_NS_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 605 return DivDiffusiveFluxBoundaryRHS_NS(ctx, Q, in, out, STATEVAR_PRIMITIVE); 606 } 607 608 CEED_QFUNCTION(DivDiffusiveFluxBoundaryRHS_NS_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 609 return DivDiffusiveFluxBoundaryRHS_NS(ctx, Q, in, out, STATEVAR_ENTROPY); 610 } 611 612 // @brief Integral for RHS of diffusive flux indirect projection 613 CEED_QFUNCTION_HELPER int DiffusiveFluxRHS_NS(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 614 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 615 const CeedScalar(*Grad_q) = in[1]; 616 const CeedScalar(*q_data) = in[2]; 617 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 618 619 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 620 const StateConservative ZeroInviscidFluxes[3] = {{0}}; 621 622 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 623 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 624 const State s = StateFromQ(context, qi, state_var); 625 CeedScalar wdetJ, dXdx[3][3]; 626 CeedScalar stress[3][3], Fe[3], Fdiff[5][3]; 627 628 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 629 { // Get stress and Fe 630 State grad_s[3]; 631 CeedScalar strain_rate[6], kmstress[6]; 632 633 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 634 KMStrainRate_State(grad_s, strain_rate); 635 NewtonianStress(context, strain_rate, kmstress); 636 KMUnpack(kmstress, stress); 637 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 638 } 639 640 FluxTotal(ZeroInviscidFluxes, stress, Fe, Fdiff); 641 642 for (CeedInt j = 1; j < 5; j++) { // Continuity has no diffusive flux, therefore skip 643 for (CeedInt k = 0; k < 3; k++) { 644 v[(j - 1) * 3 + k][i] = wdetJ * Fdiff[j][k]; 645 } 646 } 647 } 648 return 0; 649 } 650 651 CEED_QFUNCTION(DiffusiveFluxRHS_NS_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 652 return DiffusiveFluxRHS_NS(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 653 } 654 655 CEED_QFUNCTION(DiffusiveFluxRHS_NS_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 656 return DiffusiveFluxRHS_NS(ctx, Q, in, out, STATEVAR_PRIMITIVE); 657 } 658 659 CEED_QFUNCTION(DiffusiveFluxRHS_NS_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 660 return DiffusiveFluxRHS_NS(ctx, Q, in, out, STATEVAR_ENTROPY); 661 } 662