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