Skip to content

Commit

Permalink
add simple subsonic inflow/outflow BCs (#598)
Browse files Browse the repository at this point in the history 
### 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 #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.
  • Loading branch information
@BenWibking
BenWibking authored Apr 8, 2024
1 parent d13e3bc commit dd5c2fd
Show file tree
Hide file tree
Showing 2 changed files with 154 additions and 13 deletions.
53 changes: 52 additions & 1 deletion src/NSCBC_inflow.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Expand All @@ -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<problem_t>::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.;
Expand Down Expand Up @@ -145,6 +147,55 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE void setInflowX1Lower(const amrex::IntVect &
consVar(i, j, k, HydroSystem<problem_t>::scalar0_index + n) = consCell[6 + n];
}
}

template <typename problem_t>
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void setInflowX1LowerLowOrder(const amrex::IntVect &iv, amrex::Array4<amrex::Real> 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<Real, HydroSystem<problem_t>::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<problem_t>::nvar_;

/// x1 lower boundary -- subsonic inflow

// compute primitive vars from valid region
quokka::valarray<amrex::Real, N> const Q_i = HydroSystem<problem_t>::ComputePrimVars(consVar, ilo, j, k);

// compute centered ghost values
const Real rho = Q_i[0];
const Real Eint = quokka::EOS<problem_t>::ComputeEintFromTgas(rho, T_t);
quokka::valarray<amrex::Real, N> 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<problem_t>::ComputePressure(rho, Eint); // prescribe temperature
Q_im1[5] = Eint; // prescribe temperature
for (int i = 0; i < HydroSystem<problem_t>::nscalars_; ++i) {
Q_im1[6 + i] = s_t[i]; // prescribe passive scalars
}

// set cell values
quokka::valarray<amrex::Real, N> consCell{};
consCell = HydroSystem<problem_t>::ComputeConsVars(Q_im1);

consVar(i, j, k, HydroSystem<problem_t>::density_index) = consCell[0];
consVar(i, j, k, HydroSystem<problem_t>::x1Momentum_index) = consCell[1];
consVar(i, j, k, HydroSystem<problem_t>::x2Momentum_index) = consCell[2];
consVar(i, j, k, HydroSystem<problem_t>::x3Momentum_index) = consCell[3];
consVar(i, j, k, HydroSystem<problem_t>::energy_index) = consCell[4];
consVar(i, j, k, HydroSystem<problem_t>::internalEnergy_index) = consCell[5];
for (int n = 0; n < HydroSystem<problem_t>::nscalars_; ++n) {
consVar(i, j, k, HydroSystem<problem_t>::scalar0_index + n) = consCell[6 + n];
}
}
} // namespace NSCBC

#endif // NSCBC_INFLOW_HPP_
114 changes: 102 additions & 12 deletions src/NSCBC_outflow.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -249,15 +249,6 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto unpermute_vel(quokka::valarray<amrex::R
}
return newPrim;
}

template <typename problem_t>
AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto isStateValid(quokka::valarray<amrex::Real, HydroSystem<problem_t>::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 <typename problem_t, FluxDir DIR, BoundarySide SIDE>
Expand Down Expand Up @@ -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<problem_t>(Q_ip1) && detail::isStateValid<problem_t>(Q_ip2) && detail::isStateValid<problem_t>(Q_ip3) &&
detail::isStateValid<problem_t>(Q_ip4))) {
quokka::valarray<amrex::Real, N> consCell{};
if (idx[static_cast<int>(DIR)] == ip1) {
consCell = HydroSystem<problem_t>::ComputeConsVars(Q_ip1);
} else if (idx[static_cast<int>(DIR)] == ip2) {
consCell = HydroSystem<problem_t>::ComputeConsVars(Q_ip2);
} else if (idx[static_cast<int>(DIR)] == ip3) {
consCell = HydroSystem<problem_t>::ComputeConsVars(Q_ip3);
} else if (idx[static_cast<int>(DIR)] == ip4) {
consCell = HydroSystem<problem_t>::ComputeConsVars(Q_ip4);
}

consVar(i, j, k, HydroSystem<problem_t>::density_index) = consCell[0];
consVar(i, j, k, HydroSystem<problem_t>::x1Momentum_index) = consCell[1];
consVar(i, j, k, HydroSystem<problem_t>::x2Momentum_index) = consCell[2];
consVar(i, j, k, HydroSystem<problem_t>::x3Momentum_index) = consCell[3];
consVar(i, j, k, HydroSystem<problem_t>::energy_index) = consCell[4];
consVar(i, j, k, HydroSystem<problem_t>::internalEnergy_index) = consCell[5];
for (int n = 0; n < HydroSystem<problem_t>::nscalars_; ++n) {
consVar(i, j, k, HydroSystem<problem_t>::scalar0_index + n) = consCell[6 + n];
}
}

template <typename problem_t, FluxDir DIR, BoundarySide SIDE>
AMREX_GPU_DEVICE AMREX_FORCE_INLINE void setOutflowBoundaryLowOrder(const amrex::IntVect &iv, amrex::Array4<amrex::Real> 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<int, 3> idx{i, j, k};
amrex::Box const &box = geom.Domain();
constexpr int N = HydroSystem<problem_t>::nvar_;

const auto &boundary_idx = (SIDE == BoundarySide::Lower) ? box.loVect3d() : box.hiVect3d();
const int ibr = boundary_idx[static_cast<int>(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<int>(DIR));

// compute primitive vars
quokka::valarray<amrex::Real, N> Q_i{};
if constexpr (DIR == FluxDir::X1) {
Q_i = HydroSystem<problem_t>::ComputePrimVars(consVar, ibr, j, k);
} else if constexpr (DIR == FluxDir::X2) {
Q_i = HydroSystem<problem_t>::ComputePrimVars(consVar, i, ibr, k);
} else if constexpr (DIR == FluxDir::X3) {
Q_i = HydroSystem<problem_t>::ComputePrimVars(consVar, i, j, ibr);
}

quokka::valarray<amrex::Real, N> Q_ip1{};
quokka::valarray<amrex::Real, N> Q_ip2{};
quokka::valarray<amrex::Real, N> Q_ip3{};
quokka::valarray<amrex::Real, N> Q_ip4{};

// if gas is inflowing, change to reflecting B.C.
quokka::valarray<amrex::Real, N> Qr_im0 = detail::permute_vel<problem_t, DIR>(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<amrex::Real, N> Q_im1{};
quokka::valarray<amrex::Real, N> Q_im2{};
quokka::valarray<amrex::Real, N> Q_im3{};

if constexpr (DIR == FluxDir::X1) {
Q_im1 = HydroSystem<problem_t>::ComputePrimVars(consVar, im1, j, k);
Q_im2 = HydroSystem<problem_t>::ComputePrimVars(consVar, im2, j, k);
Q_im3 = HydroSystem<problem_t>::ComputePrimVars(consVar, im3, j, k);
} else if constexpr (DIR == FluxDir::X2) {
Q_im1 = HydroSystem<problem_t>::ComputePrimVars(consVar, i, im1, k);
Q_im2 = HydroSystem<problem_t>::ComputePrimVars(consVar, i, im2, k);
Q_im3 = HydroSystem<problem_t>::ComputePrimVars(consVar, i, im3, k);
} else if constexpr (DIR == FluxDir::X3) {
Q_im1 = HydroSystem<problem_t>::ComputePrimVars(consVar, i, j, im1);
Q_im2 = HydroSystem<problem_t>::ComputePrimVars(consVar, i, j, im2);
Q_im3 = HydroSystem<problem_t>::ComputePrimVars(consVar, i, j, im3);
}

// reflect velocities
quokka::valarray<amrex::Real, N> Qr_im1 = detail::permute_vel<problem_t, DIR>(Q_im1);
quokka::valarray<amrex::Real, N> Qr_im2 = detail::permute_vel<problem_t, DIR>(Q_im2);
quokka::valarray<amrex::Real, N> Qr_im3 = detail::permute_vel<problem_t, DIR>(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<problem_t, DIR>(Qr_im0);
Q_ip2 = detail::unpermute_vel<problem_t, DIR>(Qr_im1);
Q_ip3 = detail::unpermute_vel<problem_t, DIR>(Qr_im2);
Q_ip4 = detail::unpermute_vel<problem_t, DIR>(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<amrex::Real, N> consCell{};
if (idx[static_cast<int>(DIR)] == ip1) {
consCell = HydroSystem<problem_t>::ComputeConsVars(Q_ip1);
Expand Down

0 comments on commit dd5c2fd

Please sign in to comment.