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

Add leapfrog pusher for plasma #847

Merged
merged 22 commits into from
Jan 31, 2023
+627 −650
Merged

Add leapfrog pusher for plasma #847

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
22 commits
Select commit Hold shift + click to select a range
9e4fd00
leapfrog pusher
AlexanderSinn Dec 10, 2022
5ceaecd
Add second order using DualNumber Autodiff
AlexanderSinn Dec 10, 2022
daea610
include new header
AlexanderSinn Dec 10, 2022
2b3af6e
add second gather
AlexanderSinn Dec 11, 2022
87ccd5d
test with second gather
AlexanderSinn Dec 11, 2022
f07bc32
remove force terms
AlexanderSinn Dec 12, 2022
53c3f71
fix salame
AlexanderSinn Dec 12, 2022
51280dd
fix periodic issue
AlexanderSinn Dec 13, 2022
91b2d9f
leapfrog pusher with ab5 option
AlexanderSinn Jan 12, 2023
4d1e0cd
Push after ionization and collisions
AlexanderSinn Jan 12, 2023
72c62c9
fix missing variable
AlexanderSinn Jan 12, 2023
15e4cda
fix missing variable
AlexanderSinn Jan 12, 2023
ea969aa
Merge branch 'development' into leapfrog_pusher_with_ab5_option
AlexanderSinn Jan 13, 2023
3140376
merge dev
AlexanderSinn Jan 17, 2023
9685963
Merge branch 'development' into leapfrog_pusher_with_ab5_option
AlexanderSinn Jan 25, 2023
46e8d2e
push psi instead of psi_inv
AlexanderSinn Jan 25, 2023
604eee0
Add support for collisions and ionization
AlexanderSinn Jan 25, 2023
e5bf096
make leapfrog pusher the default
AlexanderSinn Jan 26, 2023
da1d688
Make AB5 pusher default and reset CI for reasons mentioned above
AlexanderSinn Jan 26, 2023
ab79215
Reset production SI benchmark manually ☹
AlexanderSinn Jan 26, 2023
730433b
Add pusher printout
AlexanderSinn Jan 30, 2023
590677d
add suggestions
AlexanderSinn Jan 31, 2023
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
11 changes: 11 additions & 0 deletions CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,13 @@ set_default_install_dirs()
option(HiPACE_MPI "Multi-node support (message-passing)" ON)
option(HiPACE_OPENPMD "openPMD I/O (HDF5, ADIOS)" ON)

set(HiPACE_PUSHER_VALUES LEAPFROG AB5)
set(HiPACE_PUSHER AB5 CACHE STRING "Plasma pusher (AB5/LEAPFROG)")
set_property(CACHE HiPACE_PUSHER PROPERTY STRINGS ${HiPACE_PUSHER_VALUES})
if(NOT HiPACE_PUSHER IN_LIST HiPACE_PUSHER_VALUES)
message(FATAL_ERROR "HiPACE_PUSHER (${HiPACE_PUSHER}) must be one of ${HiPACE_PUSHER_VALUES}")
endif()

set(HiPACE_PRECISION_VALUES SINGLE DOUBLE)
set(HiPACE_PRECISION DOUBLE CACHE STRING "Floating point precision (SINGLE/DOUBLE)")
set_property(CACHE HiPACE_PRECISION PROPERTY STRINGS ${HiPACE_PRECISION_VALUES})
Expand Down Expand Up @@ -155,6 +162,10 @@ if(HiPACE_OPENPMD)
target_link_libraries(HiPACE PUBLIC openPMD::openPMD)
endif()

if(HiPACE_PUSHER STREQUAL "AB5")
target_compile_definitions(HiPACE PUBLIC HIPACE_USE_AB5_PUSH)
endif()

if(AMReX_LINEAR_SOLVERS)
target_compile_definitions(HiPACE PUBLIC AMREX_USE_LINEAR_SOLVERS)
endif()
Expand Down
7 changes: 4 additions & 3 deletions docs/source/building/building.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -154,9 +154,9 @@ or by providing arguments to the CMake call
cmake -S . -B build -D<OPTION_A>=<VALUE_A> -D<OPTION_B>=<VALUE_B>


============================= ======================================== =====================================================
============================= ======================================== =========================================================
CMake Option Default & Values Description
----------------------------- ---------------------------------------- -----------------------------------------------------
----------------------------- ---------------------------------------- ---------------------------------------------------------
``CMAKE_BUILD_TYPE`` RelWithDebInfo/**Release**/Debug Type of build, symbols & optimizations
``HiPACE_COMPUTE`` NOACC/CUDA/SYCL/HIP/**OMP** On-node, accelerated computing backend
``HiPACE_MPI`` **ON**/OFF Multi-node support (message-passing)
Expand All @@ -165,7 +165,8 @@ or by providing arguments to the CMake call
``HiPACE_amrex_branch`` ``development`` Repository branch for ``HiPACE_amrex_repo``
``HiPACE_amrex_internal`` **ON**/OFF Needs a pre-installed AMReX library if set to ``OFF``
``HiPACE_OPENPMD`` **ON**/OFF openPMD I/O (HDF5, ADIOS2)
============================= ======================================== =====================================================
``HiPACE_PUSHER`` **AB5**/LEAPFROG Use fifth-order Adams-Bashforth or leapfrog Plasma Pusher
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@AlexanderSinn is it work making this a compile-time option?

I think pre-compiling both options and doing a runtime switch before the two specialized kernels would be more user friendly.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Discussed: needs additional (25) force terms on the particle, so that's a bit of a heavy lift.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As mentioned, it is a compile-time option not for performance of the push kernel, but rather to save on memory of the 25 Adams-Bashforth force terms (per plasma particle) if the leapfrog pusher is used. One would need to change these to runtime components to make this a runtime option.

============================= ======================================== =========================================================

HiPACE++ can be configured in further detail with options from AMReX, which are documented in the `AMReX manual <https://amrex-codes.github.io/amrex/docs_html/BuildingAMReX.html#customization-options>`__.

Expand Down
59 changes: 23 additions & 36 deletions src/Hipace.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -297,7 +297,11 @@ Hipace::InitData ()
amrex::Print() << "using CUDA version " << __CUDACC_VER_MAJOR__ << "." << __CUDACC_VER_MINOR__
<< "." << __CUDACC_VER_BUILD__ << "\n";
#endif

#ifdef HIPACE_USE_AB5_PUSH
amrex::Print() << "using the Adams-Bashforth plasma particle pusher\n";
#else
amrex::Print() << "using the leapfrog plasma particle pusher\n";
#endif

amrex::Vector<amrex::IntVect> new_max_grid_size;
for (int ilev = 0; ilev <= maxLevel(); ++ilev) {
Expand Down Expand Up @@ -462,7 +466,7 @@ Hipace::Evolve ()
// Only reset plasma after receiving time step, to use proper density
// WARNING: handling of lev is to be improved: this loops over levels, but
// lev is set to 0 above.
for (int lv=0; lv<=finestLevel(); ++lv) m_multi_plasma.ResetParticles(lv, true);
for (int lv=0; lv<=finestLevel(); ++lv) m_multi_plasma.ResetParticles(lv);
/* Store charge density of (immobile) ions into WhichSlice::RhoIons */
if (m_do_tiling) m_multi_plasma.TileSort(boxArray(lev)[0], geom[lev]);
m_multi_plasma.DepositNeutralizingBackground(
Expand Down Expand Up @@ -597,12 +601,19 @@ Hipace::SolveOneSlice (int islice_coarse, const int ibox, int step,

FillDiagnostics(lev, islice);

m_multi_plasma.doCoulombCollision(lev, bx, geom[lev]);

m_multi_plasma.DoFieldIonization(lev, geom[lev], m_fields);

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

// Push plasma particles
m_multi_plasma.AdvanceParticles(m_fields, m_multi_laser, geom[lev], false, lev);

m_multi_plasma.doCoulombCollision(lev, bx, geom[lev]);

// Push beam particles
m_multi_beam.AdvanceBeamParticlesSlice(m_fields, geom[lev], lev, islice_local, bx,
bins[lev], m_box_sorters, ibox);

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);
Expand Down Expand Up @@ -641,7 +652,7 @@ Hipace::ExplicitSolveOneSubSlice (const int lev, const int step, const int ibox,

// deposit jx, jy, jz, rho and chi for all plasmas
m_multi_plasma.DepositCurrent(
m_fields, m_multi_laser, WhichSlice::This, false, true, true, true, true, geom[lev], lev);
m_fields, m_multi_laser, WhichSlice::This, true, true, true, true, geom[lev], lev);

m_fields.setVal(0., lev, WhichSlice::Next, "jx_beam", "jy_beam");
// deposit jx_beam and jy_beam in the Next slice
Expand Down Expand Up @@ -675,23 +686,14 @@ Hipace::ExplicitSolveOneSubSlice (const int lev, const int step, const int ibox,
// Solves Bx, By using Sx, Sy and chi
ExplicitMGSolveBxBy(lev, WhichSlice::This, islice);

const bool do_salame = m_multi_beam.isSalameNow(step, islice_local, beam_bin);
if (do_salame) {
if (m_multi_beam.isSalameNow(step, islice_local, beam_bin)) {
// Modify the beam particle weights on this slice to flatten Ez.
// As the beam current is modified, Bx and By are also recomputed.
// Plasma particle force terms get shifted
SalameModule(this, m_salame_n_iter, m_salame_do_advance, m_salame_last_slice,
m_salame_overloaded, lev, step, islice, islice_local, beam_bin, ibox);
}

// shift and update force terms, push plasma particles
// don't shift force terms again if salame was used
m_multi_plasma.AdvanceParticles(m_fields, m_multi_laser, geom[lev],
false, true, true, !do_salame, lev);

// Push beam particles
m_multi_beam.AdvanceBeamParticlesSlice(m_fields, geom[lev], lev, islice_local, bx,
beam_bin, m_box_sorters, ibox);
// Push beam and plasma in SolveOneSlice
}

void
Expand All @@ -706,14 +708,11 @@ Hipace::PredictorCorrectorSolveOneSubSlice (const int lev, const int step, const
m_fields.setVal(0., lev, WhichSlice::This, "chi");
}

// push plasma
m_multi_plasma.AdvanceParticles(m_fields, m_multi_laser, geom[lev], false,
true, false, false, lev);

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

// deposit jx jy jz rho and maybe chi
m_multi_plasma.DepositCurrent(
m_fields, m_multi_laser, WhichSlice::This, false, true, true, true, m_use_laser, geom[lev], lev);
m_fields, m_multi_laser, WhichSlice::This, true, true, true, m_use_laser, geom[lev], lev);

m_fields.AddRhoIons(lev);

Expand Down Expand Up @@ -743,9 +742,7 @@ Hipace::PredictorCorrectorSolveOneSubSlice (const int lev, const int step, const
// Solves Bx and By in the current slice and modifies the force terms of the plasma particles
PredictorCorrectorLoopToSolveBxBy(islice_local, lev, step, beam_bin, ibox);

// Push beam particles
m_multi_beam.AdvanceBeamParticlesSlice(m_fields, geom[lev], lev, islice_local, bx,
beam_bin, m_box_sorters, ibox);
// Push beam and plasma in SolveOneSlice
}

void
Expand Down Expand Up @@ -995,9 +992,6 @@ Hipace::PredictorCorrectorLoopToSolveBxBy (const int islice_local, const int lev
amrex::MultiFab::Copy(By_prev_iter, m_fields.getSlices(lev, WhichSlice::This),
Comps[WhichSlice::This]["By"], 0, 1, m_fields.m_slices_nguards);

// shift force terms, update force terms using guessed Bx and By
m_multi_plasma.AdvanceParticles(m_fields, m_multi_laser, geom[lev], false, false, true, true, lev);

const int islice = islice_local + boxArray(lev)[ibox].smallEnd(Direction::z);

// Begin of predictor corrector loop
Expand All @@ -1011,12 +1005,12 @@ Hipace::PredictorCorrectorLoopToSolveBxBy (const int islice_local, const int lev
m_predcorr_avg_iterations += 1.0;

// Push particles to the next temp slice
m_multi_plasma.AdvanceParticles(m_fields, m_multi_laser, geom[lev], true, true, false, false, lev);
m_multi_plasma.AdvanceParticles(m_fields, m_multi_laser, geom[lev], true, lev);

if (m_do_tiling) m_multi_plasma.TileSort(boxArray(lev)[0], geom[lev]);
// plasmas deposit jx jy to next temp slice
m_multi_plasma.DepositCurrent(
m_fields, m_multi_laser, WhichSlice::Next, true, true, false, false, false, geom[lev], lev);
m_fields, m_multi_laser, WhichSlice::Next, true, false, false, false, geom[lev], lev);

// beams deposit jx jy to the next slice
m_multi_beam.DepositCurrentSlice(m_fields, geom, lev, step, islice_local, bins,
Expand Down Expand Up @@ -1062,13 +1056,6 @@ Hipace::PredictorCorrectorLoopToSolveBxBy (const int islice_local, const int lev
amrex::ParallelContext::pop();
}

// resetting the particle position after they have been pushed to the next slice
m_multi_plasma.ResetParticles(lev);

// Update force terms using the calculated Bx and By
m_multi_plasma.AdvanceParticles(m_fields, m_multi_laser, geom[lev],
false, false, true, false, lev);

// Shift relative_Bfield_error values
relative_Bfield_error_prev_iter = relative_Bfield_error;
} /* end of predictor corrector loop */
Expand Down
49 changes: 49 additions & 0 deletions src/fields/Fields.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -436,6 +436,55 @@ public:
}
}

/** \brief add all selected fields between slices or on the same slice
* Uses references to C-arrays as argument so that ncomps can be deduced
*
* \param[in] lev level of fields
* \param[in] islice_dst destination slice
* \param[in] comps_dst array of destination component names
* \param[in] islice_src source slice
* \param[in] comps_src array of source component names
*/
template<int ncomps>
void add (const int lev,
const int islice_dst, const char * const (&comps_dst)[ncomps],
const int islice_src, const char * const (&comps_src)[ncomps]) {
amrex::GpuArray<int, ncomps> c_idx_src = {};
amrex::GpuArray<int, ncomps> c_idx_dst = {};
for (int i=0; i<ncomps; ++i) {
c_idx_src[i] = Comps[islice_src][comps_src[i]];
c_idx_dst[i] = Comps[islice_dst][comps_dst[i]];
}

const amrex::MultiFab& mfab_src = getSlices(lev, islice_src);
amrex::MultiFab& mfab_dst = getSlices(lev, islice_dst);

for (amrex::MFIter mfi(mfab_dst, DfltMfi); mfi.isValid(); ++mfi) {
const Array3<amrex::Real> array_dst = mfab_dst.array(mfi);
// only one Kernel for all fields
if (islice_dst != islice_src) {
const Array3<amrex::Real const> array_src = mfab_src.const_array(mfi);
amrex::ParallelFor(mfi.growntilebox(),
[=] AMREX_GPU_DEVICE(int i, int j, int)
{
for (int n=0; n<ncomps; ++n) {
array_dst(i,j,c_idx_dst[n]) += array_src(i,j,c_idx_src[n]);
}
});
} else {
amrex::ParallelFor(mfi.growntilebox(),
[=] AMREX_GPU_DEVICE(int i, int j, int)
{
for (int n=0; n<ncomps; ++n) {
// array_src == array_dst here
array_dst(i,j,c_idx_dst[n]) += array_dst(i,j,c_idx_src[n]);
}
});
}

}
}

/** \brief call amrex FillBoundary for multiple fields
*
* \param[in] period periodicity of the fields
Expand Down
2 changes: 1 addition & 1 deletion src/fields/Fields.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,7 +102,7 @@ Fields::AllocData (
isl = WhichSlice::Salame;
if (any_salame) {
Comps[isl].multi_emplace(N_Comps[isl], "Ez_target", "Ez_no_salame", "Ez",
"jx", "jy", "jz_beam", "Bx", "By", "Sy", "Sx");
"jx", "jy", "jz_beam", "Bx", "By", "Sy", "Sx", "Sy_back", "Sx_back");
}

} else {
Expand Down
18 changes: 9 additions & 9 deletions src/particles/collisions/CoulombCollision.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,9 +57,9 @@ CoulombCollision::doCoulombCollision (

// Get particles SoA and AoS data
auto& soa1 = pti.GetStructOfArrays();
amrex::Real* const ux1 = soa1.GetRealData(PlasmaIdx::ux).data();
amrex::Real* const uy1 = soa1.GetRealData(PlasmaIdx::uy).data();
amrex::Real* const psi1 = soa1.GetRealData(PlasmaIdx::psi).data();
amrex::Real* const ux1 = soa1.GetRealData(PlasmaIdx::ux_half_step).data();
amrex::Real* const uy1 = soa1.GetRealData(PlasmaIdx::uy_half_step).data();
amrex::Real* const psi1 = soa1.GetRealData(PlasmaIdx::psi_half_step).data();
const amrex::Real* const w1 = soa1.GetRealData(PlasmaIdx::w).data();
PlasmaBins::index_type * const indices1 = bins1.permutationPtr();
PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr();
Expand Down Expand Up @@ -114,9 +114,9 @@ CoulombCollision::doCoulombCollision (

// Get particles SoA and AoS data for species 1
auto& soa1 = pti.GetStructOfArrays();
amrex::Real* const ux1 = soa1.GetRealData(PlasmaIdx::ux).data();
amrex::Real* const uy1 = soa1.GetRealData(PlasmaIdx::uy).data();
amrex::Real* const psi1 = soa1.GetRealData(PlasmaIdx::psi).data();
amrex::Real* const ux1 = soa1.GetRealData(PlasmaIdx::ux_half_step).data();
amrex::Real* const uy1 = soa1.GetRealData(PlasmaIdx::uy_half_step).data();
amrex::Real* const psi1 = soa1.GetRealData(PlasmaIdx::psi_half_step).data();
const amrex::Real* const w1 = soa1.GetRealData(PlasmaIdx::w).data();
PlasmaBins::index_type * const indices1 = bins1.permutationPtr();
PlasmaBins::index_type const * const offsets1 = bins1.offsetsPtr();
Expand All @@ -126,9 +126,9 @@ CoulombCollision::doCoulombCollision (
// Get particles SoA and AoS data for species 2
auto& ptile2 = species2.ParticlesAt(lev, pti.index(), pti.LocalTileIndex());
auto& soa2 = ptile2.GetStructOfArrays();
amrex::Real* const ux2 = soa2.GetRealData(PlasmaIdx::ux).data();
amrex::Real* const uy2 = soa2.GetRealData(PlasmaIdx::uy).data();
amrex::Real* const psi2= soa2.GetRealData(PlasmaIdx::psi).data();
amrex::Real* const ux2 = soa2.GetRealData(PlasmaIdx::ux_half_step).data();
amrex::Real* const uy2 = soa2.GetRealData(PlasmaIdx::uy_half_step).data();
amrex::Real* const psi2= soa2.GetRealData(PlasmaIdx::psi_half_step).data();
const amrex::Real* const w2 = soa2.GetRealData(PlasmaIdx::w).data();
PlasmaBins::index_type * const indices2 = bins2.permutationPtr();
PlasmaBins::index_type const * const offsets2 = bins2.offsetsPtr();
Expand Down
3 changes: 1 addition & 2 deletions src/particles/deposition/PlasmaDepositCurrent.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,6 @@
* \param[in,out] fields the general field class, modified by this function
* \param[in] multi_laser MultiLaser that affects the plasma during the deposition
* \param[in] which_slice defines if this or the next slice is handled
* \param[in] temp_slice if true, the temporary data (x_temp, ...) is used
* \param[in] deposit_jx_jy if true, deposit to jx and jy
* \param[in] deposit_jz if true, deposit to jz
* \param[in] deposit_rho if true, deposit to rho
Expand All @@ -32,7 +31,7 @@
*/
void
DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const MultiLaser& multi_laser,
const int which_slice, const bool temp_slice,
const int which_slice,
const bool deposit_jx_jy, const bool deposit_jz, const bool deposit_rho,
bool deposit_chi, amrex::Geometry const& gm, int const lev,
const PlasmaBins& bins, int bin_size);
Expand Down
11 changes: 4 additions & 7 deletions src/particles/deposition/PlasmaDepositCurrent.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@

void
DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const MultiLaser& multi_laser,
const int which_slice, const bool temp_slice,
const int which_slice,
const bool deposit_jx_jy, const bool deposit_jz, const bool deposit_rho,
const bool deposit_chi, amrex::Geometry const& gm, int const lev,
const PlasmaBins& bins, int bin_size)
Expand Down Expand Up @@ -74,12 +74,9 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const MultiLas
auto& soa = pti.GetStructOfArrays(); // For momenta and weights

amrex::Real * const AMREX_RESTRICT wp = soa.GetRealData(PlasmaIdx::w).data();
const amrex::Real * const AMREX_RESTRICT uxp =
soa.GetRealData(temp_slice ? PlasmaIdx::ux_temp : PlasmaIdx::ux).data();
const amrex::Real * const AMREX_RESTRICT uyp =
soa.GetRealData(temp_slice ? PlasmaIdx::uy_temp : PlasmaIdx::uy).data();
const amrex::Real * const AMREX_RESTRICT psip =
soa.GetRealData(temp_slice ? PlasmaIdx::psi_temp : PlasmaIdx::psi).data();
const amrex::Real * const AMREX_RESTRICT uxp = soa.GetRealData(PlasmaIdx::ux).data();
const amrex::Real * const AMREX_RESTRICT uyp = soa.GetRealData(PlasmaIdx::uy).data();
const amrex::Real * const AMREX_RESTRICT psip = soa.GetRealData(PlasmaIdx::psi).data();

int const * const AMREX_RESTRICT a_ion_lev =
plasma.m_can_ionize ? soa.GetIntData(PlasmaIdx::ion_lev).data() : nullptr;
Expand Down
Loading