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