Skip to content

Commit cafc744

Browse files
SeverinDiederichsSeverinDiederichs
authored
Mesh refinement: multiple Poisson solvers (#527)
* debug printouts * add IO for refined level * cleaning * clearer explanation for workaround * add const for fabs in IO, which was removed for testing purposes * more cleaning * fix typo to enable untouched IO without MR * call reset based on length of m_openPMD_series array * using multiple Poisson solvers * forgot to remove unused parameter in merge conflict
1 parent 5931ef7 commit cafc744

3 files changed

Lines changed: 25 additions & 30 deletions

File tree

‎src/Hipace.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -279,7 +279,7 @@ Hipace::MakeNewLevelFromScratch (
279279
DefineSliceGDB(ba, dm);
280280
// Note: we pass ba[0] as a dummy box, it will be resized properly in the loop over boxes in Evolve
281281
m_diags.AllocData(lev, ba[0], Comps[WhichSlice::This]["N"], Geom(lev));
282-
m_fields.AllocData(lev, Geom(lev), m_slice_ba, m_slice_dm);
282+
m_fields.AllocData(lev, Geom(), m_slice_ba, m_slice_dm);
283283
}
284284

285285
void

‎src/fields/Fields.H‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -79,11 +79,11 @@ public:
7979
* \param[in] slice_dm DistributionMapping for the slice
8080
*/
8181
void AllocData (
82-
int lev, amrex::Geometry const& geom, const amrex::BoxArray& slice_ba,
82+
int lev, amrex::Vector<amrex::Geometry> const& geom, const amrex::BoxArray& slice_ba,
8383
const amrex::DistributionMapping& slice_dm);
8484

8585
/** Class to handle transverse FFT Poisson solver on 1 slice */
86-
std::unique_ptr<FFTPoissonSolver> m_poisson_solver;
86+
amrex::Vector<std::unique_ptr<FFTPoissonSolver>> m_poisson_solver;
8787
/** get function for the 2D slices */
8888
amrex::Vector<std::array<amrex::MultiFab, m_nslices>>& getSlices () {return m_slices; }
8989
/** get function for the 2D slices

‎src/fields/Fields.cpp‎

Lines changed: 22 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -14,39 +14,34 @@ Fields::Fields (Hipace const* a_hipace)
1414

1515
void
1616
Fields::AllocData (
17-
int lev, amrex::Geometry const& geom, const amrex::BoxArray& slice_ba,
17+
int lev, amrex::Vector<amrex::Geometry> const& geom, const amrex::BoxArray& slice_ba,
1818
const amrex::DistributionMapping& slice_dm)
1919
{
2020
HIPACE_PROFILE("Fields::AllocData()");
2121
// Need at least 1 guard cell transversally for transverse derivative
2222
int nguards_xy = std::max(1, Hipace::m_depos_order_xy);
2323
m_slices_nguards = {nguards_xy, nguards_xy, 0};
2424

25-
// Only xy slices need guard cells, there is no deposition to/gather from the output array F.
26-
amrex::IntVect nguards_F = amrex::IntVect(0,0,0);
27-
2825
for (int islice=0; islice<WhichSlice::N; islice++) {
2926
m_slices[lev][islice].define(
3027
slice_ba, slice_dm, Comps[islice]["N"], m_slices_nguards,
3128
amrex::MFInfo().SetArena(amrex::The_Arena()));
3229
m_slices[lev][islice].setVal(0.0);
3330
}
3431

35-
// FIXME the Poisson solver must be constructed per level, here it is only constructed for lev 0
36-
if (lev>0) return;
3732
// The Poisson solver operates on transverse slices only.
3833
// The constructor takes the BoxArray and the DistributionMap of a slice,
3934
// so the FFTPlans are built on a slice.
4035
if (m_do_dirichlet_poisson){
41-
m_poisson_solver = std::unique_ptr<FFTPoissonSolverDirichlet>(
36+
m_poisson_solver.push_back(std::unique_ptr<FFTPoissonSolverDirichlet>(
4237
new FFTPoissonSolverDirichlet(getSlices(lev, WhichSlice::This).boxArray(),
4338
getSlices(lev, WhichSlice::This).DistributionMap(),
44-
geom));
39+
geom[lev])) );
4540
} else {
46-
m_poisson_solver = std::unique_ptr<FFTPoissonSolverPeriodic>(
41+
m_poisson_solver.push_back(std::unique_ptr<FFTPoissonSolverPeriodic>(
4742
new FFTPoissonSolverPeriodic(getSlices(lev, WhichSlice::This).boxArray(),
4843
getSlices(lev, WhichSlice::This).DistributionMap(),
49-
geom));
44+
geom[lev])) );
5045
}
5146
}
5247

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

257252
// calculating the right-hand side 1/episilon0 * -(rho-Jz/c)
258-
amrex::MultiFab::Copy(m_poisson_solver->StagingArea(), getSlices(lev, WhichSlice::This),
253+
amrex::MultiFab::Copy(m_poisson_solver[lev]->StagingArea(), getSlices(lev, WhichSlice::This),
259254
Comps[WhichSlice::This]["jz"], 0, 1, 0);
260-
m_poisson_solver->StagingArea().mult(-1./phys_const.c);
261-
amrex::MultiFab::Add(m_poisson_solver->StagingArea(), getSlices(lev, WhichSlice::This),
255+
m_poisson_solver[lev]->StagingArea().mult(-1./phys_const.c);
256+
amrex::MultiFab::Add(m_poisson_solver[lev]->StagingArea(), getSlices(lev, WhichSlice::This),
262257
Comps[WhichSlice::This]["rho"], 0, 1, 0);
263-
m_poisson_solver->StagingArea().mult(-1./phys_const.ep0);
258+
m_poisson_solver[lev]->StagingArea().mult(-1./phys_const.ep0);
264259

265-
m_poisson_solver->SolvePoissonEquation(lhs);
260+
m_poisson_solver[lev]->SolvePoissonEquation(lhs);
266261

267262
/* ---------- Transverse FillBoundary Psi ---------- */
268263
amrex::ParallelContext::push(m_comm_xy);
@@ -306,7 +301,7 @@ Fields::SolvePoissonEz (amrex::Geometry const& geom, const int lev)
306301
// from the slice MF, and store in the staging area of poisson_solver
307302
TransverseDerivative(
308303
getSlices(lev, WhichSlice::This),
309-
m_poisson_solver->StagingArea(),
304+
m_poisson_solver[lev]->StagingArea(),
310305
Direction::x,
311306
geom.CellSize(Direction::x),
312307
1./(phys_const.ep0*phys_const.c),
@@ -315,7 +310,7 @@ Fields::SolvePoissonEz (amrex::Geometry const& geom, const int lev)
315310

316311
TransverseDerivative(
317312
getSlices(lev, WhichSlice::This),
318-
m_poisson_solver->StagingArea(),
313+
m_poisson_solver[lev]->StagingArea(),
319314
Direction::y,
320315
geom.CellSize(Direction::y),
321316
1./(phys_const.ep0*phys_const.c),
@@ -324,7 +319,7 @@ Fields::SolvePoissonEz (amrex::Geometry const& geom, const int lev)
324319
// Solve Poisson equation.
325320
// The RHS is in the staging area of poisson_solver.
326321
// The LHS will be returned as lhs.
327-
m_poisson_solver->SolvePoissonEquation(lhs);
322+
m_poisson_solver[lev]->SolvePoissonEquation(lhs);
328323
}
329324

330325
void
@@ -338,7 +333,7 @@ Fields::SolvePoissonBx (amrex::MultiFab& Bx_iter, amrex::Geometry const& geom, c
338333
// and store in the staging area of poisson_solver
339334
TransverseDerivative(
340335
getSlices(lev, WhichSlice::This),
341-
m_poisson_solver->StagingArea(),
336+
m_poisson_solver[lev]->StagingArea(),
342337
Direction::y,
343338
geom.CellSize(Direction::y),
344339
-phys_const.mu0,
@@ -348,7 +343,7 @@ Fields::SolvePoissonBx (amrex::MultiFab& Bx_iter, amrex::Geometry const& geom, c
348343
LongitudinalDerivative(
349344
getSlices(lev, WhichSlice::Previous1),
350345
getSlices(lev, WhichSlice::Next),
351-
m_poisson_solver->StagingArea(),
346+
m_poisson_solver[lev]->StagingArea(),
352347
geom.CellSize(Direction::z),
353348
phys_const.mu0,
354349
SliceOperatorType::Add,
@@ -357,7 +352,7 @@ Fields::SolvePoissonBx (amrex::MultiFab& Bx_iter, amrex::Geometry const& geom, c
357352
// Solve Poisson equation.
358353
// The RHS is in the staging area of poisson_solver.
359354
// The LHS will be returned as lhs.
360-
m_poisson_solver->SolvePoissonEquation(Bx_iter);
355+
m_poisson_solver[lev]->SolvePoissonEquation(Bx_iter);
361356
}
362357

363358
void
@@ -371,7 +366,7 @@ Fields::SolvePoissonBy (amrex::MultiFab& By_iter, amrex::Geometry const& geom, c
371366
// and store in the staging area of poisson_solver
372367
TransverseDerivative(
373368
getSlices(lev, WhichSlice::This),
374-
m_poisson_solver->StagingArea(),
369+
m_poisson_solver[lev]->StagingArea(),
375370
Direction::x,
376371
geom.CellSize(Direction::x),
377372
phys_const.mu0,
@@ -381,7 +376,7 @@ Fields::SolvePoissonBy (amrex::MultiFab& By_iter, amrex::Geometry const& geom, c
381376
LongitudinalDerivative(
382377
getSlices(lev, WhichSlice::Previous1),
383378
getSlices(lev, WhichSlice::Next),
384-
m_poisson_solver->StagingArea(),
379+
m_poisson_solver[lev]->StagingArea(),
385380
geom.CellSize(Direction::z),
386381
-phys_const.mu0,
387382
SliceOperatorType::Add,
@@ -390,7 +385,7 @@ Fields::SolvePoissonBy (amrex::MultiFab& By_iter, amrex::Geometry const& geom, c
390385
// Solve Poisson equation.
391386
// The RHS is in the staging area of poisson_solver.
392387
// The LHS will be returned as lhs.
393-
m_poisson_solver->SolvePoissonEquation(By_iter);
388+
m_poisson_solver[lev]->SolvePoissonEquation(By_iter);
394389
}
395390

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

417412
TransverseDerivative(
418413
getSlices(lev, WhichSlice::This),
419-
m_poisson_solver->StagingArea(),
414+
m_poisson_solver[lev]->StagingArea(),
420415
Direction::x,
421416
geom.CellSize(Direction::x),
422417
-phys_const.mu0,
@@ -425,7 +420,7 @@ Fields::SolvePoissonBz (amrex::Geometry const& geom, const int lev)
425420
// Solve Poisson equation.
426421
// The RHS is in the staging area of m_poisson_solver.
427422
// The LHS will be returned as lhs.
428-
m_poisson_solver->SolvePoissonEquation(lhs);
423+
m_poisson_solver[lev]->SolvePoissonEquation(lhs);
429424
}
430425

431426
void

0 commit comments

Comments
 (0)