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 /// QFunctions for the `bc_freestream` and `bc_outflow` boundary conditions 10 #include "bc_freestream_type.h" 11 #include "newtonian_state.h" 12 #include "newtonian_types.h" 13 #include "riemann_solver.h" 14 15 // ***************************************************************************** 16 // Freestream Boundary Condition 17 // ***************************************************************************** 18 CEED_QFUNCTION_HELPER int Freestream(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var, 19 RiemannFluxType flux_type) { 20 const FreestreamContext context = (FreestreamContext)ctx; 21 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 22 const CeedScalar(*q_data_sur) = in[2]; 23 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 24 CeedScalar(*jac_data_sur) = context->newtonian_ctx.is_implicit ? out[1] : NULL; 25 26 const NewtonianIdealGasContext newt_ctx = &context->newtonian_ctx; 27 const bool is_implicit = newt_ctx->is_implicit; 28 29 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 30 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 31 const State s = StateFromQ(newt_ctx, qi, state_var); 32 33 CeedScalar wdetJb, normal[3]; 34 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, normal); 35 wdetJb *= is_implicit ? -1. : 1.; 36 37 StateConservative flux; 38 switch (flux_type) { 39 case RIEMANN_HLL: 40 flux = RiemannFlux_HLL(newt_ctx, s, context->S_infty, normal); 41 break; 42 case RIEMANN_HLLC: 43 flux = RiemannFlux_HLLC(newt_ctx, s, context->S_infty, normal); 44 break; 45 } 46 CeedScalar Flux[5]; 47 UnpackState_U(flux, Flux); 48 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j]; 49 50 if (is_implicit) { 51 CeedScalar zeros[6] = {0.}; 52 StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur); 53 StoredValuesPack(Q, i, 5, 6, zeros, jac_data_sur); // Every output value must be set 54 } 55 } 56 return 0; 57 } 58 59 CEED_QFUNCTION(Freestream_Conserv_HLL)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 60 return Freestream(ctx, Q, in, out, STATEVAR_CONSERVATIVE, RIEMANN_HLL); 61 } 62 63 CEED_QFUNCTION(Freestream_Prim_HLL)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 64 return Freestream(ctx, Q, in, out, STATEVAR_PRIMITIVE, RIEMANN_HLL); 65 } 66 67 CEED_QFUNCTION(Freestream_Entropy_HLL)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 68 return Freestream(ctx, Q, in, out, STATEVAR_ENTROPY, RIEMANN_HLL); 69 } 70 71 CEED_QFUNCTION(Freestream_Conserv_HLLC)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 72 return Freestream(ctx, Q, in, out, STATEVAR_CONSERVATIVE, RIEMANN_HLLC); 73 } 74 75 CEED_QFUNCTION(Freestream_Prim_HLLC)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 76 return Freestream(ctx, Q, in, out, STATEVAR_PRIMITIVE, RIEMANN_HLLC); 77 } 78 79 CEED_QFUNCTION(Freestream_Entropy_HLLC)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 80 return Freestream(ctx, Q, in, out, STATEVAR_ENTROPY, RIEMANN_HLLC); 81 } 82 83 CEED_QFUNCTION_HELPER int Freestream_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var, 84 RiemannFluxType flux_type) { 85 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 86 const CeedScalar(*q_data_sur) = in[2]; 87 const CeedScalar(*jac_data_sur) = in[4]; 88 89 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 90 91 const FreestreamContext context = (FreestreamContext)ctx; 92 const NewtonianIdealGasContext newt_ctx = &context->newtonian_ctx; 93 const bool is_implicit = newt_ctx->is_implicit; 94 const State dS_infty = {0}; 95 96 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 97 CeedScalar wdetJb, normal[3]; 98 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, normal); 99 wdetJb *= is_implicit ? -1. : 1.; 100 101 CeedScalar qi[5], dqi[5]; 102 StoredValuesUnpack(Q, i, 0, 5, jac_data_sur, qi); 103 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 104 State s = StateFromQ(newt_ctx, qi, state_var); 105 State ds = StateFromQ_fwd(newt_ctx, s, dqi, state_var); 106 107 StateConservative dflux; 108 switch (flux_type) { 109 case RIEMANN_HLL: 110 dflux = RiemannFlux_HLL_fwd(newt_ctx, s, ds, context->S_infty, dS_infty, normal); 111 break; 112 case RIEMANN_HLLC: 113 dflux = RiemannFlux_HLLC_fwd(newt_ctx, s, ds, context->S_infty, dS_infty, normal); 114 break; 115 } 116 CeedScalar dFlux[5]; 117 UnpackState_U(dflux, dFlux); 118 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j]; 119 } 120 return 0; 121 } 122 123 CEED_QFUNCTION(Freestream_Jacobian_Conserv_HLL)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 124 return Freestream_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE, RIEMANN_HLL); 125 } 126 127 CEED_QFUNCTION(Freestream_Jacobian_Prim_HLL)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 128 return Freestream_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE, RIEMANN_HLL); 129 } 130 131 CEED_QFUNCTION(Freestream_Jacobian_Entropy_HLL)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 132 return Freestream_Jacobian(ctx, Q, in, out, STATEVAR_ENTROPY, RIEMANN_HLL); 133 } 134 135 CEED_QFUNCTION(Freestream_Jacobian_Conserv_HLLC)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 136 return Freestream_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE, RIEMANN_HLLC); 137 } 138 139 CEED_QFUNCTION(Freestream_Jacobian_Prim_HLLC)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 140 return Freestream_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE, RIEMANN_HLLC); 141 } 142 143 CEED_QFUNCTION(Freestream_Jacobian_Entropy_HLLC)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 144 return Freestream_Jacobian(ctx, Q, in, out, STATEVAR_ENTROPY, RIEMANN_HLLC); 145 } 146 147 // Note the identity 148 // 149 // softplus(x) - x = log(1 + exp(x)) - x 150 // = log(1 + exp(x)) + log(exp(-x)) 151 // = log((1 + exp(x)) * exp(-x)) 152 // = log(exp(-x) + 1) 153 // = softplus(-x) 154 CEED_QFUNCTION_HELPER CeedScalar Softplus(CeedScalar x, CeedScalar width) { 155 if (x > 40 * width) return x; 156 return width * log1p(exp(x / width)); 157 } 158 159 CEED_QFUNCTION_HELPER CeedScalar Softplus_fwd(CeedScalar x, CeedScalar dx, CeedScalar width) { 160 if (x > 40 * width) return 1; 161 const CeedScalar t = exp(x / width); 162 return t / (1 + t); 163 } 164 165 // Viscous Outflow boundary condition, setting a constant exterior pressure and 166 // temperature as input for a Riemann solve. This condition is stable even in 167 // recirculating flow so long as the exterior temperature is sensible. 168 // 169 // The velocity in the exterior state has optional softplus regularization to 170 // keep it outflow. These parameters have been finnicky in practice and provide 171 // little or no benefit in the tests we've run thus far, thus we recommend 172 // skipping this feature and just allowing recirculation. 173 CEED_QFUNCTION_HELPER int RiemannOutflow(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 174 const OutflowContext outflow = (OutflowContext)ctx; 175 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 176 const CeedScalar(*Grad_q) = in[1]; 177 const CeedScalar(*q_data_sur) = in[2]; 178 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 179 CeedScalar(*jac_data_sur) = outflow->gas.is_implicit ? out[1] : NULL; 180 181 const NewtonianIdealGasContext gas = &outflow->gas; 182 const bool is_implicit = gas->is_implicit; 183 184 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 185 CeedScalar wdetJb, dXdx[2][3], normal[3]; 186 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, normal); 187 wdetJb *= is_implicit ? -1. : 1.; 188 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 189 const State s_int = StateFromQ(gas, qi, state_var); 190 191 StatePrimitive y_ext = s_int.Y; 192 y_ext.pressure = outflow->pressure; 193 y_ext.temperature = outflow->temperature; 194 const CeedScalar u_normal = Dot3(y_ext.velocity, normal); 195 const CeedScalar proj = (1 - outflow->recirc) * Softplus(-u_normal, outflow->softplus_velocity); 196 for (CeedInt j = 0; j < 3; j++) { 197 y_ext.velocity[j] += normal[j] * proj; // (I - n n^T) projects into the plane tangent to the normal 198 } 199 State s_ext = StateFromPrimitive(gas, y_ext); 200 201 State grad_s[3]; 202 StatePhysicalGradientFromReference_Boundary(Q, i, gas, s_int, state_var, Grad_q, dXdx, grad_s); 203 204 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 205 KMStrainRate_State(grad_s, strain_rate); 206 NewtonianStress(gas, strain_rate, kmstress); 207 KMUnpack(kmstress, stress); 208 ViscousEnergyFlux(gas, s_int.Y, grad_s, stress, Fe); 209 210 StateConservative F_inviscid_normal = RiemannFlux_HLLC(gas, s_int, s_ext, normal); 211 212 CeedScalar Flux[5]; 213 FluxTotal_RiemannBoundary(F_inviscid_normal, stress, Fe, normal, Flux); 214 215 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j]; 216 217 // Save values for Jacobian 218 if (is_implicit) { 219 StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur); 220 StoredValuesPack(Q, i, 5, 6, kmstress, jac_data_sur); 221 } 222 } 223 return 0; 224 } 225 226 CEED_QFUNCTION(RiemannOutflow_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 227 return RiemannOutflow(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 228 } 229 230 CEED_QFUNCTION(RiemannOutflow_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 231 return RiemannOutflow(ctx, Q, in, out, STATEVAR_PRIMITIVE); 232 } 233 234 CEED_QFUNCTION(RiemannOutflow_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 235 return RiemannOutflow(ctx, Q, in, out, STATEVAR_ENTROPY); 236 } 237 238 // ***************************************************************************** 239 // Jacobian for Riemann pressure/temperature outflow boundary condition 240 // ***************************************************************************** 241 CEED_QFUNCTION_HELPER int RiemannOutflow_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 242 StateVariable state_var) { 243 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 244 const CeedScalar(*Grad_dq) = in[1]; 245 const CeedScalar(*q_data_sur) = in[2]; 246 const CeedScalar(*jac_data_sur) = in[4]; 247 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 248 249 const OutflowContext outflow = (OutflowContext)ctx; 250 const NewtonianIdealGasContext gas = &outflow->gas; 251 const bool is_implicit = gas->is_implicit; 252 253 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 254 CeedScalar wdetJb, dXdx[2][3], normal[3]; 255 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, normal); 256 wdetJb *= is_implicit ? -1. : 1.; 257 258 CeedScalar qi[5], kmstress[6], dqi[5]; 259 StoredValuesUnpack(Q, i, 0, 5, jac_data_sur, qi); 260 StoredValuesUnpack(Q, i, 5, 6, jac_data_sur, kmstress); 261 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 262 263 State s_int = StateFromQ(gas, qi, state_var); 264 const State ds_int = StateFromQ_fwd(gas, s_int, dqi, state_var); 265 StatePrimitive y_ext = s_int.Y, dy_ext = ds_int.Y; 266 y_ext.pressure = outflow->pressure; 267 y_ext.temperature = outflow->temperature; 268 dy_ext.pressure = 0; 269 dy_ext.temperature = 0; 270 const CeedScalar u_normal = Dot3(s_int.Y.velocity, normal); 271 const CeedScalar du_normal = Dot3(ds_int.Y.velocity, normal); 272 const CeedScalar proj = (1 - outflow->recirc) * Softplus(-u_normal, outflow->softplus_velocity); 273 const CeedScalar dproj = (1 - outflow->recirc) * Softplus_fwd(-u_normal, -du_normal, outflow->softplus_velocity); 274 for (CeedInt j = 0; j < 3; j++) { 275 y_ext.velocity[j] += normal[j] * proj; 276 dy_ext.velocity[j] += normal[j] * dproj; 277 } 278 279 State s_ext = StateFromPrimitive(gas, y_ext); 280 State ds_ext = StateFromPrimitive_fwd(gas, s_ext, dy_ext); 281 282 State grad_ds[3]; 283 StatePhysicalGradientFromReference_Boundary(Q, i, gas, s_int, state_var, Grad_dq, dXdx, grad_ds); 284 285 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 286 KMStrainRate_State(grad_ds, dstrain_rate); 287 NewtonianStress(gas, dstrain_rate, dkmstress); 288 KMUnpack(dkmstress, dstress); 289 KMUnpack(kmstress, stress); 290 ViscousEnergyFlux_fwd(gas, s_int.Y, ds_int.Y, grad_ds, stress, dstress, dFe); 291 292 StateConservative dF_inviscid_normal = RiemannFlux_HLLC_fwd(gas, s_int, ds_int, s_ext, ds_ext, normal); 293 294 CeedScalar dFlux[5]; 295 FluxTotal_RiemannBoundary(dF_inviscid_normal, dstress, dFe, normal, dFlux); 296 297 for (int j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j]; 298 } 299 return 0; 300 } 301 302 CEED_QFUNCTION(RiemannOutflow_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 303 return RiemannOutflow_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 304 } 305 306 CEED_QFUNCTION(RiemannOutflow_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 307 return RiemannOutflow_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 308 } 309 310 CEED_QFUNCTION(RiemannOutflow_Jacobian_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 311 return RiemannOutflow_Jacobian(ctx, Q, in, out, STATEVAR_ENTROPY); 312 } 313 314 // ***************************************************************************** 315 // Outflow boundary condition, weakly setting a constant pressure. This is the 316 // classic outflow condition used by PHASTA-C and retained largely for 317 // comparison. In our experiments, it is never better than RiemannOutflow, and 318 // will crash if outflow ever becomes an inflow, as occurs with strong 319 // acoustics, vortices, etc. 320 // ***************************************************************************** 321 CEED_QFUNCTION_HELPER int PressureOutflow(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 322 const OutflowContext outflow = (OutflowContext)ctx; 323 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 324 const CeedScalar(*Grad_q) = in[1]; 325 const CeedScalar(*q_data_sur) = in[2]; 326 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 327 CeedScalar(*jac_data_sur) = outflow->gas.is_implicit ? out[1] : NULL; 328 329 const NewtonianIdealGasContext gas = &outflow->gas; 330 const bool is_implicit = gas->is_implicit; 331 332 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 333 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 334 State s = StateFromQ(gas, qi, state_var); 335 s.Y.pressure = outflow->pressure; 336 337 CeedScalar wdetJb, dXdx[2][3], normal[3]; 338 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, normal); 339 wdetJb *= is_implicit ? -1. : 1.; 340 341 State grad_s[3]; 342 StatePhysicalGradientFromReference_Boundary(Q, i, gas, s, state_var, Grad_q, dXdx, grad_s); 343 344 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 345 KMStrainRate_State(grad_s, strain_rate); 346 NewtonianStress(gas, strain_rate, kmstress); 347 KMUnpack(kmstress, stress); 348 ViscousEnergyFlux(gas, s.Y, grad_s, stress, Fe); 349 350 StateConservative F_inviscid[3]; 351 FluxInviscid(gas, s, F_inviscid); 352 353 CeedScalar Flux[5]; 354 FluxTotal_Boundary(F_inviscid, stress, Fe, normal, Flux); 355 356 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j]; 357 358 // Save values for Jacobian 359 if (is_implicit) { 360 StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur); 361 StoredValuesPack(Q, i, 5, 6, kmstress, jac_data_sur); 362 } 363 } 364 return 0; 365 } 366 367 CEED_QFUNCTION(PressureOutflow_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 368 return PressureOutflow(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 369 } 370 371 CEED_QFUNCTION(PressureOutflow_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 372 return PressureOutflow(ctx, Q, in, out, STATEVAR_PRIMITIVE); 373 } 374 375 CEED_QFUNCTION(PressureOutflow_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 376 return PressureOutflow(ctx, Q, in, out, STATEVAR_ENTROPY); 377 } 378 379 // ***************************************************************************** 380 // Jacobian for weak-pressure outflow boundary condition 381 // ***************************************************************************** 382 CEED_QFUNCTION_HELPER int PressureOutflow_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 383 StateVariable state_var) { 384 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 385 const CeedScalar(*Grad_dq) = in[1]; 386 const CeedScalar(*q_data_sur) = in[2]; 387 const CeedScalar(*jac_data_sur) = in[4]; 388 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 389 390 const OutflowContext outflow = (OutflowContext)ctx; 391 const NewtonianIdealGasContext gas = &outflow->gas; 392 const bool is_implicit = gas->is_implicit; 393 394 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 395 CeedScalar wdetJb, dXdx[2][3], normal[3]; 396 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, dXdx, normal); 397 wdetJb *= is_implicit ? -1. : 1.; 398 399 CeedScalar qi[5], kmstress[6], dqi[5]; 400 StoredValuesUnpack(Q, i, 0, 5, jac_data_sur, qi); 401 StoredValuesUnpack(Q, i, 5, 6, jac_data_sur, kmstress); 402 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 403 404 State s = StateFromQ(gas, qi, state_var); 405 State ds = StateFromQ_fwd(gas, s, dqi, state_var); 406 s.Y.pressure = outflow->pressure; 407 ds.Y.pressure = 0.; 408 409 State grad_ds[3]; 410 StatePhysicalGradientFromReference_Boundary(Q, i, gas, s, state_var, Grad_dq, dXdx, grad_ds); 411 412 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 413 KMStrainRate_State(grad_ds, dstrain_rate); 414 NewtonianStress(gas, dstrain_rate, dkmstress); 415 KMUnpack(dkmstress, dstress); 416 KMUnpack(kmstress, stress); 417 ViscousEnergyFlux_fwd(gas, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 418 419 StateConservative dF_inviscid[3]; 420 FluxInviscid_fwd(gas, s, ds, dF_inviscid); 421 422 CeedScalar dFlux[5]; 423 FluxTotal_Boundary(dF_inviscid, dstress, dFe, normal, dFlux); 424 425 for (int j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j]; 426 } 427 return 0; 428 } 429 430 CEED_QFUNCTION(PressureOutflow_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 431 return PressureOutflow_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 432 } 433 434 CEED_QFUNCTION(PressureOutflow_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 435 return PressureOutflow_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 436 } 437 438 CEED_QFUNCTION(PressureOutflow_Jacobian_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 439 return PressureOutflow_Jacobian(ctx, Q, in, out, STATEVAR_ENTROPY); 440 } 441