Skip to content

Commit

Permalink
Add MLMG solver for EB
Browse files Browse the repository at this point in the history
  • Loading branch information
@asalmgren
asalmgren committed Nov 21, 2024
1 parent b301795 commit 14d9f98
Show file tree
Hide file tree
Showing 10 changed files with 256 additions and 20 deletions.
13 changes: 8 additions & 5 deletions Exec/DevTests/EB_Test/ERF_prob.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,8 +67,8 @@ Problem::init_custom_pert(
Array4<Real > const& z_vel_pert,
Array4<Real > const& /*r_hse*/,
Array4<Real > const& /*p_hse*/,
Array4<Real const> const& z_nd,
Array4<Real const> const& z_cc,
Array4<Real const> const& /*z_nd*/,
Array4<Real const> const& /*z_cc*/,
GeometryData const& geomdata,
Array4<Real const> const& /*mf_m*/,
Array4<Real const> const& /*mf_u*/,
Expand Down Expand Up @@ -148,7 +148,7 @@ Problem::init_custom_pert(

// User function parameters
Real a = 0.5;
Real num = 8 * a * a * a;
Real num = 8.11 * a * a * a;
Real xcen = 0.5 * (prob_lo[0] + prob_hi[0]);

// Grown box with no z range
Expand All @@ -157,6 +157,9 @@ Problem::init_custom_pert(

amrex::Array4<amrex::Real> const& z_arr = z_phys_nd.array();

Real x_in = (-xcen);
Real height_at_inflow = num / (x_in * x_in + 4.0 * a * a);

ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int)
{
// Clip indices for ghost-cells
Expand All @@ -166,9 +169,9 @@ Problem::init_custom_pert(
Real x = (ii * dx[0] - xcen);

// WoA Hill in x-direction
Real height = num / (x*x + 4 * a * a);
Real height = num / (x*x + 4.0 * a * a);

// Populate terrain height
z_arr(i,j,k0) = height;
z_arr(i,j,k0) = height - height_at_inflow;
});
}
1 change: 0 additions & 1 deletion Exec/DevTests/EB_Test/GNUmakefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,6 @@ USE_SYCL = FALSE

# Debugging
DEBUG = FALSE
DEBUG = TRUE

TEST = TRUE
USE_ASSERTION = TRUE
Expand Down
16 changes: 10 additions & 6 deletions Exec/DevTests/EB_Test/inputs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,29 +5,33 @@ max_step = 0
amr.max_grid_size = 256 256 256

eb2.geometry = terrain
erf.anelastic = 1
erf.mg_v = 2

#erf.project_initial_velocity=true;

amrex.fpe_trap_invalid = 1

fabarray.mfiter_tile_size = 1024 1024 1024

# PROBLEM SIZE & GEOMETRY
geometry.prob_lo = 0. 0. 0.
geometry.prob_hi = 10. 1. 4.
geometry.prob_hi = 16. 4. 16.

amr.n_cell = 128 8 64
amr.n_cell = 64 16 64

geometry.is_periodic = 0 1 0

xlo.type = "Inflow"
xhi.type = "Outflow"
xlo.type = Inflow
xhi.type = Outflow

xlo.velocity = 1. 0. 0.
xlo.density = 1.16
xlo.theta = 300.
xlo.scalar = 0.

zlo.type = "SlipWall"
zhi.type = "SlipWall"
zlo.type = SlipWall
zhi.type = SlipWall

# TIME STEP CONTROL
erf.substepping_type = None
Expand Down
7 changes: 6 additions & 1 deletion Source/EB/ERF_initEB.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,12 @@ void ERF::MakeEBGeometry()
* *
******************************************************************************/

int max_coarsening_level = 0;
int max_coarsening_level;
if (solverChoice.anelastic[0] == 1) {
max_coarsening_level = 100;
} else {
max_coarsening_level = 0;
}

if(geom_type == "cylinder") {
amrex::Print() << "\n Building cylinder geometry." << std::endl;
Expand Down
11 changes: 9 additions & 2 deletions Source/ERF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -176,6 +176,10 @@ public:
void solve_with_fft (int lev, amrex::MultiFab& rhs, amrex::MultiFab& p,
amrex::Array<amrex::MultiFab,AMREX_SPACEDIM>& fluxes);
#endif
#ifdef ERF_USE_EB
void solve_with_EB_mlmg (int lev, amrex::Vector<amrex::MultiFab>& rhs, amrex::Vector<amrex::MultiFab>& p,
amrex::Vector<amrex::Array<amrex::MultiFab,AMREX_SPACEDIM>>& fluxes);
#endif

void solve_with_gmres (int lev, amrex::Vector<amrex::MultiFab>& rhs, amrex::Vector<amrex::MultiFab>& p,
amrex::Vector<amrex::Array<amrex::MultiFab,AMREX_SPACEDIM>>& fluxes);
Expand Down Expand Up @@ -1342,11 +1346,12 @@ private:
//! The filename of the ith samplelinelog file.
[[nodiscard]] std::string SampleLineLogName (int i) const noexcept { return samplelinelogname[i]; }

amrex::Vector<std::unique_ptr<amrex::FabFactory<amrex::FArrayBox> > > m_factory;
#ifdef ERF_USE_EB

amrex::Vector<std::unique_ptr<amrex::EBFArrayBoxFactory> > m_factory;

[[nodiscard]] amrex::FabFactory<amrex::FArrayBox> const&
Factory (int lev) const noexcept { return *m_factory[lev]; }
#ifdef ERF_USE_EB
[[nodiscard]] amrex::EBFArrayBoxFactory const&
EBFactory (int lev) const noexcept {
return static_cast<amrex::EBFArrayBoxFactory const&>(*m_factory[lev]);
Expand All @@ -1361,6 +1366,8 @@ private:

[[nodiscard]] static int nghost_eb_full ()
{ return 4; }
#else
amrex::Vector<std::unique_ptr<amrex::FabFactory<amrex::FArrayBox> > > m_factory;
#endif

#ifdef ERF_USE_FFT
Expand Down
4 changes: 1 addition & 3 deletions Source/ERF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -333,12 +333,10 @@ ERF::ERF_shared ()
ref_ratio[lev][0] << " " << ref_ratio[lev][1] << " " << ref_ratio[lev][2] << std::endl;
}

// We define m_factory even with no EB
// We will create each of these in MakeNewLevel.../RemakeLevel
m_factory.resize(max_level+1);

#ifdef ERF_USE_EB
// We will create each of these in MakeNewLevel.../RemakeLevel

// This is needed before initializing level MultiFabs
MakeEBGeometry();
#endif
Expand Down
27 changes: 25 additions & 2 deletions Source/TimeIntegration/ERF_TI_no_substep_fun.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,24 +106,47 @@
});

} else { // Fixed or no terrain
#ifdef ERF_USE_EB
const Array4<const Real>& dJ_old = detJ_cc[level]->const_array(mfi);
#endif
ParallelFor(bx, ncomp_fast[IntVars::cons],
[=] AMREX_GPU_DEVICE (int i, int j, int k, int nn) {
const int n = scomp_fast[IntVars::cons] + nn;
ssum[IntVars::cons](i,j,k,n) = sold[IntVars::cons](i,j,k,n) + slow_dt *
( fslow[IntVars::cons](i,j,k,n) );
#ifdef ERF_USE_EB
if (dJ_old(i,j,k) > 0.0) {
#endif
ssum[IntVars::cons](i,j,k,n) = sold[IntVars::cons](i,j,k,n) + slow_dt *
( fslow[IntVars::cons](i,j,k,n) );
#ifdef ERF_USE_EB
}
#endif
});

// Commenting out the update is a total HACK while developing the EB capability
ParallelFor(tbx, tby, tbz,
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
#ifdef ERF_USE_EB
ssum[IntVars::xmom](i,j,k) = sold[IntVars::xmom](i,j,k);
#else
ssum[IntVars::xmom](i,j,k) = sold[IntVars::xmom](i,j,k)
+ slow_dt * fslow[IntVars::xmom](i,j,k);
#endif
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
#ifdef ERF_USE_EB
ssum[IntVars::ymom](i,j,k) = sold[IntVars::ymom](i,j,k);
#else
ssum[IntVars::ymom](i,j,k) = sold[IntVars::ymom](i,j,k)
+ slow_dt * fslow[IntVars::ymom](i,j,k);
#endif
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
#ifdef ERF_USE_EB
ssum[IntVars::zmom](i,j,k) = sold[IntVars::zmom](i,j,k);
#else
ssum[IntVars::zmom](i,j,k) = sold[IntVars::zmom](i,j,k)
+ slow_dt * fslow[IntVars::zmom](i,j,k);
#endif
});
}
} // mfi
Expand Down
20 changes: 20 additions & 0 deletions Source/Utils/ERF_PoissonSolve.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,8 +14,10 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
bool l_use_terrain = SolverChoice::terrain_type != TerrainType::None;
bool use_gmres = (l_use_terrain && !SolverChoice::terrain_is_flat);

#ifndef ERF_USE_EB
auto const dom_lo = lbound(geom[lev].Domain());
auto const dom_hi = ubound(geom[lev].Domain());
#endif

// Make sure the solver only sees the levels over which we are solving
Vector<BoxArray> ba_tmp; ba_tmp.push_back(mom_mf[Vars::cons].boxArray());
Expand All @@ -28,12 +30,21 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
Vector<MultiFab> phi;
Vector<Array<MultiFab,AMREX_SPACEDIM> > fluxes;

#ifdef ERF_USE_EB
rhs.resize(1); rhs[0].define(ba_tmp[0], dm_tmp[0], 1, 0, MFInfo(), Factory(lev));
phi.resize(1); phi[0].define(ba_tmp[0], dm_tmp[0], 1, 1, MFInfo(), Factory(lev));
#else
rhs.resize(1); rhs[0].define(ba_tmp[0], dm_tmp[0], 1, 0);
phi.resize(1); phi[0].define(ba_tmp[0], dm_tmp[0], 1, 1);
#endif

fluxes.resize(1);
for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
#ifdef ERF_USE_EB
fluxes[0][idim].define(convert(ba_tmp[0], IntVect::TheDimensionVector(idim)), dm_tmp[0], 1, 0, MFInfo(), Factory(lev));
#else
fluxes[0][idim].define(convert(ba_tmp[0], IntVect::TheDimensionVector(idim)), dm_tmp[0], 1, 0);
#endif
}

auto dxInv = geom[lev].InvCellSizeArray();
Expand All @@ -49,6 +60,10 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
// Compute divergence which will form RHS
// Note that we replace "rho0w" with the contravariant momentum, Omega
// ****************************************************************************
#ifdef ERF_USE_EB
bool already_on_centroids = true;
EB_computeDivergence(rhs[0], rho0_u_const, geom_tmp[0], already_on_centroids);
#else
if (l_use_terrain)
{
for ( MFIter mfi(rhs[0],TilingIfNotGPU()); mfi.isValid(); ++mfi)
Expand Down Expand Up @@ -96,6 +111,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
computeDivergence(rhs[0], rho0_u_const, geom_tmp[0]);

}
#endif

Real rhsnorm = rhs[0].norm0();

Expand All @@ -115,6 +131,9 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult
// ****************************************************************************
// Choose the solver and solve
// ****************************************************************************
#ifdef ERF_USE_EB
solve_with_EB_mlmg(lev, rhs, phi, fluxes);
#else
#ifdef ERF_USE_FFT
if (use_fft) {

Expand All @@ -130,6 +149,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector<MultiFab>& mom_mf, Mult

solve_with_mlmg(lev, rhs, phi, fluxes);
}
#endif

// ****************************************************************************
// Subtract dt grad(phi) from the momenta (rho0u, rho0v, Omega)
Expand Down
Loading

0 comments on commit 14d9f98

Please sign in to comment.