Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Mesh Refinement with Explicit Solver #835

Merged
merged 7 commits into from
Dec 7, 2022
Merged

Mesh Refinement with Explicit Solver #835

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
65d610a
Explicit mesh refinement
AlexanderSinn Dec 2, 2022
eba27ae
account for hpmg boundary position
AlexanderSinn Dec 2, 2022
656dc87
fix doc
AlexanderSinn Dec 2, 2022
823cc8f
make amrex mg solver work
AlexanderSinn Dec 5, 2022
d183f71
add suggestions
AlexanderSinn Dec 7, 2022
827d9e0
fix missing greek question mark
AlexanderSinn Dec 7, 2022
2cead96
fix doc
AlexanderSinn Dec 7, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 5 additions & 4 deletions src/Hipace.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -209,8 +209,9 @@ public:
*
* \param[lev] MR level
* \param[in] which_slice defines if this or the salame slice is handled
* \param[in] islice longitudinal slice
*/
void ExplicitMGSolveBxBy (const int lev, const int which_slice);
void ExplicitMGSolveBxBy (const int lev, const int which_slice, const int islice);

/** \brief Reset plasma and field slice quantities to initial value.
*
Expand Down Expand Up @@ -403,12 +404,12 @@ private:

#ifdef AMREX_USE_LINEAR_SOLVERS
/** Linear operator for the explicit Bx and By solver */
std::unique_ptr<amrex::MLALaplacian> m_mlalaplacian;
amrex::Vector<std::unique_ptr<amrex::MLALaplacian>> m_mlalaplacian;
/** Geometric multigrid solver class, for the explicit Bx and By solver */
std::unique_ptr<amrex::MLMG> m_mlmg;
amrex::Vector<std::unique_ptr<amrex::MLMG>> m_mlmg;
#endif
/** hpmg solver for the explicit Bx and by solver */
std::unique_ptr<hpmg::MultiGrid> m_hpmg;
amrex::Vector<std::unique_ptr<hpmg::MultiGrid>> m_hpmg;
public:
/** Used to sort the beam particles into boxes for pipelining */
amrex::Vector<BoxSorter> m_box_sorters;
Expand Down
100 changes: 55 additions & 45 deletions src/Hipace.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -160,6 +160,8 @@ Hipace::Hipace () :
if (solver == "explicit") m_explicit = true;
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_explicit || !m_multi_beam.AnySpeciesSalame(),
"Cannot use SALAME algorithm with predictor-corrector solver");
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(maxLevel()==0 || !m_multi_beam.AnySpeciesSalame(),
"Cannot use SALAME algorithm with mesh refinement");
queryWithParser(pph, "MG_tolerance_rel", m_MG_tolerance_rel);
queryWithParser(pph, "MG_tolerance_abs", m_MG_tolerance_abs);
queryWithParser(pph, "MG_verbose", m_MG_verbose);
Expand All @@ -171,8 +173,6 @@ Hipace::Hipace () :

if (maxLevel() > 0) {
AMREX_ALWAYS_ASSERT(maxLevel() < 2);
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!m_explicit, "Mesh refinement + explicit solver is not yet"
" supported! Please use hipace.bxby_solver = predictor-corrector");
amrex::Array<amrex::Real, AMREX_SPACEDIM> loc_array;
getWithParser(pph, "patch_lo", loc_array);
for (int idim=0; idim<AMREX_SPACEDIM; ++idim) patch_lo[idim] = loc_array[idim];
Expand Down Expand Up @@ -587,18 +587,25 @@ Hipace::SolveOneSlice (int islice_coarse, const int ibox, int step,

if (m_multi_plasma.IonizationOn() && m_do_tiling) m_multi_plasma.TileSort(bx, geom[lev]);

} // end for (int isubslice = nsubslice-1; isubslice >= 0; --isubslice)
if (lev == 0) {
// Advance laser slice by 1 step and store result to 3D array
m_multi_laser.AdvanceSlice(m_fields, Geom(0), m_dt, step);
m_multi_laser.Copy(islice_coarse, true);
}

// Advance laser slice by 1 step and store result to 3D array
m_multi_laser.AdvanceSlice(m_fields, Geom(0), m_dt, step);
m_multi_laser.Copy(islice_coarse, true);
if (lev != 0) {
// shift slices of level 1
m_fields.ShiftSlices(lev, islice_coarse, Geom(0), patch_lo[2], patch_hi[2]);
}
} // end for (int isubslice = nsubslice-1; isubslice >= 0; --isubslice)

// After this, the parallel context is the full 3D communicator again
amrex::ParallelContext::pop();

} // end for (int lev = 0; lev <= finestLevel(); ++lev)

// shift slices of all levels
m_fields.ShiftSlices(finestLevel()+1, islice_coarse, Geom(0), patch_lo[2], patch_hi[2]);
// shift level 0
m_fields.ShiftSlices(0, islice_coarse, Geom(0), patch_lo[2], patch_hi[2]);
}


Expand Down Expand Up @@ -650,7 +657,7 @@ Hipace::ExplicitSolveOneSubSlice (const int lev, const int step, const int ibox,
m_multi_plasma.ExplicitDeposition(m_fields, m_multi_laser, geom[lev], lev);

// Solves Bx, By using Sx, Sy and chi
ExplicitMGSolveBxBy(lev, WhichSlice::This);
ExplicitMGSolveBxBy(lev, WhichSlice::This, islice);

const bool do_salame = m_multi_beam.isSalameNow(step, islice_local, beam_bin);
if (do_salame) {
Expand Down Expand Up @@ -820,7 +827,7 @@ Hipace::InitializeSxSyWithBeam (const int lev)


void
Hipace::ExplicitMGSolveBxBy (const int lev, const int which_slice)
Hipace::ExplicitMGSolveBxBy (const int lev, const int which_slice, const int islice)
{
HIPACE_PROFILE("Hipace::ExplicitMGSolveBxBy()");
amrex::ParallelContext::push(m_comm_xy);
Expand All @@ -839,86 +846,89 @@ Hipace::ExplicitMGSolveBxBy (const int lev, const int which_slice)
AMREX_ALWAYS_ASSERT(Comps[which_slice]["Bx"] + 1 == Comps[which_slice]["By"]);
AMREX_ALWAYS_ASSERT(Comps[which_slice]["Sy"] + 1 == Comps[which_slice]["Sx"]);

// construct slice geometry
// Set the lo and hi of domain and probdomain in the z direction
amrex::RealBox tmp_probdom({AMREX_D_DECL(Geom(lev).ProbLo(Direction::x),
Geom(lev).ProbLo(Direction::y),
Geom(lev).ProbLo(Direction::z))},
{AMREX_D_DECL(Geom(lev).ProbHi(Direction::x),
Geom(lev).ProbHi(Direction::y),
Geom(lev).ProbHi(Direction::z))});
amrex::Box tmp_dom = Geom(lev).Domain();
const amrex::Real hi = Geom(lev).ProbHi(Direction::z);
const amrex::Real lo = hi - Geom(lev).CellSize(Direction::z);
tmp_probdom.setLo(Direction::z, lo);
tmp_probdom.setHi(Direction::z, hi);
tmp_dom.setSmall(Direction::z, 0);
tmp_dom.setBig(Direction::z, 0);
amrex::Geometry slice_geom = amrex::Geometry(
tmp_dom, tmp_probdom, Geom(lev).Coord(), Geom(lev).isPeriodic());

slice_geom.setPeriodicity({0,0,0});

amrex::MultiFab& slicemf_BxBySySx = m_fields.getSlices(lev, which_slice);
amrex::MultiFab BxBy (slicemf_BxBySySx, amrex::make_alias, Comps[which_slice]["Bx"], 2);
amrex::MultiFab SySx (slicemf_BxBySySx, amrex::make_alias, Comps[which_slice]["Sy"], 2);

amrex::MultiFab& slicemf_chi = m_fields.getSlices(lev, which_slice_chi);
amrex::MultiFab Mult (slicemf_chi, amrex::make_alias, Comps[which_slice_chi]["chi"], ncomp_chi);

if (lev!=0) {
m_fields.SetBoundaryCondition(Geom(), lev, "Bx", islice,
m_fields.getField(lev, which_slice, "Sy"));
m_fields.SetBoundaryCondition(Geom(), lev, "By", islice,
m_fields.getField(lev, which_slice, "Sx"));
}

#ifdef AMREX_USE_LINEAR_SOLVERS
if (m_use_amrex_mlmg) {
// Copy chi to chi2
m_fields.duplicate(lev, which_slice_chi, {"chi2"}, which_slice_chi, {"chi"});
amrex::Gpu::streamSynchronize();
if (!m_mlalaplacian){
if (m_mlalaplacian.size()==0) {
MaxThevenet marked this conversation as resolved.
Show resolved Hide resolved
m_mlalaplacian.resize(maxLevel()+1);
m_mlmg.resize(maxLevel()+1);
}

// construct slice geometry
const amrex::RealBox slice_box{slicemf_BxBySySx.boxArray()[0], m_slice_geom[lev].CellSize(),
m_slice_geom[lev].ProbLo()};
amrex::Geometry slice_geom{slicemf_BxBySySx.boxArray()[0], slice_box,
m_slice_geom[lev].CoordInt(), {0,0,0}};

if (!m_mlalaplacian[lev]){
// If first call, initialize the MG solver
amrex::LPInfo lpinfo{};
lpinfo.setHiddenDirection(2).setAgglomeration(false).setConsolidation(false);

// make_unique requires explicit types
m_mlalaplacian = std::make_unique<amrex::MLALaplacian>(
m_mlalaplacian[lev] = std::make_unique<amrex::MLALaplacian>(
amrex::Vector<amrex::Geometry>{slice_geom},
amrex::Vector<amrex::BoxArray>{S.boxArray()},
amrex::Vector<amrex::DistributionMapping>{S.DistributionMap()},
amrex::Vector<amrex::BoxArray>{slicemf_BxBySySx.boxArray()},
amrex::Vector<amrex::DistributionMapping>{slicemf_BxBySySx.DistributionMap()},
lpinfo,
amrex::Vector<amrex::FabFactory<amrex::FArrayBox> const*>{}, 2);

m_mlalaplacian->setDomainBC(
m_mlalaplacian[lev]->setDomainBC(
{AMREX_D_DECL(amrex::LinOpBCType::Dirichlet,
amrex::LinOpBCType::Dirichlet,
amrex::LinOpBCType::Dirichlet)},
{AMREX_D_DECL(amrex::LinOpBCType::Dirichlet,
amrex::LinOpBCType::Dirichlet,
amrex::LinOpBCType::Dirichlet)});

m_mlmg = std::make_unique<amrex::MLMG>(*m_mlalaplacian);
m_mlmg->setVerbose(m_MG_verbose);
m_mlmg[lev] = std::make_unique<amrex::MLMG>(*(m_mlalaplacian[lev]));
m_mlmg[lev]->setVerbose(m_MG_verbose);
}

// BxBy is assumed to have at least one ghost cell in x and y.
// The ghost cells outside the domain should contain Dirichlet BC values.
BxBy.setDomainBndry(0.0, slice_geom); // Set Dirichlet BC to zero
m_mlalaplacian->setLevelBC(0, &BxBy);
m_mlalaplacian[lev]->setLevelBC(0, &BxBy);

m_mlalaplacian->setACoeffs(0, Mult);
m_mlalaplacian[lev]->setACoeffs(0, Mult);

// amrex solves ascalar A phi - bscalar Laplacian(phi) = rhs
// So we solve Delta BxBy - A * BxBy = S
m_mlalaplacian->setScalars(-1.0, -1.0);
m_mlalaplacian[lev]->setScalars(-1.0, -1.0);

m_mlmg->solve({&BxBy}, {&SySx}, m_MG_tolerance_rel, m_MG_tolerance_abs);
m_mlmg[lev]->solve({&BxBy}, {&SySx}, m_MG_tolerance_rel, m_MG_tolerance_abs);
} else
#endif
{
AMREX_ALWAYS_ASSERT(slicemf_BxBySySx.boxArray().size() == 1);
AMREX_ALWAYS_ASSERT(slicemf_chi.boxArray().size() == 1);
if (!m_hpmg) {
m_hpmg = std::make_unique<hpmg::MultiGrid>(slice_geom);
if (m_hpmg.size()==0) {
MaxThevenet marked this conversation as resolved.
Show resolved Hide resolved
m_hpmg.resize(maxLevel()+1);
}
if (!m_hpmg[lev]) {
m_hpmg[lev] = std::make_unique<hpmg::MultiGrid>(m_slice_geom[lev].CellSize(0),
m_slice_geom[lev].CellSize(1),
slicemf_BxBySySx.boxArray()[0]);
}
const int max_iters = 200;
m_hpmg->solve1(BxBy[0], SySx[0], Mult[0], m_MG_tolerance_rel, m_MG_tolerance_abs,
max_iters, m_MG_verbose);
m_hpmg[lev]->solve1(BxBy[0], SySx[0], Mult[0], m_MG_tolerance_rel, m_MG_tolerance_abs,
max_iters, m_MG_verbose);
}
amrex::ParallelContext::pop();
}
Expand Down
8 changes: 5 additions & 3 deletions src/fields/Fields.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -171,13 +171,13 @@ public:
* to compute each slice. The current slice is always stored in index 1.
* Hence, after one slice is computed, slices must be shifted by 1 element.
*
* \param[in] nlev number of MR level
* \param[in] lev the MR level
* \param[in] islice current slice to be shifted
* \param[in] geom geometry
* \param[in] patch_lo lower longitudinal end of refined level
* \param[in] patch_hi upper longitudinal end of refined level
*/
void ShiftSlices (int nlev, int islice, const amrex::Geometry geom, const amrex::Real patch_lo,
void ShiftSlices (int lev, int islice, const amrex::Geometry geom, const amrex::Real patch_lo,
const amrex::Real patch_hi);

/** add rho of the ions to rho (this slice) */
Expand All @@ -191,9 +191,11 @@ public:
* \param[in] lev current level
* \param[in] component which can be Psi, Ez, By, Bx ...
* \param[in] islice longitudinal slice
* \param[in,out] staging_area Target MultiFab where the boundary condition is applied
*/
void SetBoundaryCondition (amrex::Vector<amrex::Geometry> const& geom, const int lev,
std::string component, const int islice);
std::string component, const int islice,
amrex::MultiFab&& staging_area);

/** \brief Interpolate values from coarse grid to the fine grid
*
Expand Down
Loading