Skip to content

Commit

Permalink
WIP - MAC velocity looks reasonable. Needs hydro changes currently
Browse files Browse the repository at this point in the history 
at cgilet fork f2717ac
  • Loading branch information
@cgilet
cgilet committed Feb 12, 2024
1 parent 3f2dea0 commit 1529739
Show file tree
Hide file tree
Showing 4 changed files with 107 additions and 231 deletions.
167 changes: 86 additions & 81 deletions src/boundary_conditions/incflo_set_bcs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,41 +6,41 @@

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;
// 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);
}
}
}
}
// 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;
}
// return new_mask;
// }

iMultiFab
incflo::make_nodalBC_mask(int lev)
Expand All @@ -49,6 +49,9 @@ incflo::make_nodalBC_mask(int lev)
iMultiFab new_mask(amrex::convert(grids[lev], IntVect::TheNodeVector()), dmap[lev], 1, 0);
// Overset mask needs to be defined everywhere
new_mask = 1;
// NodalProj expects outflow regions to be indcatied with a 0. This is treated as a boolean
// mask, so any other regions can have any other int value
int outflow = 0;

Geometry const& gm = Geom(lev);
Box const& domain = gm.Domain();
Expand All @@ -66,10 +69,10 @@ incflo::make_nodalBC_mask(int lev)
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);
prob_set_BC_MF(olo, blo, mask_arr, lev, outflow);
}
if (bhi.ok()) {
prob_set_BC_MF(ohi, bhi, mask_arr, lev);
prob_set_BC_MF(ohi, bhi, mask_arr, lev, outflow);
}
}
}
Expand All @@ -78,54 +81,56 @@ incflo::make_nodalBC_mask(int lev)
return new_mask;
}

// void
// incflo::make_ccBC_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)
// {
// bool has_mixedBC = false;
// for (OrientationIter oit; oit; ++oit) {
// if (m_bc_type[oit()] == BC::mixed )
// {
// has_mixedBC = true;
// break;
// }
// }

// if (!has_mixedBC) { return; }


// // BC info stored in ghost cell.
// std::unique_ptr<iMultiFab> new_mask(new iMultiFab(ba, dm, 1, 1));
// *new_mask = 0;
std::unique_ptr<iMultiFab>
incflo::make_BC_MF(int lev, amrex::Gpu::DeviceVector<amrex::BCRec> bcs) //, int scomp, int ncomp)
{
int ncomp = bcs.size();
// BC info stored in ghost cells to avoid any ambiguity
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 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);

// 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);
// }
// }
// }
// }
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& mask_arr = BC_MF->array(mfi);
if (m_bc_type[olo] == BC::mixed) {
prob_set_BC_MF(olo, blo, mask_arr, lev, outflow);
} else {
amrex::ParallelFor(blo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
for (int n = 0; n < ncomp; n++) {
mask_arr(i,j,k,n) = bcs[n].lo(dir);
}
});
}
if (m_bc_type[ohi] == BC::mixed) {
prob_set_BC_MF(ohi, bhi, mask_arr, lev, outflow);
} else {
amrex::ParallelFor(bhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
for (int n = 0; n < ncomp; n++) {
mask_arr(i,j,k,n) = bcs[n].hi(dir);
}
});
}
}
}

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

// void
// incflo::make_nodalBC_mask(int lev,
Expand Down
148 changes: 7 additions & 141 deletions src/convection/incflo_compute_MAC_projected_velocities.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -233,56 +233,10 @@ incflo::compute_MAC_projected_velocities (
// full position dependent BCs, would need the MF to hold the BCType, which are enums,
// so ints...
//
// if (m_mixedBC_mask[lev]) {

// // Create MF to hold position-dependent BC info in ghost cells:
// // Low x-face is represented by indices ((i_lo-1, j_lo, k_lo)(i_lo-1, j_hi, k_hi))
// // High x-face is represented by indices ((i_hi+1, j_lo, k_lo)(i_hi+1, j_hi, k_hi))
// // Interior values don't get used
// iMultiFab bcs (grids[lev],dmap[lev],1,1,MFInfo(),Factory(lev));

// for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
// Orientation olo(dir,Orientation::low);
// Orientation ohi(dir,Orientation::high);

// Box dlo = amrex::adjCellLo(domain,dir,nghost);
// Box dhi = amrex::adjCellHi(domain,dir,nghost);
// #ifdef _OPENMP
// #pragma omp parallel if (Gpu::notInLaunchRegion())
// #endif
// for (MFIter mfi(bcs,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
// Box const& gbx = mfi.growntilebox();
// Box blo = gbx & dlo;
// Box bhi = gbx & dhi;
// Array4<Real> const& bcs_arr = bcs.array(mfi);
// Array4<int const> const& mask = m_mixedBC_mask[lev]->const_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
// if (blo.ok()) {
// // robin_b is the same as the mixedBC_mask, except with vals in
// // cell-centered ghosts rather than on BC face.
// // So, lo side i_cc->i_nd+1, hi side i_cc=i_nd
// Dim3 shift = IntVect::TheDimensionVector(dir).dim3();

// amrex::ParallelFor(blo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
// {
// b_arr(i,j,k) = mask(i+shift.x, j+shift.y, k+shift.z);
// //robin_a is the "opposite" of b; 0->1 and vice versa
// a_arr(i,j,k) = (mask(i+shift.x, j+shift.y, k+shift.z) == 1) ? 0. : 1.;
// f_arr(i,j,k) = 0.;
// });
// }
// if (bhi.ok()) {
// amrex::ParallelFor(bhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
// {
// b_arr(i,j,k) = mask(i,j,k);
// a_arr(i,j,k) = (mask(i,j,k) == 1) ? 0. : 1.;
// f_arr(i,j,k) = 0.;
// });
// }
// }
// }
std::unique_ptr<iMultiFab> BC_MF;
if (m_has_mixedBC) {
BC_MF = make_BC_MF(lev, m_bcrec_velocity_d);
}

// Predict normal velocity to faces -- note that the {u_mac, v_mac, w_mac}
// returned from this call are on face CENTROIDS
Expand All @@ -295,101 +249,13 @@ 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);
}

//fixme - fix up umac to confirm projection working as expected...
Box const& domain_box = geom[0].Domain();
int direction = 1;
int half_num_cells = domain_box.length(direction) / 2;
for (MFIter mfi(*u_mac[0],TilingIfNotGPU()); mfi.isValid(); ++mfi) {
Box const& bx = mfi.validbox();
Array4<Real > const& um = u_mac[0]->array(mfi);
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
if (i == domain_box.smallEnd(0))
{
if (j > half_num_cells) {
// Here we take the inflow BC specified in inputs file
um(i,j,k) = 1.;
}
}
else if (i == domain_box.bigEnd(0)+1)
{
if (j <= half_num_cells) {
// Here we take minus the inflow BC specified in inputs file
um(i,j,k) = -1.;
}
}
});
l_advection_type, PPM::default_limiter, BC_MF.get());
}

// if (m_mixedBC_mask[0]) {
// // Fix up the mixedBC faces which went through advection with FOEXTRAP (outflow) BCs
// // Here we overwrite the inflow portion with the Dirichlet BC
// // NOTE - there's a subtle difference in this versus what AMReX-Hydro does. See comment
// // in Utils/hydro_bcs_K.H . If weird gradients start arising at the inflow, we might
// // want to revisit this...
// // vel has filled ghost cells, this is needed for ExtrapVel
// for (int lev = 0; lev <= finest_level; ++lev)
// {
// int nghost = 1;
// Box const& domain = Geom(lev).growPeriodicDomain(nghost);
// for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
// Orientation olo(dir,Orientation::low);
// Orientation ohi(dir,Orientation::high);

// Box dlo = amrex::adjCellLo(domain,dir,nghost);
// Box dhi = amrex::adjCellHi(domain,dir,nghost);
// #ifdef _OPENMP
// #pragma omp parallel if (Gpu::notInLaunchRegion())
// #endif
// for (MFIter mfi(*vel[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) {
// Box const& gbx = amrex::grow(mfi.validbox(),1);
// Box blo = gbx & dlo;
// Box bhi = gbx & dhi;

// MultiFab* mac;
// if (dir == 0) {
// mac = u_mac[lev]; }
// else if (dir == 1) {
// mac = v_mac[lev]; }
// #if AMREX_SPACEDIM > 2
// else {
// mac = w_mac[lev];
// }
// #endif
// Array4<Real > const& mac_arr = mac->array(mfi);
// Array4<Real const> const& vel_arr = vel[lev]->array(mfi, dir);
// Array4<int const> const& mask = m_mixedBC_mask[lev]->const_array(mfi);

// if (blo.ok()) {
// Dim3 shift = IntVect::TheDimensionVector(dir).dim3();

// amrex::ParallelFor(blo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
// {
// if (mask(i+shift.x, j+shift.y, k+shift.z) == 1) // inflow
// {
// mac_arr(i+shift.x, j+shift.y, k+shift.z) = vel_arr(i,j,k);
// }
// });
// }
// if (bhi.ok()) {
// amrex::ParallelFor(bhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
// {
// if (mask(i,j,k) == 1) {
// mac_arr(i,j,k) = vel_arr(i,j,k);
// }
// });
// }
// }
// }
// }
// }
//fixme
//fixmeamrex::
VisMF::Write(*u_mac[0],"umac");
VisMF::Write(*v_mac[0],"vmac");

//
Vector<Array<MultiFab*,AMREX_SPACEDIM> > mac_vec(finest_level+1);
for (int lev=0; lev <= finest_level; ++lev)
{
Expand Down
5 changes: 3 additions & 2 deletions src/incflo.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,7 @@ public:
//
///////////////////////////////////////////////////////////////////////////

amrex::iMultiFab make_ccBC_mask (int lev);
std::unique_ptr<amrex::iMultiFab> make_BC_MF (int lev, amrex::Gpu::DeviceVector<amrex::BCRec> bcs);
amrex::iMultiFab make_nodalBC_mask (int lev);
// void make_ccBC_mask (int lev, const amrex::BoxArray& ba,
// const amrex::DistributionMapping& dm);
Expand Down Expand Up @@ -218,7 +218,8 @@ public:

void prob_init_fluid (int lev);
void prob_set_BC_MF (amrex::Orientation ori, amrex::Box const& bx,
amrex::Array4<int> const& mask, int lev);
amrex::Array4<int> const& mask, int lev,
int outflow_val);
void prob_set_MAC_robinBCs (amrex::Orientation ori, amrex::Box const& bx,
amrex::Array4<amrex::Real> const& robin_a,
amrex::Array4<amrex::Real> const& robin_b,
Expand Down
Loading

0 comments on commit 1529739

Please sign in to comment.