From 14d9f9871a1c2dae37002f5d2f95c20cf152b0a2 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Thu, 21 Nov 2024 14:08:30 -0800 Subject: [PATCH] Add MLMG solver for EB --- Exec/DevTests/EB_Test/ERF_prob.cpp | 13 +- Exec/DevTests/EB_Test/GNUmakefile | 1 - Exec/DevTests/EB_Test/inputs | 16 +- Source/EB/ERF_initEB.cpp | 7 +- Source/ERF.H | 11 +- Source/ERF.cpp | 4 +- .../TimeIntegration/ERF_TI_no_substep_fun.H | 27 ++- Source/Utils/ERF_PoissonSolve.cpp | 20 +++ Source/Utils/ERF_solve_with_EB_mlmg.cpp | 170 ++++++++++++++++++ Source/Utils/Make.package | 7 + 10 files changed, 256 insertions(+), 20 deletions(-) create mode 100644 Source/Utils/ERF_solve_with_EB_mlmg.cpp diff --git a/Exec/DevTests/EB_Test/ERF_prob.cpp b/Exec/DevTests/EB_Test/ERF_prob.cpp index e6d0a0a81..d67e3a7bb 100644 --- a/Exec/DevTests/EB_Test/ERF_prob.cpp +++ b/Exec/DevTests/EB_Test/ERF_prob.cpp @@ -67,8 +67,8 @@ Problem::init_custom_pert( Array4 const& z_vel_pert, Array4 const& /*r_hse*/, Array4 const& /*p_hse*/, - Array4 const& z_nd, - Array4 const& z_cc, + Array4 const& /*z_nd*/, + Array4 const& /*z_cc*/, GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -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 @@ -157,6 +157,9 @@ Problem::init_custom_pert( amrex::Array4 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 @@ -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; }); } diff --git a/Exec/DevTests/EB_Test/GNUmakefile b/Exec/DevTests/EB_Test/GNUmakefile index 2ba1293c6..0ad0b3ebf 100644 --- a/Exec/DevTests/EB_Test/GNUmakefile +++ b/Exec/DevTests/EB_Test/GNUmakefile @@ -20,7 +20,6 @@ USE_SYCL = FALSE # Debugging DEBUG = FALSE -DEBUG = TRUE TEST = TRUE USE_ASSERTION = TRUE diff --git a/Exec/DevTests/EB_Test/inputs b/Exec/DevTests/EB_Test/inputs index 51ab1b2dd..bdc31e471 100644 --- a/Exec/DevTests/EB_Test/inputs +++ b/Exec/DevTests/EB_Test/inputs @@ -5,6 +5,10 @@ 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 @@ -12,22 +16,22 @@ 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 diff --git a/Source/EB/ERF_initEB.cpp b/Source/EB/ERF_initEB.cpp index 4469f77be..3be6c423b 100644 --- a/Source/EB/ERF_initEB.cpp +++ b/Source/EB/ERF_initEB.cpp @@ -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; diff --git a/Source/ERF.H b/Source/ERF.H index a0330870b..aa8c7a8fa 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -176,6 +176,10 @@ public: void solve_with_fft (int lev, amrex::MultiFab& rhs, amrex::MultiFab& p, amrex::Array& fluxes); #endif +#ifdef ERF_USE_EB + void solve_with_EB_mlmg (int lev, amrex::Vector& rhs, amrex::Vector& p, + amrex::Vector>& fluxes); +#endif void solve_with_gmres (int lev, amrex::Vector& rhs, amrex::Vector& p, amrex::Vector>& fluxes); @@ -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 > > m_factory; +#ifdef ERF_USE_EB + + amrex::Vector > m_factory; [[nodiscard]] amrex::FabFactory 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(*m_factory[lev]); @@ -1361,6 +1366,8 @@ private: [[nodiscard]] static int nghost_eb_full () { return 4; } +#else + amrex::Vector > > m_factory; #endif #ifdef ERF_USE_FFT diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 581b45e9e..2d056fbbd 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -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 diff --git a/Source/TimeIntegration/ERF_TI_no_substep_fun.H b/Source/TimeIntegration/ERF_TI_no_substep_fun.H index d40368396..1e9419b6d 100644 --- a/Source/TimeIntegration/ERF_TI_no_substep_fun.H +++ b/Source/TimeIntegration/ERF_TI_no_substep_fun.H @@ -106,24 +106,47 @@ }); } else { // Fixed or no terrain +#ifdef ERF_USE_EB + const Array4& 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 diff --git a/Source/Utils/ERF_PoissonSolve.cpp b/Source/Utils/ERF_PoissonSolve.cpp index e4bb52e37..2a98ee4a7 100644 --- a/Source/Utils/ERF_PoissonSolve.cpp +++ b/Source/Utils/ERF_PoissonSolve.cpp @@ -14,8 +14,10 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& 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 ba_tmp; ba_tmp.push_back(mom_mf[Vars::cons].boxArray()); @@ -28,12 +30,21 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& mom_mf, Mult Vector phi; Vector > 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(); @@ -49,6 +60,10 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& 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) @@ -96,6 +111,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& mom_mf, Mult computeDivergence(rhs[0], rho0_u_const, geom_tmp[0]); } +#endif Real rhsnorm = rhs[0].norm0(); @@ -115,6 +131,9 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& 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) { @@ -130,6 +149,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& mom_mf, Mult solve_with_mlmg(lev, rhs, phi, fluxes); } +#endif // **************************************************************************** // Subtract dt grad(phi) from the momenta (rho0u, rho0v, Omega) diff --git a/Source/Utils/ERF_solve_with_EB_mlmg.cpp b/Source/Utils/ERF_solve_with_EB_mlmg.cpp new file mode 100644 index 000000000..6e66eaea1 --- /dev/null +++ b/Source/Utils/ERF_solve_with_EB_mlmg.cpp @@ -0,0 +1,170 @@ +#include "ERF.H" +#include "ERF_Utils.H" + +#include +#include + +using namespace amrex; + +/** + * Define the domain boundary conditions for the (optional) Poisson solve + */ + +using BCType = LinOpBCType; + +/** + * Solve the Poisson equation using EB_enabled MLMG + * Note that the level may or may not be level 0. + */ +void ERF::solve_with_EB_mlmg (int lev, Vector& rhs, Vector& phi, Vector>& fluxes) +{ + BL_PROFILE("ERF::solve_with_EB_mlmg()"); + + auto const dom_lo = lbound(geom[lev].Domain()); + auto const dom_hi = ubound(geom[lev].Domain()); + + LPInfo info; + // Allow a hidden direction if the domain is one cell wide in any lateral direction + if (dom_lo.x == dom_hi.x) { + info.setHiddenDirection(0); + } else if (dom_lo.y == dom_hi.y) { + info.setHiddenDirection(1); + } + + // Make sure the solver only sees the levels over which we are solving + Vector ba_tmp; ba_tmp.push_back(rhs[0].boxArray()); + Vector dm_tmp; dm_tmp.push_back(rhs[0].DistributionMap()); + Vector geom_tmp; geom_tmp.push_back(geom[lev]); + + auto bclo = get_projection_bc(Orientation::low); + auto bchi = get_projection_bc(Orientation::high); + + // amrex::Print() << "BCLO " << bclo[0] << " " << bclo[1] << " " << bclo[2] << std::endl; + // amrex::Print() << "BCHI " << bchi[0] << " " << bchi[1] << " " << bchi[2] << std::endl; + + Real reltol = solverChoice.poisson_reltol; + Real abstol = solverChoice.poisson_abstol; + + // **************************************************************************** + // Multigrid solve + // **************************************************************************** + + MLEBABecLap mleb (geom_tmp, ba_tmp, dm_tmp, info, {m_factory[lev].get()}); + + mleb.setMaxOrder(2); + mleb.setDomainBC(bclo, bchi); + mleb.setLevelBC(0, nullptr); + + // + // This sets A = 0, B = 1 so that + // the operator A alpha - b del dot beta grad to b + // becomes - del dot beta grad + // + mleb.setScalars(0.0, 1.0); + + Array bcoef; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + bcoef[idim].define(convert(ba_tmp[0],IntVect::TheDimensionVector(idim)), + dm_tmp[0], 1, 0, MFInfo(), *m_factory[lev]); + bcoef[idim].setVal(1.0); + } + mleb.setBCoeffs(0, amrex::GetArrOfConstPtrs(bcoef)); + + MLMG mlmg(mleb); + + int max_iter = 100; + mlmg.setMaxIter(max_iter); + mlmg.setVerbose(mg_verbose); + mlmg.setBottomVerbose(0); + + mlmg.solve(GetVecOfPtrs(phi), GetVecOfConstPtrs(rhs), reltol, abstol); + + mlmg.getFluxes(GetVecOfArrOfPtrs(fluxes)); + + phi[0].FillBoundary(geom[lev].periodicity()); + + // **************************************************************************** + // Impose bc's on pprime + // **************************************************************************** +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(phi[0],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Array4 const& pp_arr = phi[0].array(mfi); + Box const& bx = mfi.tilebox(); + auto const bx_lo = lbound(bx); + auto const bx_hi = ubound(bx); + if (bx_lo.x == dom_lo.x) { + auto bc_type = domain_bc_type[Orientation(0,Orientation::low)]; + if (bc_type == "Outflow" || bc_type == "Open") { + ParallelFor(makeSlab(bx,0,dom_lo.x), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pp_arr(i-1,j,k) = -pp_arr(i,j,k); + }); + } else { + ParallelFor(makeSlab(bx,0,dom_lo.x), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pp_arr(i-1,j,k) = pp_arr(i,j,k); + }); + } + } + if (bx_lo.y == dom_lo.y) { + auto bc_type = domain_bc_type[Orientation(1,Orientation::low)]; + if (bc_type == "Outflow" || bc_type == "Open") { + ParallelFor(makeSlab(bx,1,dom_lo.y), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pp_arr(i,j-1,k) = -pp_arr(i,j,k); + }); + } else { + ParallelFor(makeSlab(bx,1,dom_lo.y), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pp_arr(i,j-1,k) = pp_arr(i,j,k); + }); + } + } + if (bx_lo.z == dom_lo.z) { + ParallelFor(makeSlab(bx,2,dom_lo.z), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pp_arr(i,j,k-1) = pp_arr(i,j,k); + }); + } + if (bx_hi.x == dom_hi.x) { + auto bc_type = domain_bc_type[Orientation(0,Orientation::high)]; + if (bc_type == "Outflow" || bc_type == "Open") { + ParallelFor(makeSlab(bx,0,dom_hi.x), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pp_arr(i+1,j,k) = -pp_arr(i,j,k); + }); + } else { + ParallelFor(makeSlab(bx,0,dom_hi.x), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pp_arr(i+1,j,k) = pp_arr(i,j,k); + }); + } + } + if (bx_hi.y == dom_hi.y) { + auto bc_type = domain_bc_type[Orientation(1,Orientation::high)]; + if (bc_type == "Outflow" || bc_type == "Open") { + ParallelFor(makeSlab(bx,1,dom_hi.y), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pp_arr(i,j+1,k) = -pp_arr(i,j,k); + }); + } else { + ParallelFor(makeSlab(bx,1,dom_hi.y), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pp_arr(i,j+1,k) = pp_arr(i,j,k); + }); + } + } + if (bx_hi.z == dom_hi.z) { + ParallelFor(makeSlab(bx,2,dom_hi.z), [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pp_arr(i,j,k+1) = pp_arr(i,j,k); + }); + } + } // mfi + + // Now overwrite with periodic fill outside domain and fine-fine fill inside + phi[0].FillBoundary(geom[lev].periodicity()); +} diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package index 3fc5cd0ec..8667c9a3d 100644 --- a/Source/Utils/Make.package +++ b/Source/Utils/Make.package @@ -10,6 +10,9 @@ CEXE_headers += ERF_TerrainMetrics.H CEXE_headers += ERF_TileNoZ.H CEXE_headers += ERF_Utils.H +CEXE_headers += ERF_TerrainPoisson_3D_K.H +CEXE_headers += ERF_MLTerrainPoisson.H + CEXE_headers += ERF_ParFunctions.H CEXE_headers += ERF_Sat_methods.H @@ -29,3 +32,7 @@ CEXE_sources += ERF_PoissonSolve_tb.cpp CEXE_sources += ERF_solve_with_fft.cpp CEXE_sources += ERF_solve_with_gmres.cpp CEXE_sources += ERF_solve_with_mlmg.cpp + +ifeq ($(USE_EB), TRUE) +CEXE_sources += ERF_solve_with_EB_mlmg.cpp +endif