Skip to content

Commit

Permalink
Add mixed bc and update docs (#89)
Browse files Browse the repository at this point in the history 
This allows inflow and outflow on the same domain face, where the Dirichlet and Neumann regions are separated by an EB.
Depends on Hydro pr 126 and AMReX PR 3788
  • Loading branch information
@cgilet
cgilet authored Mar 15, 2024
1 parent ae8f233 commit e78cd29
Show file tree
Hide file tree
Showing 26 changed files with 1,109 additions and 118 deletions.
3 changes: 3 additions & 0 deletions Docs/sphinx_documentation/source/InputsAlgorithm.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,4 +26,7 @@ The following inputs must be preceded by "incflo."
+----------------------+-----------------------------------------------------------------------+-------------+--------------+
| mu_s | scalar diffusivity | Real(s) | 0.0 |
+----------------------+-----------------------------------------------------------------------+-------------+--------------+
| use_tensor_solve | In velocity solve, use multicomponent :math:`\nabla \cdot \tau` | bool | true |
| | otherwise use separate solves for each velocity component | | |
+----------------------+-----------------------------------------------------------------------+-------------+--------------+

10 changes: 10 additions & 0 deletions Docs/sphinx_documentation/source/InputsProblemDefinition.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,7 @@ Setting basic boundary conditions can be specified by inputs preceded by "xlo",
| | * 'po' or 'pressure_outflow' | | |
| | * 'mi' or 'mass_inflow' | | |
| | * 'nsw' or 'no_slip_wall' | | |
| | * 'mixed' | | |
+--------------------+---------------------------------------------------------------------------+-------------+-----------+
| pressure | Sets boundary pressure for pressure inflows, outflows and mass inflows | Real | None |
+--------------------+---------------------------------------------------------------------------+-------------+-----------+
Expand All @@ -60,5 +61,14 @@ Setting basic boundary conditions can be specified by inputs preceded by "xlo",
| tracer | Sets boundary tracer for mass inflows | Real | 0.0 |
+--------------------+---------------------------------------------------------------------------+-------------+-----------+

The 'mixed' boundary type allows for inflow and outflow on the same domain face.
The implementation requires that there is EB separating the inflow and outflow regions,
and only currently treats constant viscosity (i.e. requires the inputs file contains ``incflo.use_tensor_solve = false``).
It does allow for non-constant scalar diffusivity. To create a new problem setup with mixed BCs, one must additionally
specify the Dirchlet and Neumann areas in ``prob/prob_bc.cpp`` In this file there are functions to set the needed BC info for the diffusion solver, the MAC projection, and with the same function, the nodal projection and advection. Comments in the code provide a detailed explaination for each instance. For addditional details on mixed BCs, also see AMReX-Hydro's documentation (:ref:`bcs`).

..
To create a new setup:
initial conditions go in prob/prob_init_fluid.cpp
inflow boundary conditions are set in prob/prob_bc.H -- the BCType for mixed bc is set at foextrap here to allow the outflow to get appropriately filled, which would already have happened before reaching these functions.
43 changes: 40 additions & 3 deletions src/boundary_conditions/boundary_conditions.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,39 @@ void incflo::init_bcs ()
// tracer has Dirichlet bcs
pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac);
}
else if (bc_type == "mixed" )
{
amrex::Print() << bcid << " set to mixed inflow outflow.\n";
m_has_mixedBC = true;
#ifdef AMREX_USE_EB
// ReadParameters() already called
if (m_advection_type != "Godunov") { amrex::Abort("mixed BCs require Godunov"); }

ParmParse ipp("incflo");
std::string eb_geom = "null";
ipp.query("geometry", eb_geom);
eb_geom = amrex::toLower(eb_geom);
if (eb_geom == "null" || eb_geom == "all_regular")
#endif
{
Abort("For now, mixed BCs must be separated by an EB");
}
Warning("Using BC type mixed requires that the Dirichlet and Neumann regions are separated by EB.");

m_bc_type[ori] = BC::mixed;

pp.get("pressure", m_bc_pressure[ori]);

std::vector<Real> v;
if (pp.queryarr("velocity", v, 0, AMREX_SPACEDIM)) {
for (int i=0; i<AMREX_SPACEDIM; i++){
m_bc_velocity[ori][i] = v[i];
}
}

pp.query("density", m_bc_density[ori]);
pp.queryarr("tracer", m_bc_tracer[ori], 0, m_ntrac);
}
else
{
m_bc_type[ori] = BC::undefined;
Expand Down Expand Up @@ -146,7 +179,9 @@ void incflo::init_bcs ()
Orientation::Side side = ori.faceDir();
auto const bct = m_bc_type[ori];
if (bct == BC::pressure_inflow ||
bct == BC::pressure_outflow)
bct == BC::pressure_outflow ||
bct == BC::mixed ) // BC_mask will handle Dirichlet part.

{
if (side == Orientation::low) {
AMREX_D_TERM(m_bcrec_velocity[0].setLo(dir, BCType::foextrap);,
Expand Down Expand Up @@ -220,7 +255,8 @@ void incflo::init_bcs ()
auto const bct = m_bc_type[ori];
if (bct == BC::pressure_inflow ||
bct == BC::pressure_outflow ||
bct == BC::no_slip_wall)
bct == BC::no_slip_wall ||
bct == BC::mixed)
{
if (side == Orientation::low) {
m_bcrec_density[0].setLo(dir, BCType::foextrap);
Expand Down Expand Up @@ -271,7 +307,8 @@ void incflo::init_bcs ()
Orientation::Side side = ori.faceDir();
auto const bct = m_bc_type[ori];
if (bct == BC::pressure_inflow ||
bct == BC::pressure_outflow)
bct == BC::pressure_outflow ||
bct == BC::mixed)
{
if (side == Orientation::low) {
for (auto& b : m_bcrec_tracer) b.setLo(dir, BCType::foextrap);
Expand Down
171 changes: 171 additions & 0 deletions src/boundary_conditions/incflo_set_bcs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,177 @@

using namespace amrex;

// Make the iMultiFab needed by the Nodal Projection to support mixed BCs
iMultiFab
incflo::make_nodalBC_mask(int lev)
{
// We use the solver's overset mask to implement mixed BCs. The mask is nodal, and
// so exists on the boundary face, where 0 indicates a Dirichlet boundary condition
// (i,e, no solve needed for that point because we know the solution there already).
// Because it's an overset mask and not a BC, we must fill the entire MF.
iMultiFab new_mask(amrex::convert(grids[lev], IntVect::TheNodeVector()), dmap[lev], 1, 0);
new_mask = 1;
int inflow = 1;
int outflow = 0;

Geometry const& gm = Geom(lev);
Box const& domain = gm.Domain();
for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
Orientation olo(dir,Orientation::low);
Orientation ohi(dir,Orientation::high);
if (m_bc_type[olo] == BC::mixed || m_bc_type[ohi] == BC::mixed) {
Box dlo = (m_bc_type[olo] == BC::mixed) ? surroundingNodes(bdryLo(domain,dir)) : Box().setType(IndexType::TheNodeType());
Box dhi = (m_bc_type[ohi] == BC::mixed) ? surroundingNodes(bdryHi(domain,dir)) : Box().setType(IndexType::TheNodeType());
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
for (MFIter mfi(new_mask); mfi.isValid(); ++mfi) {
Box blo = mfi.validbox() & dlo;
Box bhi = mfi.validbox() & dhi;
Array4<int> const& mask_arr = new_mask.array(mfi);
if (blo.ok()) {
prob_set_BC_MF(olo, blo, mask_arr, lev, inflow, outflow, "projection");
}
if (bhi.ok()) {
prob_set_BC_MF(ohi, bhi, mask_arr, lev, inflow, outflow, "projection");
}
}
}
}

return new_mask;
}

// Make a position dependent BC MultiFab for the advection routines to support
// mixed BCs.
std::unique_ptr<iMultiFab>
incflo::make_BC_MF(int lev, amrex::Gpu::DeviceVector<amrex::BCRec> const& bcs,
std::string const& field)
{
auto ncomp = static_cast<int>(bcs.size());
// The advection routines expect that the BC type is stored in the first ghost
// cell of a cell-centered MF.
std::unique_ptr<iMultiFab> BC_MF(new iMultiFab(grids[lev], dmap[lev], ncomp, 1));
// initialize to 0 == BCType::int_dir, not sure this is needed really
*BC_MF = 0;
// For advection, we use FOEXTRAP for outflow regions per incflo::init_bcs()
int inflow = BCType::ext_dir;
int outflow = BCType::foextrap;

Box const& domain = geom[lev].Domain();
for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
Orientation olo(dir,Orientation::low);
Orientation ohi(dir,Orientation::high);

Box dlo = adjCellLo(domain,dir);
Box dhi = adjCellHi(domain,dir);
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
for (MFIter mfi(*BC_MF); mfi.isValid(); ++mfi) {
Box blo = mfi.growntilebox() & dlo;
Box bhi = mfi.growntilebox() & dhi;
Array4<int> const& bc_arr = BC_MF->array(mfi);
BCRec const* const bc_ptr = bcs.data();
if (m_bc_type[olo] == BC::mixed) {
prob_set_BC_MF(olo, blo, bc_arr, lev, inflow, outflow, field);
} else {
amrex::ParallelFor(blo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
for (int n = 0; n < ncomp; n++) {
bc_arr(i,j,k,n) = bc_ptr[n].lo(dir);
}
});
}
if (m_bc_type[ohi] == BC::mixed) {
prob_set_BC_MF(ohi, bhi, bc_arr, lev, inflow, outflow, field);
} else {
amrex::ParallelFor(bhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
for (int n = 0; n < ncomp; n++) {
bc_arr(i,j,k,n) = bc_ptr[n].hi(dir);
}
});
}
}
}

return BC_MF;
}

// Create the MutliFabs needed by the linear solver to define Robin BCs, which incflo
// uses to support mixed BCs in the MAC Projection and diffusion.
Vector<MultiFab>
incflo::make_robinBC_MFs(int lev, MultiFab* state)
{
// MLMG defines the Robin BC with
// a u + b du/dn = f
// The values are stored in the first ghost cell, although the BC is considered on
// the face. Only the ghost cells of robin BC arrays are used in MLMG.
int nghost = 1;

Vector<MultiFab> robin(3);
for (int n = 0; n < 3; n++) {
robin[n].define(grids[lev], dmap[lev], 1, nghost);
}

MultiFab& robin_a = robin[0];
MultiFab& robin_b = robin[1];
MultiFab& robin_f = robin[2];

Box const& domain = Geom(lev).Domain();
for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
Orientation olo(dir,Orientation::low);
Orientation ohi(dir,Orientation::high);

// Only need to fill Robin BC sides, MLMG will check for Robin BC first
if (m_bc_type[olo] == BC::mixed || m_bc_type[ohi] == BC::mixed) {
Box dlo = (m_bc_type[olo] == BC::mixed) ? adjCellLo(domain,dir) : Box();
Box dhi = (m_bc_type[ohi] == BC::mixed) ? adjCellHi(domain,dir) : Box();
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
// FIXME - do we want tiling here...
for (MFIter mfi(robin_a,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
Box const& gbx = amrex::grow(mfi.validbox(),nghost);
Box blo = gbx & dlo;
Box bhi = gbx & dhi;
Array4<Real> const& a_arr = robin_a.array(mfi);
Array4<Real> const& b_arr = robin_b.array(mfi);
Array4<Real> const& f_arr = robin_f.array(mfi);

if (state) {
// state has filled boundry, so can use its ghost cells for the BC values
Array4<Real const> const& bcv = state->const_array(mfi);
// Diffusion, solving for velocity, density, tracer
//
// Robin BC: a u + b du/dn = f -- inflow, Dirichlet a=1, b=0, f=bcv
// -- outflow, Neumann a=0, b=1, f=0
if (blo.ok()) {
prob_set_diffusion_robinBCs(olo, blo, a_arr, b_arr, f_arr, bcv, lev);
}
if (bhi.ok()) {
prob_set_diffusion_robinBCs(ohi, bhi, a_arr, b_arr, f_arr, bcv, lev);
}
} else {
// MAC, solving for phi (~pressure)
//
// Robin BC: a u + b du/dn = f -- inflow, Neumann a=0, b=1, f=0
// -- outflow, Dirichlet a=1, b=0, f=0
if (blo.ok()) {
prob_set_MAC_robinBCs(olo, blo, a_arr, b_arr, f_arr, lev);
}
if (bhi.ok()) {
prob_set_MAC_robinBCs(ohi, bhi, a_arr, b_arr, f_arr, lev);
}
}
}
}
}

return robin;
}

#ifdef AMREX_USE_EB
void
incflo::set_eb_velocity (int lev, amrex::Real /*time*/, MultiFab& eb_vel, int nghost)
Expand Down
27 changes: 24 additions & 3 deletions src/convection/incflo_compute_MAC_projected_velocities.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,14 +110,23 @@ incflo::compute_MAC_projected_velocities (
macproj->initProjector(lp_info, GetVecOfArrOfConstPtrs(inv_rho));
}
macproj->setDomainBC(get_mac_projection_bc(Orientation::low), get_mac_projection_bc(Orientation::high));

if ( m_has_mixedBC ) {
for (int lev = 0; lev <= finest_level; ++lev)
{
auto const robin = make_robinBC_MFs(lev);
macproj->setLevelBC(lev, nullptr,
&robin[0], &robin[1], &robin[2]);
}
}
} else {
#ifndef AMREX_USE_EB
if (m_constant_density) {
macproj->updateBeta(l_dt/m_ro_0); // unnecessary unless m_ro_0 changes.
} else
#endif
{
macproj->updateBeta(GetVecOfArrOfConstPtrs(inv_rho));
macproj->updateCoeffs(GetVecOfArrOfConstPtrs(inv_rho));
}
}

Expand All @@ -138,6 +147,19 @@ incflo::compute_MAC_projected_velocities (
l_advection_type = "Godunov";
}

std::unique_ptr<iMultiFab> BC_MF;
if (m_has_mixedBC) {
// Create a MF to hold the BCType info. Note that this is different than the
// bcs for the MAC projection because the MAC operates on phi, this is velocity.
//
// TODO? Could consider creating an incflo member variable to save the BC_MF
// amrex::Vector<std::unique_ptr<amrex::iMultiFab> > m_BC_MF;
// Could stack components as vel, density, tracer, and then could use for
// scalars' advective step as well. But not sure it really matters one way or
// the other.
BC_MF = make_BC_MF(lev, m_bcrec_velocity_d, "velocity");
}

// Predict normal velocity to faces -- note that the {u_mac, v_mac, w_mac}
// returned from this call are on face CENTROIDS
HydroUtils::ExtrapVelToFaces(*vel[lev], *vel_forces[lev],
Expand All @@ -149,7 +171,7 @@ incflo::compute_MAC_projected_velocities (
m_eb_flow.enabled ? get_velocity_eb()[lev] : nullptr,
#endif
m_godunov_ppm, m_godunov_use_forces_in_trans,
l_advection_type);
l_advection_type, PPM::default_limiter, BC_MF.get());
}

Vector<Array<MultiFab*,AMREX_SPACEDIM> > mac_vec(finest_level+1);
Expand Down Expand Up @@ -190,7 +212,6 @@ incflo::compute_MAC_projected_velocities (
} else {
macproj->project(m_mac_mg_rtol,m_mac_mg_atol);
}

// Note that the macproj->project call above ensures that the MAC velocities are averaged down --
// we don't need to do that again here
}
Loading

0 comments on commit e78cd29

Please sign in to comment.