Skip to content

Commit

Permalink
MLEBNodeFDLaplacian: Make it work with AMREX_USE_EB but with no EB
Browse files Browse the repository at this point in the history 
This will allow the use of the linear solver when AMReX is built with EB
support at compile time but no EB is built at runtime.
  • Loading branch information
@WeiqunZhang
WeiqunZhang committed Aug 15, 2024
1 parent 8535675 commit 0101d34
Showing 1 changed file with 54 additions and 32 deletions.
86 changes: 54 additions & 32 deletions Src/LinearSolvers/MLMG/AMReX_MLEBNodeFDLaplacian.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
#include <AMReX_MultiFabUtil.H>

#ifdef AMREX_USE_EB
#include <AMReX_EB2.H>
#include <AMReX_EBMultiFabUtil.H>
#endif

Expand Down Expand Up @@ -151,10 +152,14 @@ MLEBNodeFDLaplacian::define (const Vector<Geometry>& a_geom,
std::unique_ptr<FabFactory<FArrayBox> >
MLEBNodeFDLaplacian::makeFactory (int amrlev, int mglev) const
{
return makeEBFabFactory(m_geom[amrlev][mglev],
m_grids[amrlev][mglev],
m_dmap[amrlev][mglev],
{1,1,1}, EBSupport::full);
if (EB2::TopIndexSpaceIfPresent()) {
return makeEBFabFactory(m_geom[amrlev][mglev],
m_grids[amrlev][mglev],
m_dmap[amrlev][mglev],
{1,1,1}, EBSupport::full);
} else {
return MLNodeLinOp::makeFactory(amrlev, mglev);
}
}
#endif

Expand Down Expand Up @@ -264,17 +269,19 @@ MLEBNodeFDLaplacian::prepareForSolve ()
for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev) {
for (int mglev = 0; mglev < m_num_mg_levels[amrlev]; ++mglev) {
const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(m_factory[amrlev][mglev].get());
auto const& levset_mf = factory->getLevelSet();
auto const& levset_ar = levset_mf.const_arrays();
auto& dmask_mf = *m_dirichlet_mask[amrlev][mglev];
auto const& dmask_ar = dmask_mf.arrays();
amrex::ParallelFor(dmask_mf,
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
{
if (levset_ar[box_no](i,j,k) >= Real(0.0)) {
dmask_ar[box_no](i,j,k) = -1;
}
});
if (factory) {
auto const& levset_mf = factory->getLevelSet();
auto const& levset_ar = levset_mf.const_arrays();
auto& dmask_mf = *m_dirichlet_mask[amrlev][mglev];
auto const& dmask_ar = dmask_mf.arrays();
amrex::ParallelFor(dmask_mf,
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
{
if (levset_ar[box_no](i,j,k) >= Real(0.0)) {
dmask_ar[box_no](i,j,k) = -1;
}
});
}
}
}
#endif
Expand Down Expand Up @@ -359,8 +366,10 @@ MLEBNodeFDLaplacian::prepareForSolve ()
void
MLEBNodeFDLaplacian::scaleRHS (int amrlev, MultiFab& rhs) const
{
auto const& dmask = *m_dirichlet_mask[amrlev][0];
const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(m_factory[amrlev][0].get());
if (!factory) { return; }

auto const& dmask = *m_dirichlet_mask[amrlev][0];
auto const& edgecent = factory->getEdgeCent();

#ifdef AMREX_USE_OMP
Expand Down Expand Up @@ -407,9 +416,14 @@ MLEBNodeFDLaplacian::Fapply (int amrlev, int mglev, MultiFab& out, const MultiFa
#ifdef AMREX_USE_EB
const auto phieb = (m_in_solution_mode) ? m_s_phi_eb : Real(0.0);
const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(m_factory[amrlev][mglev].get());
auto const& edgecent = factory->getEdgeCent();
auto const& levset_mf = factory->getLevelSet();
auto const& volfrac = factory->getVolFrac();
Array<const MultiCutFab*,AMREX_SPACEDIM> edgecent;
MultiFab const* levset_mf = nullptr;
MultiFab const* volfrac = nullptr;
if (factory) {
edgecent = factory->getEdgeCent();
levset_mf = &(factory->getLevelSet());
volfrac = &(factory->getVolFrac());
}
#endif

#ifdef AMREX_USE_OMP
Expand All @@ -422,12 +436,12 @@ MLEBNodeFDLaplacian::Fapply (int amrlev, int mglev, MultiFab& out, const MultiFa
Array4<Real> const& yarr = out.array(mfi);
Array4<int const> const& dmarr = dmask.const_array(mfi);
#ifdef AMREX_USE_EB
bool cutfab = edgecent[0]->ok(mfi);
bool cutfab = factory && edgecent[0]->ok(mfi);
if (cutfab) {
AMREX_D_TERM(Array4<Real const> const& ecx = edgecent[0]->const_array(mfi);,
Array4<Real const> const& ecy = edgecent[1]->const_array(mfi);,
Array4<Real const> const& ecz = edgecent[2]->const_array(mfi));
auto const& levset = levset_mf.const_array(mfi);
auto const& levset = levset_mf->const_array(mfi);
if (phieb == std::numeric_limits<Real>::lowest()) {
auto const& phiebarr = m_phi_eb[amrlev].const_array(mfi);
#if (AMREX_SPACEDIM == 2)
Expand All @@ -441,7 +455,7 @@ MLEBNodeFDLaplacian::Fapply (int amrlev, int mglev, MultiFab& out, const MultiFa
#endif
if (m_has_sigma_mf) {
auto const& sigarr = m_sigma_mf[amrlev][mglev]->const_array(mfi);
auto const& vfrc = volfrac.const_array(mfi);
auto const& vfrc = volfrac->const_array(mfi);
AMREX_HOST_DEVICE_FOR_3D(box, i, j, k,
{
mlebndfdlap_sig_adotx_eb(i,j,k,yarr,xarr,levset,dmarr,AMREX_D_DECL(ecx,ecy,ecz),
Expand All @@ -466,7 +480,7 @@ MLEBNodeFDLaplacian::Fapply (int amrlev, int mglev, MultiFab& out, const MultiFa
#endif
if (m_has_sigma_mf) {
auto const& sigarr = m_sigma_mf[amrlev][mglev]->const_array(mfi);
auto const& vfrc = volfrac.const_array(mfi);
auto const& vfrc = volfrac->const_array(mfi);
AMREX_HOST_DEVICE_FOR_3D(box, i, j, k,
{
mlebndfdlap_sig_adotx_eb(i,j,k,yarr,xarr,levset,dmarr,AMREX_D_DECL(ecx,ecy,ecz),
Expand Down Expand Up @@ -533,9 +547,14 @@ MLEBNodeFDLaplacian::Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiF

#ifdef AMREX_USE_EB
const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(m_factory[amrlev][mglev].get());
auto const& edgecent = factory->getEdgeCent();
auto const& levset_mf = factory->getLevelSet();
auto const& volfrac = factory->getVolFrac();
Array<const MultiCutFab*,AMREX_SPACEDIM> edgecent;
MultiFab const* levset_mf = nullptr;
MultiFab const* volfrac = nullptr;
if (factory) {
edgecent = factory->getEdgeCent();
levset_mf = &(factory->getLevelSet());
volfrac = &(factory->getVolFrac());
}
#endif

#ifdef AMREX_USE_OMP
Expand All @@ -548,12 +567,12 @@ MLEBNodeFDLaplacian::Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiF
Array4<Real const> const& rhsarr = rhs.const_array(mfi);
Array4<int const> const& dmskarr = dmask.const_array(mfi);
#ifdef AMREX_USE_EB
bool cutfab = edgecent[0]->ok(mfi);
bool cutfab = factory && edgecent[0]->ok(mfi);
if (cutfab) {
AMREX_D_TERM(Array4<Real const> const& ecx = edgecent[0]->const_array(mfi);,
Array4<Real const> const& ecy = edgecent[1]->const_array(mfi);,
Array4<Real const> const& ecz = edgecent[2]->const_array(mfi));
auto const& levset = levset_mf.const_array(mfi);
auto const& levset = levset_mf->const_array(mfi);
#if (AMREX_SPACEDIM == 2)
if (m_rz) {
AMREX_HOST_DEVICE_FOR_3D(box, i, j, k,
Expand All @@ -565,7 +584,7 @@ MLEBNodeFDLaplacian::Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiF
#endif
if (m_has_sigma_mf) {
auto const& sigarr = m_sigma_mf[amrlev][mglev]->const_array(mfi);
auto const& vfrc = volfrac.const_array(mfi);
auto const& vfrc = volfrac->const_array(mfi);
AMREX_HOST_DEVICE_FOR_3D(box, i, j, k,
{
mlebndfdlap_sig_gsrb_eb(i,j,k,solarr,rhsarr,levset,dmskarr,AMREX_D_DECL(ecx,ecy,ecz),
Expand Down Expand Up @@ -641,8 +660,10 @@ MLEBNodeFDLaplacian::compGrad (int amrlev, const Array<MultiFab*,AMREX_SPACEDIM>
auto const& dmask = *m_dirichlet_mask[amrlev][mglev];
const auto phieb = m_s_phi_eb;
const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(m_factory[amrlev][mglev].get());
AMREX_ASSERT(factory);
auto const& edgecent = factory->getEdgeCent();
Array<const MultiCutFab*,AMREX_SPACEDIM> edgecent;
if (factory) {
edgecent = factory->getEdgeCent();
}
#endif

#ifdef AMREX_USE_OMP
Expand All @@ -659,7 +680,7 @@ MLEBNodeFDLaplacian::compGrad (int amrlev, const Array<MultiFab*,AMREX_SPACEDIM>
Array4<Real> const& gpz = grad[2]->array(mfi);)
#ifdef AMREX_USE_EB
Array4<int const> const& dmarr = dmask.const_array(mfi);
bool cutfab = edgecent[0]->ok(mfi);
bool cutfab = factory && edgecent[0]->ok(mfi);
AMREX_D_TERM(Array4<Real const> const& ecx
= cutfab ? edgecent[0]->const_array(mfi) : Array4<Real const>{};,
Array4<Real const> const& ecy
Expand Down Expand Up @@ -741,6 +762,7 @@ MLEBNodeFDLaplacian::postSolve (Vector<MultiFab>& sol) const
for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev) {
const auto phieb = m_s_phi_eb;
const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(m_factory[amrlev][0].get());
if (!factory) { return; }
auto const& levset_mf = factory->getLevelSet();
auto const& levset_ar = levset_mf.const_arrays();
MultiFab& mf = sol[amrlev];
Expand Down

0 comments on commit 0101d34

Please sign in to comment.