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