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