1 // Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 /// @file 9 /// Operator for Navier-Stokes example using PETSc 10 #include <ceed.h> 11 #include <math.h> 12 #include <stdlib.h> 13 14 #include "newtonian_state.h" 15 #include "newtonian_types.h" 16 #include "stabilization.h" 17 #include "utils.h" 18 19 CEED_QFUNCTION_HELPER void InternalDampingLayer(const NewtonianIdealGasContext context, const State s, const CeedScalar sigma, CeedScalar damp_Y[5], 20 CeedScalar damp_residual[5]) { 21 ScaleN(damp_Y, sigma, 5); 22 State damp_s = StateFromY_fwd(context, s, damp_Y); 23 24 CeedScalar U[5]; 25 UnpackState_U(damp_s.U, U); 26 for (int i = 0; i < 5; i++) damp_residual[i] += U[i]; 27 } 28 29 // ***************************************************************************** 30 // This QFunction sets a "still" initial condition for generic Newtonian IG problems 31 // ***************************************************************************** 32 CEED_QFUNCTION_HELPER int ICsNewtonianIG(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 33 CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 34 35 const SetupContext context = (SetupContext)ctx; 36 37 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 38 CeedScalar q[5] = {0.}; 39 State s = StateFromPrimitive(&context->gas, context->reference); 40 StateToQ(&context->gas, s, q, state_var); 41 for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j]; 42 } 43 return 0; 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 CEED_QFUNCTION(ICsNewtonianIG_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 50 return ICsNewtonianIG(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 51 } 52 53 CEED_QFUNCTION_HELPER void MassFunction_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 54 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.}, 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, 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 } 91 92 CEED_QFUNCTION(MassFunction_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 93 MassFunction_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 94 return 0; 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 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 147 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 148 149 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 150 const CeedScalar *g = context->g; 151 const CeedScalar dt = context->dt; 152 153 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 154 CeedScalar U[5], wdetJ, dXdx[3][3]; 155 for (int j = 0; j < 5; j++) U[j] = q[j][i]; 156 StoredValuesUnpack(Q, i, 0, 1, q_data, &wdetJ); 157 StoredValuesUnpack(Q, i, 1, 9, q_data, (CeedScalar *)dXdx); 158 State s = StateFromU(context, U); 159 160 State grad_s[3]; 161 StatePhysicalGradientFromReference(Q, i, context, s, STATEVAR_CONSERVATIVE, Grad_q, dXdx, grad_s); 162 163 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 164 KMStrainRate_State(grad_s, strain_rate); 165 NewtonianStress(context, strain_rate, kmstress); 166 KMUnpack(kmstress, stress); 167 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 168 169 StateConservative F_inviscid[3]; 170 FluxInviscid(context, s, F_inviscid); 171 172 // Total flux 173 CeedScalar Flux[5][3]; 174 FluxTotal(F_inviscid, stress, Fe, Flux); 175 176 for (CeedInt j = 0; j < 5; j++) { 177 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]); 178 } 179 180 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)}; 181 for (int j = 0; j < 5; j++) v[j][i] = wdetJ * body_force[j]; 182 183 // -- Stabilization method: none (Galerkin), SU, or SUPG 184 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}; 185 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 186 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, stab); 187 188 for (CeedInt j = 0; j < 5; j++) { 189 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]); 190 } 191 } 192 return 0; 193 } 194 195 // ***************************************************************************** 196 // This QFunction implements the Navier-Stokes equations (mentioned above) with implicit time stepping method 197 // 198 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 199 // SUPG = Galerkin + grad(v) . ( Ai^T * Tau * (q_dot + Aj q,j - body force) ) 200 // (diffusive terms will be added later) 201 // ***************************************************************************** 202 CEED_QFUNCTION_HELPER int IFunction_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 203 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 204 const CeedScalar(*Grad_q) = in[1]; 205 const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 206 const CeedScalar(*q_data) = in[3]; 207 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 208 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 209 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 210 CeedScalar(*jac_data) = out[2]; 211 212 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 213 const CeedScalar *g = context->g; 214 const CeedScalar dt = context->dt; 215 const CeedScalar P0 = context->P0; 216 217 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 218 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 219 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 220 const State s = StateFromQ(context, qi, state_var); 221 222 CeedScalar wdetJ, dXdx[3][3]; 223 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 224 State grad_s[3]; 225 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 226 227 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 228 KMStrainRate_State(grad_s, strain_rate); 229 NewtonianStress(context, strain_rate, kmstress); 230 KMUnpack(kmstress, stress); 231 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 232 233 StateConservative F_inviscid[3]; 234 FluxInviscid(context, s, F_inviscid); 235 236 // Total flux 237 CeedScalar Flux[5][3]; 238 FluxTotal(F_inviscid, stress, Fe, Flux); 239 240 for (CeedInt j = 0; j < 5; j++) { 241 for (CeedInt k = 0; k < 3; k++) { 242 Grad_v[k][j][i] = -wdetJ * (dXdx[k][0] * Flux[j][0] + dXdx[k][1] * Flux[j][1] + dXdx[k][2] * Flux[j][2]); 243 } 244 } 245 246 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)}; 247 248 // -- Stabilization method: none (Galerkin), SU, or SUPG 249 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}, qi_dot[5]; 250 for (int j = 0; j < 5; j++) qi_dot[j] = q_dot[j][i]; 251 State s_dot = StateFromQ_fwd(context, s, qi_dot, state_var); 252 UnpackState_U(s_dot.U, U_dot); 253 254 for (CeedInt j = 0; j < 5; j++) v[j][i] = wdetJ * (U_dot[j] - body_force[j]); 255 if (context->idl_enable) { 256 const CeedScalar sigma = LinearRampCoefficient(context->idl_amplitude, context->idl_length, context->idl_start, x_i[0]); 257 StoredValuesPack(Q, i, 14, 1, &sigma, jac_data); 258 CeedScalar damp_state[5] = {s.Y.pressure - P0, 0, 0, 0, 0}, idl_residual[5] = {0.}; 259 InternalDampingLayer(context, s, sigma, damp_state, idl_residual); 260 for (int j = 0; j < 5; j++) v[j][i] += wdetJ * idl_residual[j]; 261 } 262 263 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 264 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, stab); 265 266 for (CeedInt j = 0; j < 5; j++) { 267 for (CeedInt k = 0; k < 3; k++) { 268 Grad_v[k][j][i] += wdetJ * (stab[j][0] * dXdx[k][0] + stab[j][1] * dXdx[k][1] + stab[j][2] * dXdx[k][2]); 269 } 270 } 271 StoredValuesPack(Q, i, 0, 5, qi, jac_data); 272 StoredValuesPack(Q, i, 5, 6, kmstress, jac_data); 273 StoredValuesPack(Q, i, 11, 3, Tau_d, jac_data); 274 } 275 return 0; 276 } 277 278 CEED_QFUNCTION(IFunction_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 279 return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 280 } 281 282 CEED_QFUNCTION(IFunction_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 283 return IFunction_Newtonian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 284 } 285 286 // ***************************************************************************** 287 // This QFunction implements the jacobian of the Navier-Stokes equations for implicit time stepping method. 288 // ***************************************************************************** 289 CEED_QFUNCTION_HELPER int IJacobian_Newtonian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 290 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 291 const CeedScalar(*Grad_dq) = in[1]; 292 const CeedScalar(*q_data) = in[2]; 293 const CeedScalar(*jac_data) = in[3]; 294 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 295 CeedScalar(*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 296 297 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 298 const CeedScalar *g = context->g; 299 300 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 301 CeedScalar wdetJ, dXdx[3][3]; 302 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 303 304 CeedScalar qi[5], kmstress[6], Tau_d[3]; 305 StoredValuesUnpack(Q, i, 0, 5, jac_data, qi); 306 StoredValuesUnpack(Q, i, 5, 6, jac_data, kmstress); 307 StoredValuesUnpack(Q, i, 11, 3, jac_data, Tau_d); 308 State s = StateFromQ(context, qi, state_var); 309 310 CeedScalar dqi[5]; 311 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 312 State ds = StateFromQ_fwd(context, s, dqi, state_var); 313 314 State grad_ds[3]; 315 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_dq, dXdx, grad_ds); 316 317 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 318 KMStrainRate_State(grad_ds, dstrain_rate); 319 NewtonianStress(context, dstrain_rate, dkmstress); 320 KMUnpack(dkmstress, dstress); 321 KMUnpack(kmstress, stress); 322 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 323 324 StateConservative dF_inviscid[3]; 325 FluxInviscid_fwd(context, s, ds, dF_inviscid); 326 327 // Total flux 328 CeedScalar dFlux[5][3]; 329 FluxTotal(dF_inviscid, dstress, dFe, dFlux); 330 331 for (int j = 0; j < 5; j++) { 332 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]); 333 } 334 335 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)}; 336 CeedScalar dU[5] = {0.}; 337 UnpackState_U(ds.U, dU); 338 for (int j = 0; j < 5; j++) v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]); 339 340 if (context->idl_enable) { 341 const CeedScalar sigma = jac_data[14 * Q + i]; 342 CeedScalar damp_state[5] = {ds.Y.pressure, 0, 0, 0, 0}, idl_residual[5] = {0.}; 343 // This is a Picard-type linearization of the damping and could be replaced by an InternalDampingLayer_fwd that uses s and ds. 344 InternalDampingLayer(context, s, sigma, damp_state, idl_residual); 345 for (int j = 0; j < 5; j++) v[j][i] += wdetJ * idl_residual[j]; 346 } 347 348 // -- Stabilization method: none (Galerkin), SU, or SUPG 349 CeedScalar dstab[5][3], U_dot[5] = {0}; 350 for (CeedInt j = 0; j < 5; j++) U_dot[j] = context->ijacobian_time_shift * dU[j]; 351 Stabilization(context, s, Tau_d, grad_ds, U_dot, dbody_force, dstab); 352 353 for (int j = 0; j < 5; j++) { 354 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]); 355 } 356 } 357 return 0; 358 } 359 360 CEED_QFUNCTION(IJacobian_Newtonian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 361 return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 362 } 363 364 CEED_QFUNCTION(IJacobian_Newtonian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 365 return IJacobian_Newtonian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 366 } 367 368 // ***************************************************************************** 369 // Compute boundary integral (ie. for strongly set inflows) 370 // ***************************************************************************** 371 CEED_QFUNCTION_HELPER int BoundaryIntegral(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 372 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 373 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 374 const CeedScalar(*Grad_q) = in[1]; 375 const CeedScalar(*q_data_sur) = in[2]; 376 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 377 CeedScalar(*jac_data_sur) = context->is_implicit ? out[1] : NULL; 378 379 const bool is_implicit = context->is_implicit; 380 381 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 382 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 383 State s = StateFromQ(context, qi, state_var); 384 385 CeedScalar wdetJb, dXdx[2][3], norm[3]; 386 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, norm); 387 wdetJb *= is_implicit ? -1. : 1.; 388 389 State grad_s[3]; 390 StatePhysicalGradientFromReference_Boundary(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 391 392 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 393 KMStrainRate_State(grad_s, strain_rate); 394 NewtonianStress(context, strain_rate, kmstress); 395 KMUnpack(kmstress, stress); 396 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 397 398 StateConservative F_inviscid[3]; 399 FluxInviscid(context, s, F_inviscid); 400 401 CeedScalar Flux[5]; 402 FluxTotal_Boundary(F_inviscid, stress, Fe, norm, Flux); 403 404 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j]; 405 406 if (is_implicit) { 407 StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur); 408 StoredValuesPack(Q, i, 5, 6, kmstress, jac_data_sur); 409 } 410 } 411 return 0; 412 } 413 414 CEED_QFUNCTION(BoundaryIntegral_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 415 return BoundaryIntegral(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 416 } 417 418 CEED_QFUNCTION(BoundaryIntegral_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 419 return BoundaryIntegral(ctx, Q, in, out, STATEVAR_PRIMITIVE); 420 } 421 422 // ***************************************************************************** 423 // Jacobian for "set nothing" boundary integral 424 // ***************************************************************************** 425 CEED_QFUNCTION_HELPER int BoundaryIntegral_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 426 StateVariable state_var) { 427 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 428 const CeedScalar(*Grad_dq) = in[1]; 429 const CeedScalar(*q_data_sur) = in[2]; 430 const CeedScalar(*jac_data_sur) = in[4]; 431 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 432 433 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 434 const bool is_implicit = context->is_implicit; 435 436 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 437 CeedScalar wdetJb, dXdx[2][3], norm[3]; 438 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, norm); 439 wdetJb *= is_implicit ? -1. : 1.; 440 441 CeedScalar qi[5], kmstress[6], dqi[5]; 442 StoredValuesUnpack(Q, i, 0, 5, jac_data_sur, qi); 443 StoredValuesUnpack(Q, i, 5, 6, jac_data_sur, kmstress); 444 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 445 446 State s = StateFromQ(context, qi, state_var); 447 State ds = StateFromQ_fwd(context, s, dqi, state_var); 448 449 State grad_ds[3]; 450 StatePhysicalGradientFromReference_Boundary(Q, i, context, s, state_var, Grad_dq, dXdx, grad_ds); 451 452 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 453 KMStrainRate_State(grad_ds, dstrain_rate); 454 NewtonianStress(context, dstrain_rate, dkmstress); 455 KMUnpack(dkmstress, dstress); 456 KMUnpack(kmstress, stress); 457 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 458 459 StateConservative dF_inviscid[3]; 460 FluxInviscid_fwd(context, s, ds, dF_inviscid); 461 462 CeedScalar dFlux[5]; 463 FluxTotal_Boundary(dF_inviscid, dstress, dFe, norm, dFlux); 464 465 for (int j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j]; 466 } 467 return 0; 468 } 469 470 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 471 return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 472 } 473 474 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 475 return BoundaryIntegral_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 476 } 477