From dd5c2fd31c7f28e4444df4caf9eef33194050d98 Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Sun, 7 Apr 2024 20:14:51 -0400 Subject: [PATCH] add simple subsonic inflow/outflow BCs (#598) ### Description This adds simplified subsonic inflow and outflow boundary conditions. For outflows, it specifies the pressure at the outflow and extrapolates the other primitive variables (while turning into a reflecting boundary if the normal velocity is inflowing). For inflows, it specifies the velocity components, the temperature, and the passive scalars, while extrapolating the density. Unlike the naive boundary conditions, these boundary conditions are well-posed, since they specify the correct number of quantities required for the steady-state solution to be unique. The amplitude of reflected waves is greater than for the NSCBC boundary conditions, but these boundary conditions are much more stable in the presence of supersonic turbulent flow near the boundary. ### Related issues Workaround for https://github.com/quokka-astro/quokka/issues/532. ### Checklist _Before this pull request can be reviewed, all of these tasks should be completed. Denote completed tasks with an `x` inside the square brackets `[ ]` in the Markdown source below:_ - [x] I have added a description (see above). - [x] I have added a link to any related issues see (see above). - [x] I have read the [Contributing Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md). - [ ] I have added tests for any new physics that this PR adds to the code. - [x] I have tested this PR on my local computer and all tests pass. - [x] I have manually triggered the GPU tests with the magic comment `/azp run`. - [ ] I have requested a reviewer for this PR. --- src/NSCBC_inflow.hpp | 53 +++++++++++++++++++- src/NSCBC_outflow.hpp | 114 +++++++++++++++++++++++++++++++++++++----- 2 files changed, 154 insertions(+), 13 deletions(-) diff --git a/src/NSCBC_inflow.hpp b/src/NSCBC_inflow.hpp index 72e145014..823580062 100644 --- a/src/NSCBC_inflow.hpp +++ b/src/NSCBC_inflow.hpp @@ -31,6 +31,8 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto dQ_dx_inflow_x1_lower(quokka::valarray< // (This is only necessary for continuous inflows, i.e., where a shock or contact discontinuity is not desired.) // NOTE: This boundary condition is only defined for an ideal gas (with constant k_B/mu). + // TODO(bwibking): add transverse terms + const Real rho = Q[0]; const Real u = Q[1]; const Real v = Q[2]; @@ -49,7 +51,7 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto dQ_dx_inflow_x1_lower(quokka::valarray< const Real dP_dx = dQ_dx_data[4]; const Real c = quokka::EOS::ComputeSoundSpeed(rho, P); - const Real M = std::sqrt(u * u + v * v + w * w) / c; + const amrex::Real M = std::clamp(std::sqrt(u * u + v * v + w * w) / c, 0., 1.); const Real eta_2 = 2.; const Real eta_3 = 2.; @@ -145,6 +147,55 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE void setInflowX1Lower(const amrex::IntVect & consVar(i, j, k, HydroSystem::scalar0_index + n) = consCell[6 + n]; } } + +template +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void setInflowX1LowerLowOrder(const amrex::IntVect &iv, amrex::Array4 const &consVar, + amrex::GeometryData const &geom, const amrex::Real T_t, const amrex::Real u_t, + const amrex::Real v_t, const amrex::Real w_t, + amrex::GpuArray::nscalars_> const &s_t) +{ + // x1 upper boundary -- subsonic outflow + auto [i, j, k] = iv.dim3(); + amrex::Box const &box = geom.Domain(); + const auto &domain_lo = box.loVect3d(); + const int ilo = domain_lo[0]; + const Real dx = geom.CellSize(0); + const Real Lx = geom.prob_domain.length(0); + constexpr int N = HydroSystem::nvar_; + + /// x1 lower boundary -- subsonic inflow + + // compute primitive vars from valid region + quokka::valarray const Q_i = HydroSystem::ComputePrimVars(consVar, ilo, j, k); + + // compute centered ghost values + const Real rho = Q_i[0]; + const Real Eint = quokka::EOS::ComputeEintFromTgas(rho, T_t); + quokka::valarray Q_im1{}; + Q_im1[0] = rho; // extrapolate density + Q_im1[1] = u_t; // prescribe velocity + Q_im1[2] = v_t; + Q_im1[3] = w_t; + Q_im1[4] = quokka::EOS::ComputePressure(rho, Eint); // prescribe temperature + Q_im1[5] = Eint; // prescribe temperature + for (int i = 0; i < HydroSystem::nscalars_; ++i) { + Q_im1[6 + i] = s_t[i]; // prescribe passive scalars + } + + // set cell values + quokka::valarray consCell{}; + consCell = HydroSystem::ComputeConsVars(Q_im1); + + consVar(i, j, k, HydroSystem::density_index) = consCell[0]; + consVar(i, j, k, HydroSystem::x1Momentum_index) = consCell[1]; + consVar(i, j, k, HydroSystem::x2Momentum_index) = consCell[2]; + consVar(i, j, k, HydroSystem::x3Momentum_index) = consCell[3]; + consVar(i, j, k, HydroSystem::energy_index) = consCell[4]; + consVar(i, j, k, HydroSystem::internalEnergy_index) = consCell[5]; + for (int n = 0; n < HydroSystem::nscalars_; ++n) { + consVar(i, j, k, HydroSystem::scalar0_index + n) = consCell[6 + n]; + } +} } // namespace NSCBC #endif // NSCBC_INFLOW_HPP_ \ No newline at end of file diff --git a/src/NSCBC_outflow.hpp b/src/NSCBC_outflow.hpp index 35f3f1b14..55a442395 100644 --- a/src/NSCBC_outflow.hpp +++ b/src/NSCBC_outflow.hpp @@ -249,15 +249,6 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto unpermute_vel(quokka::valarray -AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto isStateValid(quokka::valarray::nvar_> const &Q) -> bool -{ - const amrex::Real rho = Q[0]; - const amrex::Real P = Q[4]; - // check whether density and pressure are positive - return ((rho > 0.) && (P > 0.)); -} } // namespace detail template @@ -334,15 +325,114 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE void setOutflowBoundary(const amrex::IntVect const int ip3 = (SIDE == BoundarySide::Lower) ? ibr - 3 : ibr + 3; const int ip4 = (SIDE == BoundarySide::Lower) ? ibr - 4 : ibr + 4; - // reset to zero-gradient outflow if any state Q_{ip1...ip4} is invalid - if (!(detail::isStateValid(Q_ip1) && detail::isStateValid(Q_ip2) && detail::isStateValid(Q_ip3) && - detail::isStateValid(Q_ip4))) { + quokka::valarray consCell{}; + if (idx[static_cast(DIR)] == ip1) { + consCell = HydroSystem::ComputeConsVars(Q_ip1); + } else if (idx[static_cast(DIR)] == ip2) { + consCell = HydroSystem::ComputeConsVars(Q_ip2); + } else if (idx[static_cast(DIR)] == ip3) { + consCell = HydroSystem::ComputeConsVars(Q_ip3); + } else if (idx[static_cast(DIR)] == ip4) { + consCell = HydroSystem::ComputeConsVars(Q_ip4); + } + + consVar(i, j, k, HydroSystem::density_index) = consCell[0]; + consVar(i, j, k, HydroSystem::x1Momentum_index) = consCell[1]; + consVar(i, j, k, HydroSystem::x2Momentum_index) = consCell[2]; + consVar(i, j, k, HydroSystem::x3Momentum_index) = consCell[3]; + consVar(i, j, k, HydroSystem::energy_index) = consCell[4]; + consVar(i, j, k, HydroSystem::internalEnergy_index) = consCell[5]; + for (int n = 0; n < HydroSystem::nscalars_; ++n) { + consVar(i, j, k, HydroSystem::scalar0_index + n) = consCell[6 + n]; + } +} + +template +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void setOutflowBoundaryLowOrder(const amrex::IntVect &iv, amrex::Array4 const &consVar, + amrex::GeometryData const &geom, const amrex::Real P_outflow) +{ + // subsonic outflow on the DIR SIDE boundary + + auto [i, j, k] = iv.dim3(); + std::array idx{i, j, k}; + amrex::Box const &box = geom.Domain(); + constexpr int N = HydroSystem::nvar_; + + const auto &boundary_idx = (SIDE == BoundarySide::Lower) ? box.loVect3d() : box.hiVect3d(); + const int ibr = boundary_idx[static_cast(DIR)]; + const int im1 = (SIDE == BoundarySide::Lower) ? ibr + 1 : ibr - 1; + const int im2 = (SIDE == BoundarySide::Lower) ? ibr + 2 : ibr - 2; + const int im3 = (SIDE == BoundarySide::Lower) ? ibr + 3 : ibr - 3; + const Real dx = geom.CellSize(static_cast(DIR)); + + // compute primitive vars + quokka::valarray Q_i{}; + if constexpr (DIR == FluxDir::X1) { + Q_i = HydroSystem::ComputePrimVars(consVar, ibr, j, k); + } else if constexpr (DIR == FluxDir::X2) { + Q_i = HydroSystem::ComputePrimVars(consVar, i, ibr, k); + } else if constexpr (DIR == FluxDir::X3) { + Q_i = HydroSystem::ComputePrimVars(consVar, i, j, ibr); + } + + quokka::valarray Q_ip1{}; + quokka::valarray Q_ip2{}; + quokka::valarray Q_ip3{}; + quokka::valarray Q_ip4{}; + + // if gas is inflowing, change to reflecting B.C. + quokka::valarray Qr_im0 = detail::permute_vel(Q_i); + amrex::Real const v_normal = Qr_im0[1]; // normal velocity component + + if (((SIDE == BoundarySide::Lower) && (v_normal > 0.)) || ((SIDE == BoundarySide::Upper) && (v_normal < 0.))) { + // compute primitive vars + quokka::valarray Q_im1{}; + quokka::valarray Q_im2{}; + quokka::valarray Q_im3{}; + + if constexpr (DIR == FluxDir::X1) { + Q_im1 = HydroSystem::ComputePrimVars(consVar, im1, j, k); + Q_im2 = HydroSystem::ComputePrimVars(consVar, im2, j, k); + Q_im3 = HydroSystem::ComputePrimVars(consVar, im3, j, k); + } else if constexpr (DIR == FluxDir::X2) { + Q_im1 = HydroSystem::ComputePrimVars(consVar, i, im1, k); + Q_im2 = HydroSystem::ComputePrimVars(consVar, i, im2, k); + Q_im3 = HydroSystem::ComputePrimVars(consVar, i, im3, k); + } else if constexpr (DIR == FluxDir::X3) { + Q_im1 = HydroSystem::ComputePrimVars(consVar, i, j, im1); + Q_im2 = HydroSystem::ComputePrimVars(consVar, i, j, im2); + Q_im3 = HydroSystem::ComputePrimVars(consVar, i, j, im3); + } + + // reflect velocities + quokka::valarray Qr_im1 = detail::permute_vel(Q_im1); + quokka::valarray Qr_im2 = detail::permute_vel(Q_im2); + quokka::valarray Qr_im3 = detail::permute_vel(Q_im3); + + Qr_im0[1] *= -1.0; + Qr_im1[1] *= -1.0; + Qr_im2[1] *= -1.0; + Qr_im3[1] *= -1.0; + + Q_ip1 = detail::unpermute_vel(Qr_im0); + Q_ip2 = detail::unpermute_vel(Qr_im1); + Q_ip3 = detail::unpermute_vel(Qr_im2); + Q_ip4 = detail::unpermute_vel(Qr_im3); + } else { + // specify pressure at boundary + Q_i[4] = P_outflow; Q_ip1 = Q_i; Q_ip2 = Q_i; Q_ip3 = Q_i; Q_ip4 = Q_i; } + // set cell values + const int ip1 = (SIDE == BoundarySide::Lower) ? ibr - 1 : ibr + 1; + const int ip2 = (SIDE == BoundarySide::Lower) ? ibr - 2 : ibr + 2; + const int ip3 = (SIDE == BoundarySide::Lower) ? ibr - 3 : ibr + 3; + const int ip4 = (SIDE == BoundarySide::Lower) ? ibr - 4 : ibr + 4; + quokka::valarray consCell{}; if (idx[static_cast(DIR)] == ip1) { consCell = HydroSystem::ComputeConsVars(Q_ip1);