Skip to content

Commit

Permalink
Make scalar diffusion work with mixed BC. Plus various mixed bc
Browse files Browse the repository at this point in the history 
fixes.
Use with cgilet hydro fdbb31e
  • Loading branch information
@cgilet
cgilet committed Feb 14, 2024
1 parent e3f9635 commit 293c196
Show file tree
Hide file tree
Showing 9 changed files with 376 additions and 114 deletions.
65 changes: 62 additions & 3 deletions src/boundary_conditions/incflo_set_bcs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,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, outflow, "velocity");
prob_set_BC_MF(olo, blo, mask_arr, lev, outflow, "projection");
}
if (bhi.ok()) {
prob_set_BC_MF(ohi, bhi, mask_arr, lev, outflow, "velocity");
prob_set_BC_MF(ohi, bhi, mask_arr, lev, outflow, "projection");
}
}
}
Expand All @@ -82,7 +82,7 @@ incflo::make_nodalBC_mask(int lev)
}

std::unique_ptr<iMultiFab>
incflo::make_BC_MF(int lev, amrex::Gpu::DeviceVector<amrex::BCRec> bcs,
incflo::make_BC_MF(int lev, amrex::Gpu::DeviceVector<amrex::BCRec>& bcs,
std::string field)
{
int ncomp = bcs.size();
Expand Down Expand Up @@ -133,6 +133,65 @@ 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)
{
int nghost = 1;

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

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;

// 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 - don't think 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);
}
}
}
}

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...
Expand Down
37 changes: 3 additions & 34 deletions src/convection/incflo_compute_MAC_projected_velocities.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,7 +115,7 @@ incflo::compute_MAC_projected_velocities (
// 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 and that's it
// 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;
Expand All @@ -132,13 +132,8 @@ incflo::compute_MAC_projected_velocities (
robin_b = 0;
robin_f = 0;

// I don't believe interior values get used,
// only ghost cells of robin BC arrays are used in MLMG
// bc in ghost cells that are outside the domain.
// What does MLMG expect for periodic? Should have been taken
// care of in creating the mask
//amrex::Geometry const& gm = Geom(lev);
//Box const& domain = gm.growPeriodicDomain(nghost);
Box const& domain = Geom(lev).Domain();
for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
Orientation olo(dir,Orientation::low);
Expand All @@ -158,45 +153,19 @@ incflo::compute_MAC_projected_velocities (
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<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
// Only need to fill Robin BC sides, MLMG will check for Robin BC first
// Only need to fill Robin BC sides, MLMG will check for Robin BC first
if (blo.ok()) {
// robin_b is the same as the nodalBC_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();

prob_set_MAC_robinBCs(olo, blo, a_arr, b_arr, f_arr, lev);

// 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()) {

prob_set_MAC_robinBCs(ohi, bhi, a_arr, b_arr, f_arr, lev);

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

VisMF::Write(robin_a, "ra");
VisMF::Write(robin_b, "rb");
VisMF::Write(robin_f, "rf");
macproj->setLevelBC(lev, nullptr, &robin_a, &robin_b, &robin_f);
}
}
Expand Down Expand Up @@ -239,7 +208,7 @@ incflo::compute_MAC_projected_velocities (
}
//fixme
VisMF::Write(*vel[0],"vin");
amrex::Write(*BC_MF,"bcmf");
//amrex::Write(*BC_MF,"bcmf");
//

// Predict normal velocity to faces -- note that the {u_mac, v_mac, w_mac}
Expand Down
11 changes: 7 additions & 4 deletions src/convection/incflo_compute_advection_term.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -340,7 +340,8 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
int face_comp = 0;
int ncomp = AMREX_SPACEDIM;
bool is_velocity = true;
Array4<int const> const& velbc_arr = (*velBC_MF).const_array(mfi);
Array4<int const> const& velbc_arr = velBC_MF ? (*velBC_MF).const_array(mfi)
: Array4<int const>{};
HydroUtils::ComputeFluxesOnBoxFromState( bx, ncomp, mfi,
(m_advect_momentum) ? rhovel[lev].array(mfi) : vel[lev]->const_array(mfi),
AMREX_D_DECL(flux_x[lev].array(mfi,face_comp),
Expand Down Expand Up @@ -380,7 +381,8 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
face_comp = AMREX_SPACEDIM;
ncomp = 1;
is_velocity = false;
Array4<int const> const& densbc_arr = (*densBC_MF).const_array(mfi);
Array4<int const> const& densbc_arr = densBC_MF ? (*densBC_MF).const_array(mfi)
: Array4<int const>{};
HydroUtils::ComputeFluxesOnBoxFromState( bx, ncomp, mfi,
density[lev]->const_array(mfi),
AMREX_D_DECL(flux_x[lev].array(mfi,face_comp),
Expand Down Expand Up @@ -434,7 +436,8 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
face_comp = AMREX_SPACEDIM+1;
ncomp = m_ntrac;
is_velocity = false;
Array4<int const> const& tracbc_arr = (*tracBC_MF).const_array(mfi);
Array4<int const> const& tracbc_arr = tracBC_MF ? (*tracBC_MF).const_array(mfi)
: Array4<int const>{};
HydroUtils::ComputeFluxesOnBoxFromState( bx, ncomp, mfi, ro_trac,
AMREX_D_DECL(flux_x[lev].array(mfi,face_comp),
flux_y[lev].array(mfi,face_comp),
Expand Down Expand Up @@ -705,5 +708,5 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
} // lev
//fixme
VisMF::Write(*conv_u[0],"conv");
//
VisMF::Write(*conv_r[0],"conr");
}
30 changes: 27 additions & 3 deletions src/diffusion/DiffusionScalarOp.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -366,8 +366,19 @@ DiffusionScalarOp::diffuse_vel_components (Vector<MultiFab*> const& vel,
// For when we use the stencil for centroid values
// m_eb_vel_solve_op->setPhiOnCentroid();

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

// VisMF::Write(robin_a, "dra");
// VisMF::Write(robin_b, "drb");
// VisMF::Write(robin_f, "drf");

m_eb_vel_solve_op->setLevelBC(lev, &phi[lev],
robin[0].get(), robin[1].get(), robin[2].get());
}
else {
m_eb_vel_solve_op->setLevelBC(lev, &phi[lev]);
}
} else
#endif
{
Expand Down Expand Up @@ -577,8 +588,21 @@ void DiffusionScalarOp::compute_divtau (Vector<MultiFab*> const& a_divtau,

for (int lev = 0; lev <= finest_level; ++lev) {
divtau_single.emplace_back(divtau_tmp[lev],amrex::make_alias,comp,1);
vel_single.emplace_back( vel[lev],amrex::make_alias,comp,1);
m_eb_vel_apply_op->setLevelBC(lev, &vel_single[lev]);
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]);

// VisMF::Write(robin_a, "dra");
// VisMF::Write(robin_b, "drb");
// VisMF::Write(robin_f, "drf");

m_eb_vel_apply_op->setLevelBC(lev, &vel_single[lev],
robin[0].get(), robin[1].get(), robin[2].get());
}
else {
m_eb_vel_apply_op->setLevelBC(lev, &vel_single[lev]);
}

Array<MultiFab,AMREX_SPACEDIM> b =
m_incflo->average_scalar_eta_to_faces(lev, eta_comp, *a_eta[lev]);
Expand Down
12 changes: 12 additions & 0 deletions src/diffusion/incflo_diffusion.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -193,6 +193,13 @@ incflo::get_diffuse_velocity_bc (Orientation::Side side, int comp) const noexcep
r[dir][dir] = LinOpBCType::Dirichlet;
break;
}
case BC::mixed:
{
AMREX_D_TERM(r[0][dir] = LinOpBCType::Robin;,
r[1][dir] = LinOpBCType::Robin;,
r[2][dir] = LinOpBCType::Robin;);
break;
}
default:
amrex::Abort("get_diffuse_tensor_bc: undefined BC type");
};
Expand Down Expand Up @@ -237,6 +244,11 @@ incflo::get_diffuse_scalar_bc (Orientation::Side side) const noexcept
r[dir] = LinOpBCType::Dirichlet;
break;
}
case BC::mixed:
{
r[dir] = LinOpBCType::Robin;
break;
}
default:
amrex::Abort("get_diffuse_scalar_bc: undefined BC type");
};
Expand Down
10 changes: 9 additions & 1 deletion src/incflo.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,9 +111,11 @@ public:
///////////////////////////////////////////////////////////////////////////

std::unique_ptr<amrex::iMultiFab> make_BC_MF (int lev,
amrex::Gpu::DeviceVector<amrex::BCRec> bcs,
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);
// 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 Expand Up @@ -227,6 +229,12 @@ public:
amrex::Array4<amrex::Real> const& robin_b,
amrex::Array4<amrex::Real> const& robin_f,
int lev);
void prob_set_diffusion_robinBCs (amrex::Orientation ori, amrex::Box const& bx,
amrex::Array4<amrex::Real> const& robin_a,
amrex::Array4<amrex::Real> const& robin_b,
amrex::Array4<amrex::Real> const& robin_f,
amrex::Array4<amrex::Real const> const& bcval,
int lev);

#include "incflo_prob_I.H"
#include "incflo_prob_usr_I.H"
Expand Down
21 changes: 13 additions & 8 deletions src/incflo_apply_predictor.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -140,21 +140,26 @@ void incflo::ApplyPredictor (bool incremental_projection)
if (need_divtau() || use_tensor_correction )
{
// FIXME - deal with viscosity later...
auto divtau = get_divtau_old();
for (int lev = 0; lev <= finest_level; ++lev) {
*divtau[lev] = 0.;
}
// compute_divtau(get_divtau_old(),get_velocity_old_const(),
// get_density_old_const(),GetVecOfConstPtrs(vel_eta));
// auto divtau = get_divtau_old();
// for (int lev = 0; lev <= finest_level; ++lev) {
// *divtau[lev] = 0.;
// }
compute_divtau(get_divtau_old(),get_velocity_old_const(),
get_density_old_const(),GetVecOfConstPtrs(vel_eta));
}

// *************************************************************************************
// Compute explicit diffusive terms
// *************************************************************************************
if (m_advect_tracer && need_divtau())
{
compute_laps(get_laps_old(), get_tracer_old_const(), get_density_old_const(),
GetVecOfConstPtrs(tra_eta));
// FIXME - deal with diffusion later...
auto laps = get_laps_old();
for (int lev = 0; lev <= finest_level; ++lev) {
*laps[lev] = 0.;
}
// compute_laps(get_laps_old(), get_tracer_old_const(), get_density_old_const(),
// GetVecOfConstPtrs(tra_eta));
}

// **********************************************************************************************
Expand Down
Loading

0 comments on commit 293c196

Please sign in to comment.