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 int MassFunction_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 55 const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 56 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 57 const CeedScalar(*q_data) = in[2]; 58 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 59 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 60 61 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 62 63 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 64 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 65 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]}; 66 const State s = StateFromQ(context, qi, state_var); 67 const State s_dot = StateFromQ(context, qi_dot, state_var); 68 CeedScalar wdetJ, dXdx[3][3]; 69 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 70 71 // Standard mass matrix term 72 for (CeedInt f = 0; f < 5; f++) { 73 v[f][i] = wdetJ * qi_dot[f]; 74 } 75 76 // Stabilization method: none (Galerkin), SU, or SUPG 77 State grad_s[3] = {{{0.}}}; 78 CeedScalar Tau_d[3], stab[5][3], body_force[5] = {0.}, divFdiff[5] = {0.}, U_dot[5]; 79 UnpackState_U(s_dot.U, U_dot); 80 Tau_diagPrim(context, s, dXdx, context->dt, Tau_d); 81 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, divFdiff, stab); 82 83 // Stabilized mass term 84 for (CeedInt j = 0; j < 5; j++) { 85 for (CeedInt k = 0; k < 3; k++) { 86 Grad_v[k][j][i] = wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXdx[k][1] + stab[j][2] * dXdx[k][2]); 87 } 88 } 89 } 90 return 0; 91 } 92 93 CEED_QFUNCTION(MassFunction_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 94 return MassFunction_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 95 } 96 97 // ***************************************************************************** 98 // This QFunction implements the following formulation of Navier-Stokes with explicit time stepping method 99 // 100 // This is 3D compressible Navier-Stokes in conservation form with state variables of density, momentum density, and total energy density. 101 // 102 // State Variables: q = ( rho, U1, U2, U3, E ) 103 // rho - Mass Density 104 // Ui - Momentum Density, Ui = rho ui 105 // E - Total Energy Density, E = rho (cv T + (u u)/2 + g z) 106 // 107 // Navier-Stokes Equations: 108 // drho/dt + div( U ) = 0 109 // dU/dt + div( rho (u x u) + P I3 ) + rho g khat = div( Fu ) 110 // dE/dt + div( (E + P) u ) = div( Fe ) 111 // 112 // Viscous Stress: 113 // Fu = mu (grad( u ) + grad( u )^T + lambda div ( u ) I3) 114 // 115 // Thermal Stress: 116 // Fe = u Fu + k grad( T ) 117 // Equation of State 118 // P = (gamma - 1) (E - rho (u u) / 2 - rho g z) 119 // 120 // Stabilization: 121 // Tau = diag(TauC, TauM, TauM, TauM, TauE) 122 // f1 = rho sqrt(ui uj gij) 123 // gij = dXi/dX * dXi/dX 124 // TauC = Cc f1 / (8 gii) 125 // TauM = min( 1 , 1 / f1 ) 126 // TauE = TauM / (Ce cv) 127 // 128 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 129 // 130 // Constants: 131 // lambda = - 2 / 3, From Stokes hypothesis 132 // mu , Dynamic viscosity 133 // k , Thermal conductivity 134 // cv , Specific heat, constant volume 135 // cp , Specific heat, constant pressure 136 // g , Gravity 137 // gamma = cp / cv, Specific heat ratio 138 // 139 // 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 140 // gradu ) 141 // ***************************************************************************** 142 CEED_QFUNCTION(RHSFunction_Newtonian)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 143 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 144 const CeedScalar(*Grad_q) = in[1]; 145 const CeedScalar(*q_data) = in[2]; 146 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 147 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 148 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 149 150 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 151 const CeedScalar *g = context->g; 152 const CeedScalar dt = context->dt; 153 const CeedScalar P0 = context->idl_pressure; 154 155 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 156 CeedScalar U[5], wdetJ, dXdx[3][3]; 157 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 158 for (int j = 0; j < 5; j++) U[j] = q[j][i]; 159 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 160 State s = StateFromU(context, U); 161 162 State grad_s[3]; 163 StatePhysicalGradientFromReference(Q, i, context, s, STATEVAR_CONSERVATIVE, Grad_q, dXdx, grad_s); 164 165 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 166 KMStrainRate_State(grad_s, strain_rate); 167 NewtonianStress(context, strain_rate, kmstress); 168 KMUnpack(kmstress, stress); 169 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 170 171 StateConservative F_inviscid[3]; 172 FluxInviscid(context, s, F_inviscid); 173 174 // Total flux 175 CeedScalar Flux[5][3]; 176 FluxTotal(F_inviscid, stress, Fe, Flux); 177 178 for (CeedInt j = 0; j < 5; j++) { 179 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]); 180 } 181 182 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)}; 183 for (int j = 0; j < 5; j++) v[j][i] = wdetJ * body_force[j]; 184 185 if (context->idl_enable) { 186 const CeedScalar sigma = LinearRampCoefficient(context->idl_amplitude, context->idl_length, context->idl_start, x_i[0]); 187 CeedScalar damp_state[5] = {s.Y.pressure - P0, 0, 0, 0, 0}, idl_residual[5] = {0.}; 188 InternalDampingLayer(context, s, sigma, damp_state, idl_residual); 189 for (int j = 0; j < 5; j++) v[j][i] -= wdetJ * idl_residual[j]; 190 } 191 192 // -- Stabilization method: none (Galerkin), SU, or SUPG 193 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}, zeroFlux[5] = {0.}; 194 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 195 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, zeroFlux, stab); 196 197 for (CeedInt j = 0; j < 5; j++) { 198 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]); 199 } 200 } 201 return 0; 202 } 203 204 // ***************************************************************************** 205 // This QFunction implements the Navier-Stokes equations (mentioned above) with implicit time stepping method 206 // 207 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 208 // SUPG = Galerkin + grad(v) . ( Ai^T * Tau * (q_dot + Aj q,j - body force) ) 209 // (diffusive terms will be added later) 210 // ***************************************************************************** 211 CEED_QFUNCTION_HELPER int IFunction_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 212 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 213 214 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 215 const CeedScalar(*Grad_q) = in[1]; 216 const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 217 const CeedScalar(*q_data) = in[3]; 218 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 219 const CeedScalar(*divFdiff)[CEED_Q_VLA] = context->divFdiff_method != DIV_DIFF_FLUX_PROJ_NONE ? (const CeedScalar(*)[CEED_Q_VLA])in[5] : NULL; 220 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 221 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 222 CeedScalar(*jac_data) = out[2]; 223 224 const CeedScalar *g = context->g; 225 const CeedScalar dt = context->dt; 226 const CeedScalar P0 = context->idl_pressure; 227 228 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 229 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 230 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 231 const State s = StateFromQ(context, qi, state_var); 232 233 CeedScalar wdetJ, dXdx[3][3]; 234 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 235 State grad_s[3]; 236 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 237 238 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 239 KMStrainRate_State(grad_s, strain_rate); 240 NewtonianStress(context, strain_rate, kmstress); 241 KMUnpack(kmstress, stress); 242 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 243 244 StateConservative F_inviscid[3]; 245 FluxInviscid(context, s, F_inviscid); 246 247 // Total flux 248 CeedScalar Flux[5][3]; 249 FluxTotal(F_inviscid, stress, Fe, Flux); 250 251 for (CeedInt j = 0; j < 5; j++) { 252 for (CeedInt k = 0; k < 3; k++) { 253 Grad_v[k][j][i] = -wdetJ * (dXdx[k][0] * Flux[j][0] + dXdx[k][1] * Flux[j][1] + dXdx[k][2] * Flux[j][2]); 254 } 255 } 256 257 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)}; 258 259 // -- Stabilization method: none (Galerkin), SU, or SUPG 260 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}, qi_dot[5]; 261 for (int j = 0; j < 5; j++) qi_dot[j] = q_dot[j][i]; 262 State s_dot = StateFromQ_fwd(context, s, qi_dot, state_var); 263 UnpackState_U(s_dot.U, U_dot); 264 265 for (CeedInt j = 0; j < 5; j++) v[j][i] = wdetJ * (U_dot[j] - body_force[j]); 266 if (context->idl_enable) { 267 const CeedScalar sigma = LinearRampCoefficient(context->idl_amplitude, context->idl_length, context->idl_start, x_i[0]); 268 StoredValuesPack(Q, i, 14, 1, &sigma, jac_data); 269 CeedScalar damp_state[5] = {s.Y.pressure - P0, 0, 0, 0, 0}, idl_residual[5] = {0.}; 270 InternalDampingLayer(context, s, sigma, damp_state, idl_residual); 271 for (int j = 0; j < 5; j++) v[j][i] += wdetJ * idl_residual[j]; 272 } 273 274 CeedScalar divFdiff_i[5] = {0.}; 275 if (context->divFdiff_method != DIV_DIFF_FLUX_PROJ_NONE) { 276 for (int j = 1; j < 5; j++) divFdiff_i[j] = divFdiff[j - 1][i]; 277 } 278 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 279 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, divFdiff_i, stab); 280 281 for (CeedInt j = 0; j < 5; j++) { 282 for (CeedInt k = 0; k < 3; k++) { 283 Grad_v[k][j][i] += wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXdx[k][1] + stab[j][2] * dXdx[k][2]); 284 } 285 } 286 StoredValuesPack(Q, i, 0, 5, qi, jac_data); 287 StoredValuesPack(Q, i, 5, 6, kmstress, jac_data); 288 StoredValuesPack(Q, i, 11, 3, Tau_d, jac_data); 289 } 290 return 0; 291 } 292 293 CEED_QFUNCTION(IFunction_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 294 return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 295 } 296 297 CEED_QFUNCTION(IFunction_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 298 return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 299 } 300 301 CEED_QFUNCTION(IFunction_Newtonian_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 302 return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_ENTROPY); 303 } 304 305 // ***************************************************************************** 306 // This QFunction implements the jacobian of the Navier-Stokes equations for implicit time stepping method. 307 // ***************************************************************************** 308 CEED_QFUNCTION_HELPER int IJacobian_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 309 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 310 const CeedScalar(*Grad_dq) = in[1]; 311 const CeedScalar(*q_data) = in[2]; 312 const CeedScalar(*jac_data) = in[3]; 313 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 314 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 315 316 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 317 const CeedScalar *g = context->g; 318 319 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 320 CeedScalar wdetJ, dXdx[3][3]; 321 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 322 323 CeedScalar qi[5], kmstress[6], Tau_d[3]; 324 StoredValuesUnpack(Q, i, 0, 5, jac_data, qi); 325 StoredValuesUnpack(Q, i, 5, 6, jac_data, kmstress); 326 StoredValuesUnpack(Q, i, 11, 3, jac_data, Tau_d); 327 State s = StateFromQ(context, qi, state_var); 328 329 CeedScalar dqi[5]; 330 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 331 State ds = StateFromQ_fwd(context, s, dqi, state_var); 332 333 State grad_ds[3]; 334 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_dq, dXdx, grad_ds); 335 336 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 337 KMStrainRate_State(grad_ds, dstrain_rate); 338 NewtonianStress(context, dstrain_rate, dkmstress); 339 KMUnpack(dkmstress, dstress); 340 KMUnpack(kmstress, stress); 341 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 342 343 StateConservative dF_inviscid[3]; 344 FluxInviscid_fwd(context, s, ds, dF_inviscid); 345 346 // Total flux 347 CeedScalar dFlux[5][3]; 348 FluxTotal(dF_inviscid, dstress, dFe, dFlux); 349 350 for (int j = 0; j < 5; j++) { 351 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]); 352 } 353 354 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)}; 355 CeedScalar dU[5] = {0.}; 356 UnpackState_U(ds.U, dU); 357 for (int j = 0; j < 5; j++) v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]); 358 359 if (context->idl_enable) { 360 const CeedScalar sigma = jac_data[14 * Q + i]; 361 CeedScalar damp_state[5] = {ds.Y.pressure, 0, 0, 0, 0}, idl_residual[5] = {0.}; 362 // This is a Picard-type linearization of the damping and could be replaced by an InternalDampingLayer_fwd that uses s and ds. 363 InternalDampingLayer(context, s, sigma, damp_state, idl_residual); 364 for (int j = 0; j < 5; j++) v[j][i] += wdetJ * idl_residual[j]; 365 } 366 367 // -- Stabilization method: none (Galerkin), SU, or SUPG 368 CeedScalar dstab[5][3], U_dot[5] = {0}; 369 for (CeedInt j = 0; j < 5; j++) U_dot[j] = context->ijacobian_time_shift * dU[j]; 370 const CeedScalar zeroFlux[5] = {0.}; 371 Stabilization(context, s, Tau_d, grad_ds, U_dot, dbody_force, zeroFlux, dstab); 372 373 for (int j = 0; j < 5; j++) { 374 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]); 375 } 376 } 377 return 0; 378 } 379 380 CEED_QFUNCTION(IJacobian_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 381 return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 382 } 383 384 CEED_QFUNCTION(IJacobian_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 385 return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 386 } 387 388 CEED_QFUNCTION(IJacobian_Newtonian_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 389 return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_ENTROPY); 390 } 391 392 // ***************************************************************************** 393 // Compute boundary integral (ie. for strongly set inflows) 394 // ***************************************************************************** 395 CEED_QFUNCTION_HELPER int BoundaryIntegral(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 396 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 397 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 398 const CeedScalar(*Grad_q) = in[1]; 399 const CeedScalar(*q_data_sur) = in[2]; 400 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 401 CeedScalar(*jac_data_sur) = context->is_implicit ? out[1] : NULL; 402 403 const bool is_implicit = context->is_implicit; 404 405 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 406 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 407 State s = StateFromQ(context, qi, state_var); 408 409 CeedScalar wdetJb, dXdx[2][3], normal[3]; 410 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, normal); 411 wdetJb *= is_implicit ? -1. : 1.; 412 413 State grad_s[3]; 414 StatePhysicalGradientFromReference_Boundary(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 415 416 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 417 KMStrainRate_State(grad_s, strain_rate); 418 NewtonianStress(context, strain_rate, kmstress); 419 KMUnpack(kmstress, stress); 420 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 421 422 StateConservative F_inviscid[3]; 423 FluxInviscid(context, s, F_inviscid); 424 425 CeedScalar Flux[5]; 426 FluxTotal_Boundary(F_inviscid, stress, Fe, normal, Flux); 427 428 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j]; 429 430 if (is_implicit) { 431 StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur); 432 StoredValuesPack(Q, i, 5, 6, kmstress, jac_data_sur); 433 } 434 } 435 return 0; 436 } 437 438 CEED_QFUNCTION(BoundaryIntegral_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 439 return BoundaryIntegral(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 440 } 441 442 CEED_QFUNCTION(BoundaryIntegral_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 443 return BoundaryIntegral(ctx, Q, in, out, STATEVAR_PRIMITIVE); 444 } 445 446 CEED_QFUNCTION(BoundaryIntegral_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 447 return BoundaryIntegral(ctx, Q, in, out, STATEVAR_ENTROPY); 448 } 449 450 // ***************************************************************************** 451 // Jacobian for "set nothing" boundary integral 452 // ***************************************************************************** 453 CEED_QFUNCTION_HELPER int BoundaryIntegral_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 454 StateVariable state_var) { 455 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 456 const CeedScalar(*Grad_dq) = in[1]; 457 const CeedScalar(*q_data_sur) = in[2]; 458 const CeedScalar(*jac_data_sur) = in[4]; 459 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 460 461 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 462 const bool is_implicit = context->is_implicit; 463 464 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 465 CeedScalar wdetJb, dXdx[2][3], normal[3]; 466 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, normal); 467 wdetJb *= is_implicit ? -1. : 1.; 468 469 CeedScalar qi[5], kmstress[6], dqi[5]; 470 StoredValuesUnpack(Q, i, 0, 5, jac_data_sur, qi); 471 StoredValuesUnpack(Q, i, 5, 6, jac_data_sur, kmstress); 472 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 473 474 State s = StateFromQ(context, qi, state_var); 475 State ds = StateFromQ_fwd(context, s, dqi, state_var); 476 477 State grad_ds[3]; 478 StatePhysicalGradientFromReference_Boundary(Q, i, context, s, state_var, Grad_dq, dXdx, grad_ds); 479 480 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 481 KMStrainRate_State(grad_ds, dstrain_rate); 482 NewtonianStress(context, dstrain_rate, dkmstress); 483 KMUnpack(dkmstress, dstress); 484 KMUnpack(kmstress, stress); 485 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 486 487 StateConservative dF_inviscid[3]; 488 FluxInviscid_fwd(context, s, ds, dF_inviscid); 489 490 CeedScalar dFlux[5]; 491 FluxTotal_Boundary(dF_inviscid, dstress, dFe, normal, dFlux); 492 493 for (int j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j]; 494 } 495 return 0; 496 } 497 498 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 499 return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 500 } 501 502 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 503 return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 504 } 505 506 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 507 return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_ENTROPY); 508 } 509 510 // @brief Volume integral for RHS of divergence of diffusive flux direct projection 511 CEED_QFUNCTION_HELPER int DivDiffusiveFluxVolumeRHS_NS(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 512 StateVariable state_var) { 513 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 514 const CeedScalar(*Grad_q) = in[1]; 515 const CeedScalar(*q_data) = in[2]; 516 CeedScalar(*Grad_v)[4][CEED_Q_VLA] = (CeedScalar(*)[4][CEED_Q_VLA])out[0]; 517 518 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 519 const StateConservative ZeroInviscidFluxes[3] = {{0}}; 520 521 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 522 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 523 const State s = StateFromQ(context, qi, state_var); 524 CeedScalar wdetJ, dXdx[3][3]; 525 CeedScalar stress[3][3], Fe[3], Fdiff[5][3]; 526 527 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 528 { // Get stress and Fe 529 State grad_s[3]; 530 CeedScalar strain_rate[6], kmstress[6]; 531 532 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 533 KMStrainRate_State(grad_s, strain_rate); 534 NewtonianStress(context, strain_rate, kmstress); 535 KMUnpack(kmstress, stress); 536 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 537 } 538 539 FluxTotal(ZeroInviscidFluxes, stress, Fe, Fdiff); 540 541 for (CeedInt j = 1; j < 5; j++) { // Continuity has no diffusive flux, therefore skip 542 for (CeedInt k = 0; k < 3; k++) { 543 Grad_v[k][j - 1][i] = -wdetJ * Dot3(dXdx[k], Fdiff[j]); 544 } 545 } 546 } 547 return 0; 548 } 549 550 CEED_QFUNCTION(DivDiffusiveFluxVolumeRHS_NS_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 551 return DivDiffusiveFluxVolumeRHS_NS(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 552 } 553 554 CEED_QFUNCTION(DivDiffusiveFluxVolumeRHS_NS_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 555 return DivDiffusiveFluxVolumeRHS_NS(ctx, Q, in, out, STATEVAR_PRIMITIVE); 556 } 557 558 CEED_QFUNCTION(DivDiffusiveFluxVolumeRHS_NS_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 559 return DivDiffusiveFluxVolumeRHS_NS(ctx, Q, in, out, STATEVAR_ENTROPY); 560 } 561 562 // @brief Boundary integral for RHS of divergence of diffusive flux direct projection 563 CEED_QFUNCTION_HELPER int DivDiffusiveFluxBoundaryRHS_NS(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 564 StateVariable state_var) { 565 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 566 const CeedScalar(*Grad_q) = in[1]; 567 const CeedScalar(*q_data) = in[2]; 568 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 569 570 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 571 const StateConservative ZeroInviscidFluxes[3] = {{0}}; 572 573 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 574 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 575 const State s = StateFromQ(context, qi, state_var); 576 CeedScalar wdetJ, dXdx[3][3], normal[3]; 577 CeedScalar stress[3][3], Fe[3], Fdiff[5]; 578 579 QdataBoundaryGradientUnpack_3D(Q, i, q_data, &wdetJ, dXdx, normal); 580 { // Get stress and Fe 581 State grad_s[3]; 582 CeedScalar strain_rate[6], kmstress[6]; 583 584 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 585 KMStrainRate_State(grad_s, strain_rate); 586 NewtonianStress(context, strain_rate, kmstress); 587 KMUnpack(kmstress, stress); 588 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 589 } 590 591 FluxTotal_Boundary(ZeroInviscidFluxes, stress, Fe, normal, Fdiff); 592 593 // Continuity has no diffusive flux, therefore skip 594 for (CeedInt j = 1; j < 5; j++) v[j - 1][i] = wdetJ * Fdiff[j]; 595 } 596 return 0; 597 } 598 599 CEED_QFUNCTION(DivDiffusiveFluxBoundaryRHS_NS_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 600 return DivDiffusiveFluxBoundaryRHS_NS(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 601 } 602 603 CEED_QFUNCTION(DivDiffusiveFluxBoundaryRHS_NS_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 604 return DivDiffusiveFluxBoundaryRHS_NS(ctx, Q, in, out, STATEVAR_PRIMITIVE); 605 } 606 607 CEED_QFUNCTION(DivDiffusiveFluxBoundaryRHS_NS_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 608 return DivDiffusiveFluxBoundaryRHS_NS(ctx, Q, in, out, STATEVAR_ENTROPY); 609 } 610 611 // @brief Integral for RHS of diffusive flux indirect projection 612 CEED_QFUNCTION_HELPER int DiffusiveFluxRHS_NS(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 613 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 614 const CeedScalar(*Grad_q) = in[1]; 615 const CeedScalar(*q_data) = in[2]; 616 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 617 618 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 619 const StateConservative ZeroInviscidFluxes[3] = {{0}}; 620 621 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 622 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 623 const State s = StateFromQ(context, qi, state_var); 624 CeedScalar wdetJ, dXdx[3][3]; 625 CeedScalar stress[3][3], Fe[3], Fdiff[5][3]; 626 627 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 628 { // Get stress and Fe 629 State grad_s[3]; 630 CeedScalar strain_rate[6], kmstress[6]; 631 632 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 633 KMStrainRate_State(grad_s, strain_rate); 634 NewtonianStress(context, strain_rate, kmstress); 635 KMUnpack(kmstress, stress); 636 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 637 } 638 639 FluxTotal(ZeroInviscidFluxes, stress, Fe, Fdiff); 640 641 for (CeedInt j = 1; j < 5; j++) { // Continuity has no diffusive flux, therefore skip 642 for (CeedInt k = 0; k < 3; k++) { 643 v[(j - 1) * 3 + k][i] = wdetJ * Fdiff[j][k]; 644 } 645 } 646 } 647 return 0; 648 } 649 650 CEED_QFUNCTION(DiffusiveFluxRHS_NS_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 651 return DiffusiveFluxRHS_NS(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 652 } 653 654 CEED_QFUNCTION(DiffusiveFluxRHS_NS_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 655 return DiffusiveFluxRHS_NS(ctx, Q, in, out, STATEVAR_PRIMITIVE); 656 } 657 658 CEED_QFUNCTION(DiffusiveFluxRHS_NS_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 659 return DiffusiveFluxRHS_NS(ctx, Q, in, out, STATEVAR_ENTROPY); 660 } 661