Skip to content

Commit

Permalink
More clean up. Improve a few comments.
Browse files Browse the repository at this point in the history
  • Loading branch information
@cgilet
cgilet committed Feb 22, 2024
1 parent 55f38c9 commit 9e82cf4
Show file tree
Hide file tree
Showing 7 changed files with 59 additions and 226 deletions.
149 changes: 37 additions & 112 deletions src/boundary_conditions/incflo_set_bcs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,42 +6,6 @@

using namespace amrex;

// iMultiFab
// incflo::make_ccBC_mask(int lev)
// {
// // BC info stored in ghost cell.
// iMultiFab new_mask(grids[lev], dmap[lev], 1, 1);
// // No expectation in MLMG for there to be valid data where there is no Robin BC
// //new_mask = 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) ? 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
// for (MFIter mfi(new_mask); mfi.isValid(); ++mfi) {
// Box blo = mfi.growntilebox() & dlo;
// Box bhi = mfi.growntilebox() & dhi;
// Array4<int> const& mask_arr = new_mask.array(mfi);
// if (blo.ok()) {
// prob_set_BC_MF(olo, blo, mask_arr, lev);
// }
// if (bhi.ok()) {
// prob_set_BC_MF(ohi, bhi, mask_arr, lev);
// }
// }
// }
// }

// return new_mask;
// }

iMultiFab
incflo::make_nodalBC_mask(int lev)
{
Expand Down Expand Up @@ -81,6 +45,7 @@ incflo::make_nodalBC_mask(int lev)
return new_mask;
}

// This function makes a position dependent BC MF for the advection routines
std::unique_ptr<iMultiFab>
incflo::make_BC_MF(int lev, amrex::Gpu::DeviceVector<amrex::BCRec>& bcs,
std::string field)
Expand Down Expand Up @@ -133,57 +98,68 @@ incflo::make_BC_MF(int lev, amrex::Gpu::DeviceVector<amrex::BCRec>& bcs,
return BC_MF;
}

Vector<std::unique_ptr<MultiFab>>
incflo::make_diffusion_robinBC_MFs(int lev, MultiFab& phi)
// TODO There's probably a way to combine this with above function rather than repeating
// much of the same code...
Vector<MultiFab>
incflo::make_robinBC_MFs(int lev, MultiFab* phi)
{
int nghost = 1;

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

MultiFab& robin_a = *robin[0];
MultiFab& robin_b = *robin[1];
MultiFab& robin_f = *robin[2];
// fixme - for vis, init robin MFs, although i don't know that this is generally needed
robin_a = 0;
robin_b = 0;
robin_f = 0;
MultiFab& robin_a = robin[0];
MultiFab& robin_b = robin[1];
MultiFab& robin_f = robin[2];

// only ghost cells of robin BC arrays are used in MLMG
// bc in ghost cells that are outside the domain.
// define and fill Robin BC a, b, and f
// BC values are stored in the ghost cells, although the BC is considered on face.
// Only the ghost cells of robin BC arrays are used in MLMG.
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 - don't think we want tiling here...
// 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);
Array4<Real const> const& bcv = phi.const_array(mfi);

// Robin BC: a u + b du/dn = f -- inflow, Dirichlet a=1, b=0, f=bcv
// -- outflow, Neumann a=0, b=1, f=0
// Only need to fill Robin BC sides, MLMG will check for Robin BC first
if (blo.ok()) {
// fixme -- here i think we need to go to physbcfunct here
// phi has already had boundry filled, so can just use that here
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);
if (phi) {
// phi has filled boundry, so can use its ghost cells for the BC values
Array4<Real const> const& bcv = phi->const_array(mfi);
// Diffusion
// 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
// 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);
}
}
}
}
Expand All @@ -192,57 +168,6 @@ incflo::make_diffusion_robinBC_MFs(int lev, MultiFab& phi)
return robin;
}

// void
// incflo::make_nodalBC_mask(int lev,
// const BoxArray& ba, // do we really need to pass these, use member vars? - i think we need to pass to be safe when calling from RemakeLevel...
// const DistributionMapping& dm)
// {
// // FIXME do we really want this check here?
// bool has_mixedBC = false;
// for (OrientationIter oit; oit; ++oit) {
// if (m_bc_type[oit()] == BC::mixed )
// {
// has_mixedBC = true;
// break;
// }
// }

// if (!has_mixedBC) { return; }


// // MLNodeLap does not require any ghost cells...
// std::unique_ptr<iMultiFab> new_mask(new iMultiFab(amrex::convert(ba,IntVect::TheNodeVector()),
// dm, 1, 0));
// *new_mask = 1;

// 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();
// Box dhi = (m_bc_type[ohi] == BC::mixed) ? surroundingNodes(bdryHi(domain,dir)) : Box();
// #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);
// }
// if (bhi.ok()) {
// prob_set_BC_MF(ohi, bhi, mask_arr, lev);
// }
// }
// }
// }

// m_BC_MF[lev] = std::move(new_mask);
// }

#ifdef AMREX_USE_EB
void
incflo::set_eb_velocity (int lev, amrex::Real /*time*/, MultiFab& eb_vel, int nghost)
Expand Down
68 changes: 11 additions & 57 deletions src/convection/incflo_compute_MAC_projected_velocities.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,6 @@ incflo::compute_MAC_projected_velocities (
//
// Initialize (or redefine the beta in) the MacProjector
//
Print()<<"Mac proj needs init? "<<macproj->needInitialization()<<std::endl;
if (macproj->needInitialization())
{
LPInfo lp_info;
Expand All @@ -112,61 +111,12 @@ incflo::compute_MAC_projected_velocities (
}
macproj->setDomainBC(get_mac_projection_bc(Orientation::low), get_mac_projection_bc(Orientation::high));

// Not sure exactly where this goes. Does it only need initialization? set every time???
// initProj only take as many levels as 1/rho, so if a level is added, would have to
// re-init proj anyway.
// robin_a, etc needs to exist to finest_level
// Values need to be in the ghost cells, although the BC is considered on face.
if ( m_has_mixedBC ) {
int nghost = 1;

for (int lev = 0; lev <= finest_level; ++lev)
{
// define and fill Robin BC a, b, and f
MultiFab robin_a (grids[lev],dmap[lev],1,nghost,MFInfo(),Factory(lev));
MultiFab robin_b (grids[lev],dmap[lev],1,nghost,MFInfo(),Factory(lev));
MultiFab robin_f (grids[lev],dmap[lev],1,nghost,MFInfo(),Factory(lev));

// fixme? - for vis, init robin MFs, although i don't know that this is generally needed
robin_a = 0;
robin_b = 0;
robin_f = 0;

// only ghost cells of robin BC arrays are used in MLMG
// bc in ghost cells that are outside the domain.
Box const& domain = Geom(lev).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) ? 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);

// Robin BC: a u + b du/dn = f -- inflow, Neumann a=0, b=1, f=0
// -- outflow, Dirichlet a=1, b=0, f=0
// Only need to fill Robin BC sides, MLMG will check for Robin BC first
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);
}
}
}
}
macproj->setLevelBC(lev, nullptr, &robin_a, &robin_b, &robin_f);
auto const robin = make_robinBC_MFs(lev);
macproj->setLevelBC(lev, nullptr,
&robin[0], &robin[1], &robin[2]);
}
}
} else {
Expand Down Expand Up @@ -197,12 +147,16 @@ incflo::compute_MAC_projected_velocities (
l_advection_type = "Godunov";
}

// FIXME? - need to create the BC MF here. This is a iMF vs a Real MF, so can't use
// the robin BC info from MAC. This is because we need this to hold the BCType,
// which are enums, so ints
//
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");
}

Expand Down
1 change: 0 additions & 1 deletion src/convection/incflo_compute_advection_term.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -372,7 +372,6 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
is_velocity, fluxes_are_area_weighted,
m_advection_type, PPM::default_limiter, velbc_arr);


// ************************************************************************
// Density
// ************************************************************************
Expand Down
16 changes: 8 additions & 8 deletions src/diffusion/DiffusionScalarOp.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -221,10 +221,10 @@ DiffusionScalarOp::diffuse_scalar (Vector<MultiFab*> const& tracer,
#ifdef AMREX_USE_EB
if (m_eb_scal_solve_op) {
if ( m_incflo->m_has_mixedBC ) {
auto const robin = m_incflo->make_diffusion_robinBC_MFs(lev, phi[lev]);
auto const robin = m_incflo->make_robinBC_MFs(lev, &phi[lev]);

m_eb_scal_solve_op->setLevelBC(lev, &phi[lev],
robin[0].get(), robin[1].get(), robin[2].get());
&robin[0], &robin[1], &robin[2]);
}
else {
m_eb_scal_solve_op->setLevelBC(lev, &phi[lev]);
Expand Down Expand Up @@ -376,10 +376,10 @@ DiffusionScalarOp::diffuse_vel_components (Vector<MultiFab*> const& vel,
// m_eb_vel_solve_op->setPhiOnCentroid();

if ( m_incflo->m_has_mixedBC ) {
auto const robin = m_incflo->make_diffusion_robinBC_MFs(lev, phi[lev]);
auto const robin = m_incflo->make_robinBC_MFs(lev, &phi[lev]);

m_eb_vel_solve_op->setLevelBC(lev, &phi[lev],
robin[0].get(), robin[1].get(), robin[2].get());
&robin[0], &robin[1], &robin[2]);
}
else {
m_eb_vel_solve_op->setLevelBC(lev, &phi[lev]);
Expand Down Expand Up @@ -483,10 +483,10 @@ void DiffusionScalarOp::compute_laps (Vector<MultiFab*> const& a_laps,
m_eb_scal_apply_op->setBCoeffs(lev, GetArrOfConstPtrs(b), MLMG::Location::FaceCentroid);

if ( m_incflo->m_has_mixedBC ) {
auto const robin = m_incflo->make_diffusion_robinBC_MFs(lev, scalar_comp[lev]);
auto const robin = m_incflo->make_robinBC_MFs(lev, &scalar_comp[lev]);

m_eb_scal_apply_op->setLevelBC(lev, &scalar_comp[lev],
robin[0].get(), robin[1].get(), robin[2].get());
&robin[0], &robin[1], &robin[2]);
}
else {
m_eb_scal_apply_op->setLevelBC(lev, &scalar_comp[lev]);
Expand Down Expand Up @@ -605,10 +605,10 @@ void DiffusionScalarOp::compute_divtau (Vector<MultiFab*> const& a_divtau,
vel_single.emplace_back( vel[lev],amrex::make_alias,comp,1);

if ( m_incflo->m_has_mixedBC ) {
auto const robin = m_incflo->make_diffusion_robinBC_MFs(lev, vel_single[lev]);
auto const robin = m_incflo->make_robinBC_MFs(lev, &vel_single[lev]);

m_eb_vel_apply_op->setLevelBC(lev, &vel_single[lev],
robin[0].get(), robin[1].get(), robin[2].get());
&robin[0], &robin[1], &robin[2]);
}
else {
m_eb_vel_apply_op->setLevelBC(lev, &vel_single[lev]);
Expand Down
3 changes: 1 addition & 2 deletions src/incflo.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,8 +114,7 @@ public:
amrex::Gpu::DeviceVector<amrex::BCRec>& bcs,
std::string field);
amrex::iMultiFab make_nodalBC_mask (int lev);
amrex::Vector<std::unique_ptr<amrex::MultiFab>> make_diffusion_robinBC_MFs(int lev,
amrex::MultiFab& phi);
amrex::Vector<amrex::MultiFab> make_robinBC_MFs(int lev, amrex::MultiFab* phi = nullptr);
// void make_ccBC_mask (int lev, const amrex::BoxArray& ba,
// const amrex::DistributionMapping& dm);
// void make_nodalBC_mask (int lev, const amrex::BoxArray& ba,
Expand Down
4 changes: 2 additions & 2 deletions src/prob/prob_bc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,13 +88,13 @@ void incflo::prob_set_BC_MF (Orientation ori, Box const& bx,
}
}

// For MAC, the BC is on phi, so always 0 or 1
void incflo::prob_set_MAC_robinBCs (Orientation ori, Box const& bx,
Array4<Real> const& robin_a,
Array4<Real> const& robin_b,
Array4<Real> const& robin_f,
int lev)
{
// For MAC, the BC is on phi, so always 0 or 1
// Robin BC: a u + b du/dn = f -- inflow, Neumann a=0, b=1, f=0
// -- outflow, Dirichlet a=1, b=0, f=0

Expand Down Expand Up @@ -186,14 +186,14 @@ void incflo::prob_set_MAC_robinBCs (Orientation ori, Box const& bx,
}
}

// For diffusion, we also pass in the dirichlet bc
void incflo::prob_set_diffusion_robinBCs (Orientation ori, Box const& bx,
Array4<Real> const& robin_a,
Array4<Real> const& robin_b,
Array4<Real> const& robin_f,
Array4<Real const> const& bcval,
int lev)
{
// For diffusion, we also pass in the dirichlet bc (bcval)
// Robin BC: a u + b du/dn = f -- inflow, Dirichlet a=1, b=0, f=bcval
// -- outflow, Neumann a=0, b=1, f=0

Expand Down
Loading

0 comments on commit 9e82cf4

Please sign in to comment.