Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Enforcing inflow-outflow boundary solvability #129

Merged
merged 23 commits into from
Jul 21, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
23 commits
Select commit Hold shift + click to select a range
e97b7fa
adding vscode config file to gitignore
mukul1992 Feb 27, 2024
a0a1b7a
initial commit for the outflow correction
mukul1992 Feb 29, 2024
5107115
function for setting the mask values for flux calculation at the boun…
mukul1992 Mar 1, 2024
4d5c939
working on the routine to set inflow and outflow masks
mukul1992 Apr 20, 2024
5e83b8d
added the in- and out-flux computation
mukul1992 Apr 23, 2024
0d64ad3
added the outflow velocity correction routine
mukul1992 Apr 24, 2024
c43f075
replaced user_1 with new direction_dependent BCType in amrex
mukul1992 Jun 6, 2024
a672107
put in MPI reduction and corrected a bug in set_masks
mukul1992 Jun 6, 2024
b8ea32b
clean up for draft PR
mukul1992 Jun 7, 2024
3115f03
Added ability to enforce inflow-outflow solvability for cell-centered…
mukul1992 Jun 23, 2024
9cb268a
cleanup and comments
mukul1992 Jun 23, 2024
741c157
changing variable names to make functions generic rather than being s…
mukul1992 Jul 1, 2024
1f2c5f5
added the loop over levels
mukul1992 Jul 1, 2024
e2fec50
cleaned up mac-specific variable names, comments, and debug printouts
mukul1992 Jul 1, 2024
c4bd659
fixed unnecessary include
mukul1992 Jul 2, 2024
9ce8942
added new file to make.package for gnu make users
mukul1992 Jul 2, 2024
e057065
added docs for in-out solvability
mukul1992 Jul 2, 2024
6d567d5
replacing separate inflow and outflow masks with a single mask
mukul1992 Jul 3, 2024
1a529b8
adds clarification to the docs regarding the function not accounting …
mukul1992 Jul 3, 2024
54ee91e
Merge branch 'development' into correct_vels
asalmgren Jul 21, 2024
3bb2f8d
Update hydro_utils.H
asalmgren Jul 21, 2024
71ba07b
Update hydro_enforce_inout_solvability.cpp
asalmgren Jul 21, 2024
b5c73ad
Update hydro_enforce_inout_solvability.cpp
asalmgren Jul 21, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -54,3 +54,5 @@ Tutorials/EB/Donut/Exec/ls_plt*

# local config
Tools/GNUMake/Make.local

.vscode*
28 changes: 28 additions & 0 deletions Docs/source/InOutSolvability.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@
.. include:: CustomCommands.rst

Enforcing inflow-outflow solvability
------------------------------------

This routine enforces solvability for inflow-outflow boundaries,
which have both inflow and outflow cells.
A flux correction factor, :math:`\alpha_\text{fcf}` is introduced
which scales the outflow velocities to match the inflow:

.. math::

\sum_{\partial\Omega_\text{in}} {\bf u} \cdot {\bf \area} + \alpha_\text{fcf} \sum_{\partial\Omega_\text{out}} {\bf u} \cdot {\bf \area} = 0


The new flux-conserving velocities to be used for the MAC/nodal projections,
:math:`{\bf u}_\text{fc}`, are calculated from the correction factor as:

.. math::
\alpha_\text{fcf} = \frac{-\sum_{\partial\Omega_\text{in}} {\bf u} \cdot {\bf \area}}{\sum_{\partial\Omega_\text{out}} {\bf u} \cdot {\bf \area}},

.. math::
{\bf u}_\text{fc} = \alpha_\text{fcf} \cdot {\bf u}, \ \forall {\bf x} \in \partial\Omega_\text{out}.

It must be noted that this routine currently only accounts for boundaries
with the math BC ``BCType::direction_dependent``, which is to be used for
an inflow-outflow boundary. It does not compute or correct the boundary velocities
over other math BC types such as those representing pure inflow or pure outflow.
1 change: 1 addition & 0 deletions Docs/source/Utilities.rst
Original file line number Diff line number Diff line change
Expand Up @@ -38,3 +38,4 @@ Note that EB is only an option for Cartesian geometries as of this writing.

.. include:: AdvectiveTerm.rst

.. include:: InOutSolvability.rst
1 change: 1 addition & 0 deletions Utils/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -14,4 +14,5 @@ target_sources(
hydro_utils.cpp
hydro_constants.H
hydro_bcs_K.H
hydro_enforce_inout_solvability.cpp
cgilet marked this conversation as resolved.
Show resolved Hide resolved
)
1 change: 1 addition & 0 deletions Utils/Make.package
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
CEXE_sources += hydro_utils.cpp
CEXE_sources += hydro_compute_edgestate_and_flux.cpp
CEXE_sources += hydro_extrap_vel_to_faces.cpp
CEXE_sources += hydro_enforce_inout_solvability.cpp
CEXE_headers += hydro_bcs_K.H
CEXE_headers += hydro_utils.H

Expand Down
300 changes: 300 additions & 0 deletions Utils/hydro_enforce_inout_solvability.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,300 @@
/** \addtogroup Utilities
* @{
*/

#include <hydro_utils.H>

using namespace amrex;

namespace HydroUtils {

namespace {

void set_inout_masks(
const int lev,
const Vector<Array<MultiFab*, AMREX_SPACEDIM>>& vels_vec,
Array<iMultiFab, AMREX_SPACEDIM>& inout_masks,
const BCRec* bc_type,
const Box& domain,
const bool corners)
{
for (OrientationIter oit; oit != nullptr; ++oit) {
const auto ori = oit();
const auto side = ori.faceDir();
const int dir = ori.coordDir();
const auto islow = ori.isLow();
const auto ishigh = ori.isHigh();

// Multifab for normal velocity
const auto& vel_mf = vels_vec[lev][dir];

// mask iMF for the respective velocity direction
auto& inout_mask = inout_masks[dir];

IndexType::CellIndex dir_index_type = (vel_mf->ixType()).ixType(dir);
// domain extent indices for the velocities
int dlo;
if (dir_index_type == IndexType::CellIndex::CELL) {
// lower boundary is at -1 for cell-centered velocity
dlo = domain.smallEnd(dir) - 1;
} else {
// lower boundary is at 0 for face-centered velocity
dlo = domain.smallEnd(dir);
}
int dhi = domain.bigEnd(dir) + 1;

// get BCs for the normal velocity and set the boundary index
// based on low or high side
const BCRec ibcrec = bc_type[dir];
int bc, bndry;
if (side == Orientation::low) {
bc = ibcrec.lo(dir);
bndry = dlo;
} else {
bc = ibcrec.hi(dir);
bndry = dhi;
}

// limit influx/outflux calculations to the in-out boundaries only
if (bc == BCType::direction_dependent) {
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When I first saw "direction_dependent", it wasn't clear to me what this meant. The rest of the BCTypes are mathematical BCs. Why is "direction_dependent" the mathematical BC here instead of "mixed" or "Robin"?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We could remove it from amrex::BCType and use user_*.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

mixed might be okay too. But What mixed means is likely user dependent, whereas other types (e.g., reflect_odd, foextrap, etc.) are clear on what they mean and AMReX can fill them for the user.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We added this math BC type in AMReX (AMReX-Codes/amrex#3965) where the cells are to be treated differently based on in/out flow. Perhaps @asalmgren could comment on the naming.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To amrex, that's a user type that it does not know how to handle.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

On the user side, the user could define auto my_preferred_name = amrex::BCType::user_1;.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the reason for adding a special math type rather than just using user_1 was so that it can be used consistently across dependent codes, say, amr-wind and amrex-hydro, and so it does not conflict with other custom types that may be added in the future.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

amr-wind uses amrex-hydro. So it could use the types defined in amrex-hydro. (If it does not use amrex-hydro, then it probably does not matter.)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If this is meant to be used with both the Nodal and MAC projections, then since the BC types for those linops differ and I don't think that can be used for this purpose.

If this could be used more generally, in additional application codes, then maybe instead of having the BCRec it could be something more like a flag that indicates whether to include the face in the sum or not?

for (MFIter mfi(*vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) {

Box box = mfi.validbox();

// include ghost cells along normal velocity for cell-centered
// not for face-centered as boundary lies in valid region
if (dir_index_type == IndexType::CellIndex::CELL) {
box.grow(dir, 1);
}

// include boundary corners if specified
if (corners) {
box.grow((dir+1)%AMREX_SPACEDIM, 1);
box.grow((dir+2)%AMREX_SPACEDIM, 1);
}

// Enter further only if the box bndry is at the domain bndry
if ((islow && (box.smallEnd(dir) == dlo))
|| (ishigh && (box.bigEnd(dir) == dhi))) {

// create a 2D box normal to dir at the low/high bndry
Box box2d(box); box2d.setRange(dir, bndry);

auto vel_arr = vel_mf->array(mfi);
auto inout_mask_arr = inout_mask.array(mfi);

// tag cells as inflow or outflow by checking vel direction
ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
if ((side == Orientation::low && vel_arr(i,j,k) >= 0)
|| (side == Orientation::high && vel_arr(i,j,k) <= 0)) {
inout_mask_arr(i,j,k) = -1;
} else {
inout_mask_arr(i,j,k) = +1;
}
});
}
}
}
}
}

void compute_influx_outflux(
const int lev,
const Vector<Array<MultiFab*, AMREX_SPACEDIM>>& vels_vec,
const Array<iMultiFab, AMREX_SPACEDIM>& inout_masks,
const Real* a_dx,
Real& influx,
Real& outflux,
const bool corners)
{
influx = 0.0, outflux = 0.0;

for (int idim = 0; idim < AMREX_SPACEDIM; idim++) {

// normal face area
const Real ds =
a_dx[(idim+1) % AMREX_SPACEDIM] * a_dx[(idim+2) % AMREX_SPACEDIM];

// Multifab for normal velocity
const auto& vel_mf = vels_vec[lev][idim];

// grow in the respective direction if vel is cell-centered
IndexType index_type = vel_mf->ixType();
index_type.flip(idim); IntVect ngrow = index_type.ixType();

// grow in the transverse direction to include boundary corners
if (corners) {
ngrow[(idim+1)%AMREX_SPACEDIM] = 1;
ngrow[(idim+2)%AMREX_SPACEDIM] = 1;
}

// mask iMF for the respective velocity direction
const auto& inout_mask = inout_masks[idim];

// define "multi-arrays" and perform reduction using the mask
auto const& vel_ma = vel_mf->const_arrays();
auto const& inout_mask_ma = inout_mask.const_arrays();

influx += ds *
ParReduce(TypeList<ReduceOpSum>{},
TypeList<Real>{},
*vel_mf, ngrow,
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
noexcept -> GpuTuple<Real>
{
if (inout_mask_ma[box_no](i,j,k) == -1) {
return { std::abs(vel_ma[box_no](i,j,k)) };
} else {
return { 0. };
}
});

outflux += ds *
ParReduce(TypeList<ReduceOpSum>{},
TypeList<Real>{},
*vel_mf, ngrow,
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
noexcept -> GpuTuple<Real>
{
if (inout_mask_ma[box_no](i,j,k) == 1) {
return { std::abs(vel_ma[box_no](i,j,k)) };
} else {
return { 0. };
}
});
}
ParallelDescriptor::ReduceRealSum(influx);
ParallelDescriptor::ReduceRealSum(outflux);
}

void correct_outflow(
const int lev,
const Vector<Array<MultiFab*, AMREX_SPACEDIM>>& vels_vec,
const Array<iMultiFab, AMREX_SPACEDIM>& inout_masks,
const BCRec* bc_type,
const Box& domain,
const Real alpha_fcf,
const bool corners)
{
for (OrientationIter oit; oit != nullptr; ++oit) {
const auto ori = oit();
const auto side = ori.faceDir();
const int dir = ori.coordDir();
const auto islow = ori.isLow();
const auto ishigh = ori.isHigh();

// Multifab for normal velocity
const auto& vel_mf = vels_vec[lev][dir];

// mask iMF for the respective velocity direction
const auto& inout_mask = inout_masks[dir];

IndexType::CellIndex dir_index_type = (vel_mf->ixType()).ixType(dir);
// domain extent indices for the velocities
int dlo;
if (dir_index_type == IndexType::CellIndex::CELL) {
dlo = domain.smallEnd(dir) - 1; // cell-centered boundary
} else {
dlo = domain.smallEnd(dir); // face-centered boundary
}
int dhi = domain.bigEnd(dir) + 1;

// get BCs for the normal velocity and set the boundary index
const BCRec ibcrec = bc_type[dir];
int bc, bndry;
if (side == Orientation::low) {
bc = ibcrec.lo(dir);
bndry = dlo;
} else {
bc = ibcrec.hi(dir);
bndry = dhi;
}

if (bc == BCType::direction_dependent) {
for (MFIter mfi(*vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) {

Box box = mfi.validbox();
if (dir_index_type == IndexType::CellIndex::CELL) {
box.grow(dir, 1);
}
if (corners) {
box.grow((dir+1)%AMREX_SPACEDIM, 1);
box.grow((dir+2)%AMREX_SPACEDIM, 1);
}

// Enter further only if the box boundary is at the domain boundary
if ((islow && (box.smallEnd(dir) == dlo))
|| (ishigh && (box.bigEnd(dir) == dhi))) {

// create a 2D box normal to dir at the low/high boundary
Box box2d(box); box2d.setRange(dir, bndry);

auto vel_arr = vel_mf->array(mfi);
auto inout_mask_arr = inout_mask.array(mfi);

ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
if (inout_mask_arr(i,j,k) == 1) {
vel_arr(i,j,k) *= alpha_fcf;
}
});
}
}
}
}
}

} // file-local namespace

void enforceInOutSolvability (
const Vector<Array<MultiFab*, AMREX_SPACEDIM>>& vels_vec,
const BCRec* bc_type,
const Vector<Geometry>& geom,
bool include_bndry_corners
)
{
const Box domain = geom[0].Domain();

const auto nlevs = int(vels_vec.size());
for (int lev = 0; lev < nlevs; ++lev) {

// masks to tag inflow/outflow at the boundaries
// separate iMultifab for each velocity direction
// 0 for interior cells,
// -1 for inflow bndry cells, +1 for outflow bndry cells
Array<iMultiFab, AMREX_SPACEDIM> inout_masks;

// defining the mask iMultifabs in each direction
for (int idim = 0; idim < AMREX_SPACEDIM; idim++)
{
const auto& vel_mf = vels_vec[lev][idim]; // normal velocity multifab

// grow in the respective direction if vel is cell-centered
// to include the boundary cells
IndexType index_type = vel_mf->ixType();
index_type.flip(idim); IntVect ngrow = index_type.ixType();

// grow in the transverse direction to include boundary corners
if (include_bndry_corners) {
ngrow[(idim+1)%AMREX_SPACEDIM] = 1;
ngrow[(idim+2)%AMREX_SPACEDIM] = 1;
}

inout_masks[idim].define(vel_mf->boxArray(), vel_mf->DistributionMap(), 1, ngrow);
inout_masks[idim].setVal(0);
}

set_inout_masks(lev, vels_vec, inout_masks, bc_type, domain, include_bndry_corners);

const Real* a_dx = geom[lev].CellSize();
Real influx = 0.0, outflux = 0.0;
compute_influx_outflux(lev, vels_vec, inout_masks, a_dx, influx, outflux, include_bndry_corners);

const Real alpha_fcf = influx/outflux; // flux correction factor
correct_outflow(lev, vels_vec, inout_masks, bc_type, domain, alpha_fcf, include_bndry_corners);

} // levels loop
}

}
11 changes: 11 additions & 0 deletions Utils/hydro_utils.H
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,17 @@
*/

namespace HydroUtils {

/**
* \brief Enforces solvablity by scaling outflow to match with inflow.
*
*/
void enforceInOutSolvability (
const amrex::Vector<amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>>& vels_vec,
amrex::BCRec const* bc_type,
const amrex::Vector<amrex::Geometry>& geom,
bool include_bndry_corners = false);

/**
* \brief Compute edge state and flux. Most general version for use with multilevel synchronization.
* All other versions ultimately call this one.
Expand Down