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: multiple Poisson solvers #527

Merged
merged 11 commits into from
Jun 14, 2021
Merged

Mesh refinement: multiple Poisson solvers #527

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
7c51dad
debug printouts
SeverinDiederichs Jun 10, 2021
3f003db
add IO for refined level
SeverinDiederichs Jun 11, 2021
d9d7189
cleaning
SeverinDiederichs Jun 11, 2021
0d5d958
clearer explanation for workaround
SeverinDiederichs Jun 11, 2021
6b7619a
add const for fabs in IO, which was removed for testing purposes
SeverinDiederichs Jun 11, 2021
13f2edb
more cleaning
SeverinDiederichs Jun 11, 2021
c2dcab6
fix typo to enable untouched IO without MR
SeverinDiederichs Jun 11, 2021
e1c6828
call reset based on length of m_openPMD_series array
SeverinDiederichs Jun 11, 2021
17b0d8d
using multiple Poisson solvers
SeverinDiederichs Jun 12, 2021
fe94d7d
solve merge conflicts
SeverinDiederichs Jun 14, 2021
822c06c
forgot to remove unused parameter in merge conflict
SeverinDiederichs Jun 14, 2021
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
2 changes: 1 addition & 1 deletion src/Hipace.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -279,7 +279,7 @@ Hipace::MakeNewLevelFromScratch (
DefineSliceGDB(ba, dm);
// Note: we pass ba[0] as a dummy box, it will be resized properly in the loop over boxes in Evolve
m_diags.AllocData(lev, ba[0], Comps[WhichSlice::This]["N"], Geom(lev));
m_fields.AllocData(lev, Geom(lev), m_slice_ba, m_slice_dm);
m_fields.AllocData(lev, Geom(), m_slice_ba, m_slice_dm);
}

void
Expand Down
4 changes: 2 additions & 2 deletions src/fields/Fields.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -79,11 +79,11 @@ public:
* \param[in] slice_dm DistributionMapping for the slice
*/
void AllocData (
int lev, amrex::Geometry const& geom, const amrex::BoxArray& slice_ba,
int lev, amrex::Vector<amrex::Geometry> const& geom, const amrex::BoxArray& slice_ba,
const amrex::DistributionMapping& slice_dm);

/** Class to handle transverse FFT Poisson solver on 1 slice */
std::unique_ptr<FFTPoissonSolver> m_poisson_solver;
amrex::Vector<std::unique_ptr<FFTPoissonSolver>> m_poisson_solver;
/** get function for the 2D slices */
amrex::Vector<std::array<amrex::MultiFab, m_nslices>>& getSlices () {return m_slices; }
/** get function for the 2D slices
Expand Down
49 changes: 22 additions & 27 deletions src/fields/Fields.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,39 +14,34 @@ Fields::Fields (Hipace const* a_hipace)

void
Fields::AllocData (
int lev, amrex::Geometry const& geom, const amrex::BoxArray& slice_ba,
int lev, amrex::Vector<amrex::Geometry> const& geom, const amrex::BoxArray& slice_ba,
const amrex::DistributionMapping& slice_dm)
{
HIPACE_PROFILE("Fields::AllocData()");
// Need at least 1 guard cell transversally for transverse derivative
int nguards_xy = std::max(1, Hipace::m_depos_order_xy);
m_slices_nguards = {nguards_xy, nguards_xy, 0};

// Only xy slices need guard cells, there is no deposition to/gather from the output array F.
amrex::IntVect nguards_F = amrex::IntVect(0,0,0);

for (int islice=0; islice<WhichSlice::N; islice++) {
m_slices[lev][islice].define(
slice_ba, slice_dm, Comps[islice]["N"], m_slices_nguards,
amrex::MFInfo().SetArena(amrex::The_Arena()));
m_slices[lev][islice].setVal(0.0);
}

// FIXME the Poisson solver must be constructed per level, here it is only constructed for lev 0
if (lev>0) return;
// The Poisson solver operates on transverse slices only.
// The constructor takes the BoxArray and the DistributionMap of a slice,
// so the FFTPlans are built on a slice.
if (m_do_dirichlet_poisson){
m_poisson_solver = std::unique_ptr<FFTPoissonSolverDirichlet>(
m_poisson_solver.push_back(std::unique_ptr<FFTPoissonSolverDirichlet>(
new FFTPoissonSolverDirichlet(getSlices(lev, WhichSlice::This).boxArray(),
getSlices(lev, WhichSlice::This).DistributionMap(),
geom));
geom[lev])) );
} else {
m_poisson_solver = std::unique_ptr<FFTPoissonSolverPeriodic>(
m_poisson_solver.push_back(std::unique_ptr<FFTPoissonSolverPeriodic>(
new FFTPoissonSolverPeriodic(getSlices(lev, WhichSlice::This).boxArray(),
getSlices(lev, WhichSlice::This).DistributionMap(),
geom));
geom[lev])) );
}
}

Expand Down Expand Up @@ -255,14 +250,14 @@ Fields::SolvePoissonExmByAndEypBx (amrex::Geometry const& geom, const MPI_Comm&
Comps[WhichSlice::This]["Psi"], 1);

// calculating the right-hand side 1/episilon0 * -(rho-Jz/c)
amrex::MultiFab::Copy(m_poisson_solver->StagingArea(), getSlices(lev, WhichSlice::This),
amrex::MultiFab::Copy(m_poisson_solver[lev]->StagingArea(), getSlices(lev, WhichSlice::This),
Comps[WhichSlice::This]["jz"], 0, 1, 0);
m_poisson_solver->StagingArea().mult(-1./phys_const.c);
amrex::MultiFab::Add(m_poisson_solver->StagingArea(), getSlices(lev, WhichSlice::This),
m_poisson_solver[lev]->StagingArea().mult(-1./phys_const.c);
amrex::MultiFab::Add(m_poisson_solver[lev]->StagingArea(), getSlices(lev, WhichSlice::This),
Comps[WhichSlice::This]["rho"], 0, 1, 0);
m_poisson_solver->StagingArea().mult(-1./phys_const.ep0);
m_poisson_solver[lev]->StagingArea().mult(-1./phys_const.ep0);

m_poisson_solver->SolvePoissonEquation(lhs);
m_poisson_solver[lev]->SolvePoissonEquation(lhs);

/* ---------- Transverse FillBoundary Psi ---------- */
amrex::ParallelContext::push(m_comm_xy);
Expand Down Expand Up @@ -306,7 +301,7 @@ Fields::SolvePoissonEz (amrex::Geometry const& geom, const int lev)
// from the slice MF, and store in the staging area of poisson_solver
TransverseDerivative(
getSlices(lev, WhichSlice::This),
m_poisson_solver->StagingArea(),
m_poisson_solver[lev]->StagingArea(),
Direction::x,
geom.CellSize(Direction::x),
1./(phys_const.ep0*phys_const.c),
Expand All @@ -315,7 +310,7 @@ Fields::SolvePoissonEz (amrex::Geometry const& geom, const int lev)

TransverseDerivative(
getSlices(lev, WhichSlice::This),
m_poisson_solver->StagingArea(),
m_poisson_solver[lev]->StagingArea(),
Direction::y,
geom.CellSize(Direction::y),
1./(phys_const.ep0*phys_const.c),
Expand All @@ -324,7 +319,7 @@ Fields::SolvePoissonEz (amrex::Geometry const& geom, const int lev)
// Solve Poisson equation.
// The RHS is in the staging area of poisson_solver.
// The LHS will be returned as lhs.
m_poisson_solver->SolvePoissonEquation(lhs);
m_poisson_solver[lev]->SolvePoissonEquation(lhs);
}

void
Expand All @@ -338,7 +333,7 @@ Fields::SolvePoissonBx (amrex::MultiFab& Bx_iter, amrex::Geometry const& geom, c
// and store in the staging area of poisson_solver
TransverseDerivative(
getSlices(lev, WhichSlice::This),
m_poisson_solver->StagingArea(),
m_poisson_solver[lev]->StagingArea(),
Direction::y,
geom.CellSize(Direction::y),
-phys_const.mu0,
Expand All @@ -348,7 +343,7 @@ Fields::SolvePoissonBx (amrex::MultiFab& Bx_iter, amrex::Geometry const& geom, c
LongitudinalDerivative(
getSlices(lev, WhichSlice::Previous1),
getSlices(lev, WhichSlice::Next),
m_poisson_solver->StagingArea(),
m_poisson_solver[lev]->StagingArea(),
geom.CellSize(Direction::z),
phys_const.mu0,
SliceOperatorType::Add,
Expand All @@ -357,7 +352,7 @@ Fields::SolvePoissonBx (amrex::MultiFab& Bx_iter, amrex::Geometry const& geom, c
// Solve Poisson equation.
// The RHS is in the staging area of poisson_solver.
// The LHS will be returned as lhs.
m_poisson_solver->SolvePoissonEquation(Bx_iter);
m_poisson_solver[lev]->SolvePoissonEquation(Bx_iter);
}

void
Expand All @@ -371,7 +366,7 @@ Fields::SolvePoissonBy (amrex::MultiFab& By_iter, amrex::Geometry const& geom, c
// and store in the staging area of poisson_solver
TransverseDerivative(
getSlices(lev, WhichSlice::This),
m_poisson_solver->StagingArea(),
m_poisson_solver[lev]->StagingArea(),
Direction::x,
geom.CellSize(Direction::x),
phys_const.mu0,
Expand All @@ -381,7 +376,7 @@ Fields::SolvePoissonBy (amrex::MultiFab& By_iter, amrex::Geometry const& geom, c
LongitudinalDerivative(
getSlices(lev, WhichSlice::Previous1),
getSlices(lev, WhichSlice::Next),
m_poisson_solver->StagingArea(),
m_poisson_solver[lev]->StagingArea(),
geom.CellSize(Direction::z),
-phys_const.mu0,
SliceOperatorType::Add,
Expand All @@ -390,7 +385,7 @@ Fields::SolvePoissonBy (amrex::MultiFab& By_iter, amrex::Geometry const& geom, c
// Solve Poisson equation.
// The RHS is in the staging area of poisson_solver.
// The LHS will be returned as lhs.
m_poisson_solver->SolvePoissonEquation(By_iter);
m_poisson_solver[lev]->SolvePoissonEquation(By_iter);
}

void
Expand All @@ -407,7 +402,7 @@ Fields::SolvePoissonBz (amrex::Geometry const& geom, const int lev)
// from the slice MF, and store in the staging area of m_poisson_solver
TransverseDerivative(
getSlices(lev, WhichSlice::This),
m_poisson_solver->StagingArea(),
m_poisson_solver[lev]->StagingArea(),
Direction::y,
geom.CellSize(Direction::y),
phys_const.mu0,
Expand All @@ -416,7 +411,7 @@ Fields::SolvePoissonBz (amrex::Geometry const& geom, const int lev)

TransverseDerivative(
getSlices(lev, WhichSlice::This),
m_poisson_solver->StagingArea(),
m_poisson_solver[lev]->StagingArea(),
Direction::x,
geom.CellSize(Direction::x),
-phys_const.mu0,
Expand All @@ -425,7 +420,7 @@ Fields::SolvePoissonBz (amrex::Geometry const& geom, const int lev)
// Solve Poisson equation.
// The RHS is in the staging area of m_poisson_solver.
// The LHS will be returned as lhs.
m_poisson_solver->SolvePoissonEquation(lhs);
m_poisson_solver[lev]->SolvePoissonEquation(lhs);
}

void
Expand Down