From a26b66ebe505fd71e2c53fc7f2d9d242f94d55cf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maxence=20Th=C3=A9venet?= Date: Thu, 24 Nov 2022 17:55:51 +0100 Subject: [PATCH] Enable multiple laser pulses (#823) * Rename laser multilaser * option to have two lasers on 1 grid * fix input files to specify lasers * eol * init lasers * fix * safer for GPU when laser class becomes bigger, although a few more kernel launches * try fix Juwels compile issue * update doc for multiple laser pulses Co-authored-by: AlexanderSinn --- docs/source/run/parameters.rst | 66 +- examples/beam_in_vacuum/inputs_SI | 2 + examples/beam_in_vacuum/inputs_normalized | 2 + .../inputs_normalized_transverse | 3 + examples/blowout_wake/inputs_SI | 2 + examples/blowout_wake/inputs_ionization_SI | 2 + examples/blowout_wake/inputs_normalized | 2 + examples/gaussian_weight/inputs_SI | 2 + examples/gaussian_weight/inputs_normalized | 2 + examples/get_started/inputs_normalized | 4 + examples/laser/inputs_SI | 7 +- examples/linear_wake/inputs_SI | 2 + examples/linear_wake/inputs_ion_motion_SI | 2 + examples/linear_wake/inputs_normalized | 2 + src/Hipace.H | 4 +- src/Hipace.cpp | 62 +- src/fields/Fields.H | 8 +- src/fields/Fields.cpp | 6 +- src/laser/CMakeLists.txt | 1 + src/laser/Laser.H | 237 +---- src/laser/Laser.cpp | 762 +--------------- src/laser/MultiLaser.H | 244 +++++ src/laser/MultiLaser.cpp | 830 ++++++++++++++++++ src/particles/deposition/ExplicitDeposition.H | 4 +- .../deposition/ExplicitDeposition.cpp | 8 +- .../deposition/PlasmaDepositCurrent.H | 6 +- .../deposition/PlasmaDepositCurrent.cpp | 8 +- src/particles/plasma/MultiPlasma.H | 18 +- src/particles/plasma/MultiPlasma.cpp | 22 +- src/particles/pusher/PlasmaParticleAdvance.H | 6 +- .../pusher/PlasmaParticleAdvance.cpp | 6 +- src/salame/Salame.cpp | 8 +- tests/laser_blowout_wake_explicit.1Rank.sh | 6 +- tests/laser_blowout_wake_explicit.SI.1Rank.sh | 6 +- tests/laser_evolution.SI.2Rank.sh | 4 +- 35 files changed, 1278 insertions(+), 1078 deletions(-) create mode 100644 src/laser/MultiLaser.H create mode 100644 src/laser/MultiLaser.cpp diff --git a/docs/source/run/parameters.rst b/docs/source/run/parameters.rst index 1da720eb04..835e4ebe3b 100644 --- a/docs/source/run/parameters.rst +++ b/docs/source/run/parameters.rst @@ -508,57 +508,61 @@ Laser parameters ---------------- The laser profile is defined by :math:`a(x,y,z) = a_0 * \mathrm{exp}[-(x^2/w0_x^2 + y^2/w0_y^2 + z^2/L0^2)]`. -The laser pulse length :math:`L0 = \tau / c_0` can be specified via the pulse duration ``laser.tau``. The model implemented is the one from [C. Benedetti et al. Plasma Phys. Control. Fusion 60.1: 014002 (2017)]. +The model implemented is the one from [C. Benedetti et al. Plasma Phys. Control. Fusion 60.1: 014002 (2017)]. +Unlike for ``beams`` and ``plasmas``, all the laser pulses are currently stored on the same array, +which you can find in the output openPMD file as `laser_real` (for the real part of the envelope) and `laser_imag` for its imaginary part. +Parameters starting with ``lasers.`` apply to all laser pulses, parameters starting with ```` apply to a single laser pulse. -* ``laser.use_laser`` (`0` or `1`) optional (default `0`) - Whether to activate the laser envelope solver. +* ``lasers.names`` (list of `string`) + The names of the laser pulses, separated by a space. + To run without a laser, choose the name ``no_laser``. -* ``laser.a0`` (`float`) optional (default `0`) - Peak normalized vector potential of the laser pulse. - -* ``laser.position_mean`` (3 `float`) optional (default `0 0 0`) - The mean position of the laser in `x, y, z`. - -* ``laser.w0`` (2 `float`) optional (default `0 0`) - The laser waist in `x, y`. - -* ``laser.L0`` (`float`) optional (default `0`) - The laser pulse length in `z`. Use either the pulse length or the pulse duration. - -* ``laser.tau`` (`float`) optional (default `0`) - The laser pulse duration. The pulse length will be set to `laser.tau`:math:`/c_0`. - Use either the pulse length or the pulse duration. - -* ``laser.lambda0`` (`float`) - The laser pulse wavelength. - -* ``laser.focal_distance`` (`float`) - Distance at which the laser pulse if focused (in the z direction, counted from laser initial position). +* ``lasers.lambda0`` (`float`) + Wavelength of the laser pulses. Currently, all pulses must have the same wavelength. -* ``laser.use_phase`` (`bool`) optional (default `true`) +* ``lasers.use_phase`` (`bool`) optional (default `true`) Whether the phase terms (:math:`\theta` in Eq. (6) of [C. Benedetti et al. Plasma Phys. Control. Fusion 60.1: 014002 (2017)]) are computed and used in the laser envelope advance. Keeping the phase should be more accurate, but can cause numerical issues in the presence of strong depletion/frequency shift. -* ``laser.solver_type`` (`string`) optional (default `multigrid`) +* ``lasers.solver_type`` (`string`) optional (default `multigrid`) Type of solver for the laser envelope solver, either ``fft`` or ``multigrid``. Currently, the approximation that the phase is evaluated on-axis only is made with both solvers. With the multigrid solver, we could drop this assumption. For now, the fft solver should be faster, more accurate and more stable, so only use the multigrid one with care. -* ``laser.MG_tolerance_rel`` (`float`) optional (default `1e-4`) +* ``lasers.MG_tolerance_rel`` (`float`) optional (default `1e-4`) Relative error tolerance of the multigrid solver used for the laser pulse. -* ``laser.MG_tolerance_abs`` (`float`) optional (default `0.`) +* ``lasers.MG_tolerance_abs`` (`float`) optional (default `0.`) Absolute error tolerance of the multigrid solver used for the laser pulse. -* ``laser.MG_verbose`` (`int`) optional (default `0`) +* ``lasers.MG_verbose`` (`int`) optional (default `0`) Level of verbosity of the multigrid solver used for the laser pulse. -* ``laser.MG_average_rhs`` (`0` or `1`) optional (default `1`) +* ``lasers.MG_average_rhs`` (`0` or `1`) optional (default `1`) Whether to use the most stable discretization for the envelope solver. -* ``laser.3d_on_host`` (`0` or `1`) optional (default `0`) +* ``lasers.3d_on_host`` (`0` or `1`) optional (default `0`) When running on GPU: whether the 3D array containing the laser envelope is stored in host memory (CPU, slower but large memory available) or in device memory (GPU, faster but less memory available). +* ``.a0`` (`float`) optional (default `0`) + Peak normalized vector potential of the laser pulse. + +* ``.position_mean`` (3 `float`) optional (default `0 0 0`) + The mean position of the laser in `x, y, z`. + +* ``.w0`` (2 `float`) optional (default `0 0`) + The laser waist in `x, y`. + +* ``.L0`` (`float`) optional (default `0`) + The laser pulse length in `z`. Use either the pulse length or the pulse duration ``.tau``. + +* ``.tau`` (`float`) optional (default `0`) + The laser pulse duration. The pulse length will be set to `laser.tau`:math:`/c_0`. + Use either the pulse length or the pulse duration. + +* ``.focal_distance`` (`float`) + Distance at which the laser pulse if focused (in the z direction, counted from laser initial position). + Diagnostic parameters --------------------- diff --git a/examples/beam_in_vacuum/inputs_SI b/examples/beam_in_vacuum/inputs_SI index d93440fbd1..09490b5598 100644 --- a/examples/beam_in_vacuum/inputs_SI +++ b/examples/beam_in_vacuum/inputs_SI @@ -39,4 +39,6 @@ beam2.ppc = 2 2 1 plasmas.names = no_plasma +lasers.names = no_laser + diagnostic.diag_type = xyz diff --git a/examples/beam_in_vacuum/inputs_normalized b/examples/beam_in_vacuum/inputs_normalized index 13faea779e..ce204d208d 100644 --- a/examples/beam_in_vacuum/inputs_normalized +++ b/examples/beam_in_vacuum/inputs_normalized @@ -30,4 +30,6 @@ beam.ppc = 2 2 1 plasmas.names = no_plasma +lasers.names = no_laser + diagnostic.diag_type = xyz diff --git a/examples/beam_in_vacuum/inputs_normalized_transverse b/examples/beam_in_vacuum/inputs_normalized_transverse index 6a1eea4f04..f18f35e4e5 100644 --- a/examples/beam_in_vacuum/inputs_normalized_transverse +++ b/examples/beam_in_vacuum/inputs_normalized_transverse @@ -39,4 +39,7 @@ beam2.position_mean = 0. 0. 0 beam2.position_std = 8. 0.3 1.41 plasmas.names = no_plasma + +lasers.names = no_laser + diagnostic.diag_type = xz diff --git a/examples/blowout_wake/inputs_SI b/examples/blowout_wake/inputs_SI index fb69a21500..4c2bb62748 100644 --- a/examples/blowout_wake/inputs_SI +++ b/examples/blowout_wake/inputs_SI @@ -44,4 +44,6 @@ plasma.ppc = 1 1 plasma.u_mean = 0.0 0.0 0. plasma.element = electron +lasers.names = no_laser + diagnostic.diag_type = xyz diff --git a/examples/blowout_wake/inputs_ionization_SI b/examples/blowout_wake/inputs_ionization_SI index 5611476255..0d8027a586 100644 --- a/examples/blowout_wake/inputs_ionization_SI +++ b/examples/blowout_wake/inputs_ionization_SI @@ -53,4 +53,6 @@ ion.mass_Da = 1.008 ion.initial_ion_level = 0 ion.ionization_product = elec +lasers.names = no_laser + diagnostic.diag_type = xyz diff --git a/examples/blowout_wake/inputs_normalized b/examples/blowout_wake/inputs_normalized index c973902a3e..d54837479f 100644 --- a/examples/blowout_wake/inputs_normalized +++ b/examples/blowout_wake/inputs_normalized @@ -40,4 +40,6 @@ plasma.ppc = 1 1 plasma.u_mean = 0.0 0.0 0. plasma.element = electron +lasers.names = no_laser + diagnostic.diag_type = xyz diff --git a/examples/gaussian_weight/inputs_SI b/examples/gaussian_weight/inputs_SI index 8f81250191..d3378f752f 100644 --- a/examples/gaussian_weight/inputs_SI +++ b/examples/gaussian_weight/inputs_SI @@ -29,4 +29,6 @@ beam.ppc = 2 2 1 plasmas.names = no_plasma +lasers.names = no_laser + diagnostic.diag_type = xyz diff --git a/examples/gaussian_weight/inputs_normalized b/examples/gaussian_weight/inputs_normalized index 359d08c799..bab5fb6695 100644 --- a/examples/gaussian_weight/inputs_normalized +++ b/examples/gaussian_weight/inputs_normalized @@ -29,4 +29,6 @@ beam.u_std = 3. 4. 0. plasmas.names = no_plasma +lasers.names = no_laser + diagnostic.diag_type = xyz diff --git a/examples/get_started/inputs_normalized b/examples/get_started/inputs_normalized index 6a1abfe1df..0b3f3fdb7b 100644 --- a/examples/get_started/inputs_normalized +++ b/examples/get_started/inputs_normalized @@ -36,6 +36,10 @@ plasma.element = electron # type of plasma: electron, proton, or a plasma.density(x,y,z) = 1. # density function in the respective unit systems plasma.ppc = 1 1 # particles per cell in x,y +### Laser +################################ +lasers.names = no_laser # name(s) of the laser pulses. Here, no laser pulse used. + ### Diagnostics ################################ diagnostic.diag_type = xz # 2D xz slice output. Options: xyz, xz, yz diff --git a/examples/laser/inputs_SI b/examples/laser/inputs_SI index 4fb6bd71eb..c4c3db6723 100644 --- a/examples/laser/inputs_SI +++ b/examples/laser/inputs_SI @@ -12,12 +12,13 @@ hipace.do_tiling = 0 geometry.prob_lo = -6.*kp_inv -6.*kp_inv -8.*kp_inv geometry.prob_hi = 6.*kp_inv 6.*kp_inv 6.*kp_inv -laser.use_laser = 1 +lasers.names = laser +lasers.lambda0 = .8e-6 + laser.a0 = 1 laser.position_mean = 0. 0. 0 -laser.w0 = 2.*kp_inv 2.*kp_inv +laser.w0 = 2.*kp_inv laser.L0 = 2.*kp_inv -laser.lambda0 = .8e-6 laser.focal_distance = 0.001 hipace.bxby_solver=explicit diff --git a/examples/linear_wake/inputs_SI b/examples/linear_wake/inputs_SI index c5fbe65c6b..f064340bca 100644 --- a/examples/linear_wake/inputs_SI +++ b/examples/linear_wake/inputs_SI @@ -38,4 +38,6 @@ plasma.ppc = 1 1 plasma.u_mean = 0.0 0.0 0. plasma.element = electron +lasers.names = no_laser + diagnostic.diag_type = xyz diff --git a/examples/linear_wake/inputs_ion_motion_SI b/examples/linear_wake/inputs_ion_motion_SI index ee1574ad14..ba5a8f201f 100644 --- a/examples/linear_wake/inputs_ion_motion_SI +++ b/examples/linear_wake/inputs_ion_motion_SI @@ -57,4 +57,6 @@ ions.neutralize_background = false "ions.density(x,y,z)" = ne ions.ppc = 1 1 +lasers.names = no_laser + diagnostic.diag_type = xyz diff --git a/examples/linear_wake/inputs_normalized b/examples/linear_wake/inputs_normalized index 52be059a4c..ac442354f1 100644 --- a/examples/linear_wake/inputs_normalized +++ b/examples/linear_wake/inputs_normalized @@ -36,3 +36,5 @@ plasma.ppc = 1 1 plasma.u_mean = 0.0 0.0 0. plasma.element = electron diagnostic.diag_type = xyz + +lasers.names = no_laser diff --git a/src/Hipace.H b/src/Hipace.H index 9321646bc4..0adb132d93 100644 --- a/src/Hipace.H +++ b/src/Hipace.H @@ -16,7 +16,7 @@ #include "particles/beam/BeamParticleContainer.H" #include "utils/AdaptiveTimeStep.H" #include "utils/GridCurrent.H" -#include "laser/Laser.H" +#include "laser/MultiLaser.H" #include "utils/Constants.H" #include "utils/Parser.H" #include "diagnostics/Diagnostic.H" @@ -361,7 +361,7 @@ public: /** Adaptive time step instance */ AdaptiveTimeStep m_adaptive_time_step; /** Laser instance (soon to be multi laser container) */ - Laser m_laser; + MultiLaser m_multi_laser; /** GridCurrent instance */ GridCurrent m_grid_current; #ifdef HIPACE_USE_OPENPMD diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 3a7edbf498..f9a7194ce4 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -99,7 +99,7 @@ Hipace::Hipace () : m_multi_beam(this), m_multi_plasma(this), m_adaptive_time_step(m_multi_beam.get_nbeams()), - m_laser(), + m_multi_laser(), m_diags(this->maxLevel()+1) { amrex::ParmParse pp;// Traditionally, max_step and stop_time do not have prefix. @@ -189,7 +189,7 @@ Hipace::Hipace () : MPI_Comm_split(amrex::ParallelDescriptor::Communicator(), m_rank_xy, myproc, &m_comm_z); #endif - m_use_laser = m_laser.m_use_laser; + m_use_laser = m_multi_laser.m_use_laser; } Hipace::~Hipace () @@ -347,8 +347,8 @@ Hipace::MakeNewLevelFromScratch ( DefineSliceGDB(lev, ba, dm); m_fields.AllocData(lev, Geom(), m_slice_ba[lev], m_slice_dm[lev], m_multi_plasma.m_sort_bin_size); - m_laser.InitData(m_slice_ba[0], m_slice_dm[0]); // laser inits only on level 0 - m_diags.Initialize(lev, m_laser.m_use_laser); + m_multi_laser.InitData(m_slice_ba[0], m_slice_dm[0]); // laser inits only on level 0 + m_diags.Initialize(lev, m_multi_laser.m_use_laser); } void @@ -436,7 +436,7 @@ Hipace::Evolve () /* 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(m_fields, m_laser, WhichSlice::RhoIons, geom[lev], + m_multi_plasma.DepositNeutralizingBackground(m_fields, m_multi_laser, WhichSlice::RhoIons, geom[lev], finestLevel()+1); // Loop over longitudinal boxes on this rank, from head to tail @@ -445,11 +445,11 @@ Hipace::Evolve () { const amrex::Box& bx = boxArray(lev)[it]; - if (m_laser.m_use_laser) { + if (m_multi_laser.m_use_laser) { AMREX_ALWAYS_ASSERT(!m_adaptive_time_step.m_do_adaptive_time_step); AMREX_ALWAYS_ASSERT(m_multi_plasma.GetNPlasmas() <= 1); // Before that, the 3D fields of the envelope are not initialized (not even allocated). - m_laser.Init3DEnvelope(step, bx, Geom(0)); + m_multi_laser.Init3DEnvelope(step, bx, Geom(0)); if (it == n_boxes-1) ResetLaser(); } @@ -547,7 +547,7 @@ Hipace::SolveOneSlice (int islice_coarse, const int ibox, int step, boxArray(0)[ibox].smallEnd(Direction::z), m_box_sorters, ibox); // Get this laser slice from the 3D array - m_laser.Copy(islice_coarse, false); + m_multi_laser.Copy(islice_coarse, false); for (int lev = 0; lev <= finestLevel(); ++lev) { @@ -591,8 +591,8 @@ Hipace::SolveOneSlice (int islice_coarse, const int ibox, int step, } // end for (int isubslice = nsubslice-1; isubslice >= 0; --isubslice) // Advance laser slice by 1 step and store result to 3D array - m_laser.AdvanceSlice(m_fields, Geom(0), m_dt, step); - m_laser.Copy(islice_coarse, true); + m_multi_laser.AdvanceSlice(m_fields, Geom(0), m_dt, step); + m_multi_laser.Copy(islice_coarse, true); // After this, the parallel context is the full 3D communicator again amrex::ParallelContext::pop(); @@ -619,7 +619,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_laser, WhichSlice::This, false, true, true, true, true, geom[lev], lev); + m_fields, m_multi_laser, WhichSlice::This, false, 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 @@ -648,7 +648,7 @@ Hipace::ExplicitSolveOneSubSlice (const int lev, const int step, const int ibox, InitializeSxSyWithBeam(lev); // Deposit Sx and Sy for every plasma species - m_multi_plasma.ExplicitDeposition(m_fields, m_laser, geom[lev], lev); + m_multi_plasma.ExplicitDeposition(m_fields, m_multi_laser, geom[lev], lev); // Solves Bx, By using Sx, Sy and chi ExplicitMGSolveBxBy(lev, WhichSlice::This); @@ -664,7 +664,7 @@ Hipace::ExplicitSolveOneSubSlice (const int lev, const int step, const int 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_laser, geom[lev], + m_multi_plasma.AdvanceParticles(m_fields, m_multi_laser, geom[lev], false, true, true, !do_salame, lev); // Push beam particles @@ -685,13 +685,13 @@ Hipace::PredictorCorrectorSolveOneSubSlice (const int lev, const int step, const } // push plasma - m_multi_plasma.AdvanceParticles(m_fields, m_laser, geom[lev], false, + 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_laser, WhichSlice::This, false, true, true, true, m_use_laser, geom[lev], lev); + m_fields, m_multi_laser, WhichSlice::This, false, true, true, true, m_use_laser, geom[lev], lev); m_fields.AddRhoIons(lev); @@ -747,7 +747,7 @@ Hipace::ResetLaser () HIPACE_PROFILE("Hipace::ResetLaser()"); for (int sl=WhichLaserSlice::nm1j00; slalloc(buffer_size); @@ -1198,9 +1198,9 @@ Hipace::Wait (const int step, int it, bool only_ghost) // Receive laser { if (only_ghost) return; - if (!m_laser.m_use_laser) return; + if (!m_multi_laser.m_use_laser) return; AMREX_ALWAYS_ASSERT(nz_laser > 0); - amrex::FArrayBox& laser_fab = m_laser.getFAB(); + amrex::FArrayBox& laser_fab = m_multi_laser.getFAB(); amrex::Array4 laser_arr = laser_fab.array(); const amrex::Box& bx = laser_fab.box(); // does not include ghost cells AMREX_ALWAYS_ASSERT_WITH_MESSAGE( @@ -1209,7 +1209,7 @@ Hipace::Wait (const int step, int it, bool only_ghost) const std::size_t nreals = bx.numPts()*laser_fab.nComp(); MPI_Status lstatus; - if (m_laser.is3dOnHost()) { + if (m_multi_laser.is3dOnHost()) { // Directly receive envelope in laser fab MPI_Recv(laser_fab.dataPtr(), nreals, amrex::ParallelDescriptor::Mpi_typemap::type(), @@ -1249,8 +1249,8 @@ Hipace::Notify (const int step, const int it, const bool only_time = m_physical_time >= m_max_time; NotifyFinish(it, only_ghost, only_time); // finish the previous send int nz_laser = 0; - if (m_laser.m_use_laser){ - const amrex::Box& laser_bx = m_laser.getFAB().box(); + if (m_multi_laser.m_use_laser){ + const amrex::Box& laser_bx = m_multi_laser.getFAB().box(); nz_laser = laser_bx.bigEnd(2) - laser_bx.smallEnd(2) + 1; } const int nbeams = m_multi_beam.get_nbeams(); @@ -1309,7 +1309,7 @@ Hipace::Notify (const int step, const int it, // Send beam particles. Currently only one tile. { const amrex::Long np_total = std::accumulate(np_snd.begin(), np_snd.begin()+nbeams, 0); - if (np_total == 0 && !m_laser.m_use_laser) return; + if (np_total == 0 && !m_multi_laser.m_use_laser) return; const amrex::Long psize = sizeof(BeamParticleContainer::SuperParticleType); const amrex::Long buffer_size = psize*np_total; char*& psend_buffer = only_ghost ? m_psend_buffer_ghost : m_psend_buffer; @@ -1402,14 +1402,14 @@ Hipace::Notify (const int step, const int it, // Send laser data { if (only_ghost) return; - if (!m_laser.m_use_laser) return; - const amrex::FArrayBox& laser_fab = m_laser.getFAB(); + if (!m_multi_laser.m_use_laser) return; + const amrex::FArrayBox& laser_fab = m_multi_laser.getFAB(); amrex::Array4 const& laser_arr = laser_fab.array(); const amrex::Box& lbx = laser_fab.box(); // does not include ghost cells const std::size_t nreals = lbx.numPts()*laser_fab.nComp(); m_lsend_buffer = (amrex::Real*)amrex::The_Pinned_Arena()->alloc (sizeof(amrex::Real)*nreals); - if (m_laser.is3dOnHost()) { + if (m_multi_laser.is3dOnHost()) { amrex::Gpu::streamSynchronize(); // Copy from laser envelope 3D array (on host) to MPI buffer (on host) laser_fab.copyToMem(lbx, 0, laser_fab.nComp(), m_lsend_buffer); @@ -1521,7 +1521,7 @@ Hipace::FillDiagnostics (const int lev, int i_slice) if (m_diags.hasField()[lev]) { m_fields.Copy(lev, i_slice, m_diags.getGeom()[lev], m_diags.getF(lev), m_diags.getF(lev).box(), Geom(lev), - m_diags.getCompsIdx(), m_diags.getNFields(), m_laser); + m_diags.getCompsIdx(), m_diags.getNFields(), m_multi_laser); } } diff --git a/src/fields/Fields.H b/src/fields/Fields.H index d5cfda3040..45cffd82a7 100644 --- a/src/fields/Fields.H +++ b/src/fields/Fields.H @@ -11,7 +11,7 @@ #include "fft_poisson_solver/FFTPoissonSolver.H" #include "diagnostics/Diagnostic.H" -#include "laser/Laser.H" +#include "laser/MultiLaser.H" #include "utils/GPUUtil.H" #include @@ -21,7 +21,7 @@ #include class Hipace; -class Laser; +class MultiLaser; /** \brief describes which slice with respect to the currently calculated is used */ struct WhichSlice { @@ -158,12 +158,12 @@ public: * \param[in] calc_geom main geometry * \param[in] diag_comps_vect the field components to copy * \param[in] ncomp number of components to copy - * \param[in] laser Laser object + * \param[in] multi_laser MultiLaser object */ void Copy (const int lev, const int i_slice, const amrex::Geometry& diag_geom, amrex::FArrayBox& diag_fab, amrex::Box diag_box, const amrex::Geometry& calc_geom, const amrex::Gpu::DeviceVector& diag_comps_vect, const int ncomp, - Laser& laser); + MultiLaser& multi_laser); /** \brief Shift slices by 1 element: slices (1,2) are then stored in (2,3). * diff --git a/src/fields/Fields.cpp b/src/fields/Fields.cpp index 0c76401aa9..1e1e9d4a20 100644 --- a/src/fields/Fields.cpp +++ b/src/fields/Fields.cpp @@ -406,7 +406,7 @@ void Fields::Copy (const int lev, const int i_slice, const amrex::Geometry& diag_geom, amrex::FArrayBox& diag_fab, amrex::Box diag_box, const amrex::Geometry& calc_geom, const amrex::Gpu::DeviceVector& diag_comps_vect, const int ncomp, - Laser& laser) + MultiLaser& multi_laser) { HIPACE_PROFILE("Fields::Copy()"); constexpr int depos_order_xy = 1; @@ -461,7 +461,7 @@ Fields::Copy (const int lev, const int i_slice, const amrex::Geometry& diag_geom if (diag_box.isEmpty()) return; auto& slice_mf = m_slices[lev][WhichSlice::This]; auto slice_func = interpolated_field_xy{{slice_mf}, calc_geom}; - auto& laser_mf = laser.getSlices(WhichLaserSlice::n00j00); + auto& laser_mf = multi_laser.getSlices(WhichLaserSlice::n00j00); auto laser_func = interpolated_field_xy{{laser_mf}, calc_geom}; #ifdef AMREX_USE_GPU @@ -493,7 +493,7 @@ Fields::Copy (const int lev, const int i_slice, const amrex::Geometry& diag_geom diag_array(i,j,k,n) += rel_z_data[k-k_min] * slice_array(x,y,m); }); - if (!laser.m_use_laser) continue; + if (!multi_laser.m_use_laser) continue; auto laser_array = laser_func.array(mfi); amrex::ParallelFor(diag_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept diff --git a/src/laser/CMakeLists.txt b/src/laser/CMakeLists.txt index 2e94fedfc7..99a73ed011 100644 --- a/src/laser/CMakeLists.txt +++ b/src/laser/CMakeLists.txt @@ -1,4 +1,5 @@ target_sources(HiPACE PRIVATE + MultiLaser.cpp Laser.cpp ) diff --git a/src/laser/Laser.H b/src/laser/Laser.H index c99678bc1f..8d51cf2d31 100644 --- a/src/laser/Laser.H +++ b/src/laser/Laser.H @@ -1,226 +1,43 @@ +/* Copyright 2022 + * + * This file is part of HiPACE++. + * + * Authors: MaxThevenet, AlexanderSinn + * Severin Diederichs, atmyers, Angel Ferran Pousa + * License: BSD-3-Clause-LBNL + */ + #ifndef LASER_H_ #define LASER_H_ -#include "fields/Fields.H" -#include "mg_solver/HpMultiGrid.H" - -#include #include -#include -#include - -#ifdef AMREX_USE_CUDA -# include -#elif defined(AMREX_USE_HIP) -# include -# include -#else -# include -#endif - -namespace LaserFFT { -#ifdef AMREX_USE_CUDA - using VendorFFT = cufftHandle; - const auto VendorCreate = cufftPlan2d; - const auto VendorDestroy = cufftDestroy; -# ifdef AMREX_USE_FLOAT - const auto VendorExecute = cufftExecC2C; - const auto cufft_type = CUFFT_C2C; - using cufftComplex = cuComplex; -# else - const auto VendorExecute = cufftExecZ2Z; - const auto cufft_type = CUFFT_Z2Z; - using cufftComplex = cuDoubleComplex; -# endif -#elif defined(AMREX_USE_HIP) - using VendorFFT = rocfft_plan; - const auto VendorDestroy = rocfft_plan_destroy; -// TODO -#else -# ifdef AMREX_USE_FLOAT - using VendorFFT = fftwf_plan; - using FFTWComplex = fftwf_complex; - const auto VendorCreate = fftwf_plan_dft_2d; - const auto VendorExecute = fftwf_execute; - const auto VendorDestroy = fftwf_destroy_plan; -# else - using VendorFFT = fftw_plan; - using FFTWComplex = fftw_complex; - const auto VendorCreate = fftw_plan_dft_2d; - const auto VendorExecute = fftw_execute; - const auto VendorDestroy = fftw_destroy_plan; -# endif -#endif -} - -/** \brief describes which slice with respect to the currently calculated is used */ -struct WhichLaserSlice { - // n stands for the time step, j for the longitudinal slice. - // n00 is time step n, nm1 is n-1 and np1 is n+1. Similar notation for slice j. - enum slice { nm1j00, nm1jp1, nm1jp2, n00j00, n00jp1, n00jp2, np1j00, np1jp1, np1jp2, N }; -}; - -class Fields; +#include class Laser { - - using SpectralFieldLoc = amrex::BaseFab >; - -private: - - static constexpr int m_nslices = WhichLaserSlice::N; - public: - - /** Constructor */ - explicit Laser () - { - ReadParameters(); - } - - ~Laser () - { - if (!m_use_laser) return; - if (m_solver_type == "fft") { - LaserFFT::VendorDestroy( m_plan_fwd ); - LaserFFT::VendorDestroy( m_plan_bkw ); - } - } - - void ReadParameters (); - - /** get function for the 2D slices - * \param[in] islice slice index - */ - amrex::MultiFab& getSlices (int islice) {return m_slices[islice]; } - - /** get function for the 2D slices (const version) - * \param[in] islice slice index - */ - const amrex::MultiFab& getSlices (int islice) const {return m_slices[islice]; } - - /** Return the 3D FArrayBox containing the laser envelope, const version */ - const amrex::FArrayBox& getFAB () const {return m_F;} - - /** Retun the 3D FArrayBox containing the laser envelope, non-const version */ - amrex::FArrayBox& getFAB () {return m_F;} - - /** getter function, whether 3D envelope is store in host memory */ - int is3dOnHost () const {return m_3d_on_host;} - - /** \brief Allocate laser multifab - * \param[in] slice_ba box array of the slice - * \param[in] slice_dm corresponding distribution mapping - */ - void InitData (const amrex::BoxArray& slice_ba, - const amrex::DistributionMapping& slice_dm); - - /** \brief Initialize 3D laser field on current box. - * - * \param[in] step time step of the simulation - * \param[in] bx Box on which the laser field is initialized - * \param[in] gm Geometry of the problem - */ - void Init3DEnvelope (int step, amrex::Box bx, const amrex::Geometry& gm); - - /** \brief Copy from 2D slice on device to 3D array on host, and vice-versa - * - * \param[in] isl slice index, referring to the 3D slice - * \param[in] to3d if true, copy from 2D slice to 3D array. Otherwise, the other way. - */ - void Copy (int isl, bool to3d); - - /** Wrapper function to advance a laser slice by 1 time step. - * \param[in] fields Field object - * \param[in] geom Geometry object - * \param[in] dt time step of the simulation - * \param[in] step current iteration. Needed because step 0 needs a specific treatment. - */ - void AdvanceSlice (const Fields& fields, const amrex::Geometry& geom, amrex::Real dt, int step); - - /** Advance a laser slice by 1 time step using a multigrid solver. - * The complex phase of the envelope is evaluated on-axis only, but can be generalized to everywhere. - * - * \param[in] fields Field object - * \param[in] geom Geometry object - * \param[in] dt time step of the simulation - * \param[in] step current iteration. Needed because step 0 needs a specific treatment. - */ - void AdvanceSliceMG (const Fields& fields, const amrex::Geometry& geom, amrex::Real dt, int step); - - /** Advance a laser slice by 1 time step using a FFT solver. - * The complex phase of the envelope is evaluated on-axis only. - * - * \param[in] fields Field object - * \param[in] geom Geometry object - * \param[in] dt time step of the simulation - * \param[in] step current iteration. Needed because step 0 needs a specific treatment. - */ - void AdvanceSliceFFT (const Fields& fields, const amrex::Geometry& geom, amrex::Real dt, int step); - - /** Initialize 1 longitudinal slice of the laser, and store it in n00j00 (current time step) - * and nm1j00 (previous time step). - * - * \param[in] geom Geometry object for the slice - * \param[in] islice slice index - */ - void InitLaserSlice (const amrex::Geometry& geom, const int islice); - - bool m_use_laser {false}; /**< whether a laser is used or not */ - -private: - - std::string m_name = "laser"; /**< name of the laser */ + Laser (std::string name); +/* + amrex::Real getA0 () {return m_a0}; + amrex::Real getK0 () {return m_k0}; + amrex::Real getW0 () {return m_w0}; + amrex::Real getL0 () {return m_l0}; + amrex::Real getTau () {return m_tau}; + amrex::Real getLambda () {return m_lambda0}; + amrex::Real getFoc () {return m_focal_distance}; + amrex::Real getX0 () {return m_position_mean[0]}; + amrex::Real getY0 () {return m_position_mean[1]}; + amrex::Real getZ0 () {return m_position_mean[2]}; +*/ + std::string m_name {""}; amrex::Real m_a0 {0.}; /**< Laser peak normalized amplitude */ - amrex::RealVect m_w0 {0., 0., 0.}; /**< Laser waist in x and y (the third value is omitted) */ + amrex::Real m_w0 {0.}; /**< Laser waist */ amrex::Real m_L0 {0.}; /**< Laser length (HW 1/e in amplitude) */ amrex::Real m_tau {0.}; /**< Laser duration (HW 1/e in amplitude) */ - amrex::Real m_lambda0 {0.}; /**< Laser central wavelength */ + /** Focal distance of the laser pulse */ + amrex::Real m_focal_distance {0.}; /** Average position of the Gaussian laser pulse */ amrex::RealVect m_position_mean {0., 0., 0.}; - int m_3d_on_host {0};/** Whether the 3D laser envelope is stored in host or device memory */ - amrex::Real m_focal_distance {0.}; - /** Number of guard cells for slices MultiFab */ - amrex::IntVect m_slices_nguards = {-1, -1, -1}; - std::string m_solver_type = "multigrid"; - bool m_use_phase {true}; - amrex::Box m_slice_box; - - /** Nb fields in 3D array: new_real, new_imag, old_real, old_imag */ - int m_nfields_3d {4}; - /** 3D laser data. Could be a vector over levels if MR needed for laser */ - amrex::FArrayBox m_F; - /** Array of N slices required to compute current slice */ - std::array m_slices; - amrex::Real m_MG_tolerance_rel = 1.e-4; - amrex::Real m_MG_tolerance_abs = 0.; - int m_MG_verbose = 0; - /** Whether to use time-averaged RHS in envelope solver. */ - bool m_MG_average_rhs = true; - /** hpmg solver for the envelope solver */ - std::unique_ptr m_mg; - - // Elements for the FFT-based laser envelope solver - // This could belong to AnyFFT etc., with 2 caveats: - // - This solver operates on a FArrayBox instead of a MultiFab, which is more adequate - // - The array in position space must be Complex rather than real, which takes up quite some - // rewriting, see https://github.com/MaxThevenet/hipace/tree/laser_solve, - // not sure what the best way to proceed is. - /** FFTW plan for forward C2C transform to solve Complex Poisson equation */ - LaserFFT::VendorFFT m_plan_fwd; - /** FFTW plan for backward C2C transform to solve Complex Poisson equation */ - LaserFFT::VendorFFT m_plan_bkw; - /** Complex FAB to store the solution (e.g. laser envelope on current slice) */ - SpectralFieldLoc m_sol; - /** Complex FAB to store the RHS in position space */ - SpectralFieldLoc m_rhs; - /** Complex FAB to store the RHS in Fourier space */ - SpectralFieldLoc m_rhs_fourier; -#ifdef AMREX_USE_CUDA - cufftResult m_result_fwd; - cufftResult m_result_bkw; -#endif }; #endif // LASER_H_ diff --git a/src/laser/Laser.cpp b/src/laser/Laser.cpp index dfbe4b2695..19041f968c 100644 --- a/src/laser/Laser.cpp +++ b/src/laser/Laser.cpp @@ -1,759 +1,35 @@ +/* Copyright 2022 + * + * This file is part of HiPACE++. + * + * Authors: MaxThevenet, AlexanderSinn + * Severin Diederichs, atmyers, Angel Ferran Pousa + * License: BSD-3-Clause-LBNL + */ + #include "Laser.H" -#include "utils/Constants.H" +#include "utils/Parser.H" #include "Hipace.H" -#include "utils/HipaceProfilerWrapper.H" -#ifdef AMREX_USE_CUDA -# include "fields/fft_poisson_solver/fft/CuFFTUtils.H" -#endif -#include -#ifdef AMREX_USE_CUDA -# include -#elif defined(AMREX_USE_HIP) -# include -# include -#else -# include -#endif +#include +#include -void -Laser::ReadParameters () +Laser::Laser (std::string name) { + m_name = name; amrex::ParmParse pp(m_name); - queryWithParser(pp, "use_laser", m_use_laser); queryWithParser(pp, "a0", m_a0); - if (!m_use_laser) return; -#if defined(AMREX_USE_HIP) - amrex::Abort("Laser solver not implemented with HIP"); -#endif - amrex::Vector tmp_vector; - if (queryWithParser(pp, "w0", tmp_vector)){ - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(tmp_vector.size() == 2, - "The laser waist w0 must be provided in x and y, " - "so laser.w0 should contain 2 values"); - for (int i=0; i<2; i++) m_w0[i] = tmp_vector[i]; - } - + queryWithParser(pp, "w0", m_w0); bool length_is_specified = queryWithParser(pp, "L0", m_L0);; bool duration_is_specified = queryWithParser(pp, "tau", m_tau); AMREX_ALWAYS_ASSERT_WITH_MESSAGE( length_is_specified + duration_is_specified == 1, "Please specify exlusively either the pulse length L0 or the duration tau of the laser"); if (duration_is_specified) m_L0 = m_tau/get_phys_const().c; - getWithParser(pp, "lambda0", m_lambda0); - amrex::Array loc_array; queryWithParser(pp, "focal_distance", m_focal_distance); - queryWithParser(pp, "position_mean", loc_array); - queryWithParser(pp, "3d_on_host", m_3d_on_host); - queryWithParser(pp, "use_phase", m_use_phase); - queryWithParser(pp, "solver_type", m_solver_type); - AMREX_ALWAYS_ASSERT(m_solver_type == "multigrid" || m_solver_type == "fft"); - - bool mg_param_given = queryWithParser(pp, "MG_tolerance_rel", m_MG_tolerance_rel); - mg_param_given += queryWithParser(pp, "MG_tolerance_abs", m_MG_tolerance_abs); - mg_param_given += queryWithParser(pp, "MG_verbose", m_MG_verbose); - mg_param_given += queryWithParser(pp, "MG_average_rhs", m_MG_average_rhs); - - // Raise warning if user specifies MG parameters without using the MG solver - if (mg_param_given && (m_solver_type != "multigrid")) { - amrex::Print()<<"WARNING: parameters laser.MG_... only active if laser.solver_type = multigrid\n"; - } + amrex::Array loc_array; + queryWithParser(pp, "position_mean", loc_array); for (int idim=0; idim < AMREX_SPACEDIM; ++idim) m_position_mean[idim] = loc_array[idim]; -} - - -void -Laser::InitData (const amrex::BoxArray& slice_ba, - const amrex::DistributionMapping& slice_dm) -{ - if (!m_use_laser) return; - - HIPACE_PROFILE("Laser::InitData()"); - - // Alloc 2D slices - // 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}; - AMREX_ALWAYS_ASSERT(WhichLaserSlice::N == m_nslices); - for (int islice=0; islice(m_rhs.dataPtr()), - reinterpret_cast(m_rhs_fourier.dataPtr()), - FFTW_FORWARD, FFTW_ESTIMATE); - // Backward FFT plan - m_plan_bkw = LaserFFT::VendorCreate( - fft_size[1], fft_size[0], - reinterpret_cast(m_rhs_fourier.dataPtr()), - reinterpret_cast(m_sol.dataPtr()), - FFTW_BACKWARD, FFTW_ESTIMATE); -#endif - } -} - -void -Laser::Init3DEnvelope (int step, amrex::Box bx, const amrex::Geometry& gm) -{ - - if (!m_use_laser) return; - - HIPACE_PROFILE("Laser::Init3DEnvelope()"); - // Allocate the 3D field on this box - // Note: box has no guard cells - m_F.resize(bx, m_nfields_3d, m_3d_on_host ? amrex::The_Pinned_Arena() : amrex::The_Arena()); - - if (step > 0) return; - - // In order to use the normal Copy function, we use slice np1j00 as a tmp array - // to initialize the laser in the loop over slices below. - // We need to keep the value in np1j00, as it is shifted to np1jp1 and used to compute - // the following slices. This is relevant for the first slices at step 0 of every box - // (except for the head box). - amrex::FArrayBox store_np1j00; - store_np1j00.resize(m_slice_box, 2, amrex::The_Arena()); - store_np1j00.copy(m_slices[WhichLaserSlice::np1j00][0]); - - // Loop over slices - for (int isl = bx.bigEnd(Direction::z); isl >= bx.smallEnd(Direction::z); --isl){ - // Compute initial field on the current (device) slice np1j00 - InitLaserSlice(gm, isl); - // Copy: (device) slice np1j00 to the right location - // in the (host) 3D array of the current time. - Copy(isl, true); - } - - // Reset np1j00 to its original value. - m_slices[WhichLaserSlice::np1j00][0].copy(store_np1j00); -} - -void -Laser::Copy (int isl, bool to3d) -{ - if (!m_use_laser) return; - - using namespace amrex::literals; - - HIPACE_PROFILE("Laser::Copy()"); - - amrex::MultiFab& nm1j00 = m_slices[WhichLaserSlice::nm1j00]; - amrex::MultiFab& nm1jp1 = m_slices[WhichLaserSlice::nm1jp1]; - amrex::MultiFab& nm1jp2 = m_slices[WhichLaserSlice::nm1jp2]; - amrex::MultiFab& n00j00 = m_slices[WhichLaserSlice::n00j00]; - amrex::MultiFab& n00jp1 = m_slices[WhichLaserSlice::n00jp1]; - amrex::MultiFab& n00jp2 = m_slices[WhichLaserSlice::n00jp2]; - amrex::MultiFab& np1j00 = m_slices[WhichLaserSlice::np1j00]; - amrex::MultiFab& np1jp1 = m_slices[WhichLaserSlice::np1jp1]; - amrex::MultiFab& np1jp2 = m_slices[WhichLaserSlice::np1jp2]; - - for ( amrex::MFIter mfi(n00j00, DfltMfi); mfi.isValid(); ++mfi ){ - const amrex::Box& bx = mfi.tilebox(); - Array3 nm1j00_arr = nm1j00.array(mfi); - Array3 nm1jp1_arr = nm1jp1.array(mfi); - Array3 nm1jp2_arr = nm1jp2.array(mfi); - Array3 n00j00_arr = n00j00.array(mfi); - Array3 n00jp1_arr = n00jp1.array(mfi); - Array3 n00jp2_arr = n00jp2.array(mfi); - Array3 np1j00_arr = np1j00.array(mfi); - Array3 np1jp1_arr = np1jp1.array(mfi); - Array3 np1jp2_arr = np1jp2.array(mfi); - amrex::Array4 host_arr = m_F.array(); - amrex::ParallelFor( - bx, 2, - [=] AMREX_GPU_DEVICE(int i, int j, int, int n) noexcept - { - // +2 in 3D array means old. - // 2 components for complex numbers. - if (to3d){ - // this slice into old host - host_arr(i,j,isl,n+2) = n00j00_arr(i,j,n); - // next time slice into new host - host_arr(i,j,isl,n ) = np1j00_arr(i,j,n); - } else { - // Shift slices of step n-1, and get current slice from 3D array - nm1jp2_arr(i,j,n) = nm1jp1_arr(i,j,n); - nm1jp1_arr(i,j,n) = nm1j00_arr(i,j,n); - nm1j00_arr(i,j,n) = host_arr(i,j,isl,n+2); - // Shift slices of step n, and get current slice from 3D array - n00jp2_arr(i,j,n) = n00jp1_arr(i,j,n); - n00jp1_arr(i,j,n) = n00j00_arr(i,j,n); - n00j00_arr(i,j,n) = host_arr(i,j,isl,n); - // Shift slices of step n+1. Current slice will be computed - np1jp2_arr(i,j,n) = np1jp1_arr(i,j,n); - np1jp1_arr(i,j,n) = np1j00_arr(i,j,n); - } - }); - } -} - -void -Laser::AdvanceSlice (const Fields& fields, const amrex::Geometry& geom, amrex::Real dt, int step) -{ - - if (!m_use_laser) return; - - if (m_solver_type == "multigrid") { - AdvanceSliceMG(fields, geom, dt, step); - } else if (m_solver_type == "fft") { - AdvanceSliceFFT(fields, geom, dt, step); - } else { - amrex::Abort("laser.solver_type must be fft or multigrid"); - } -} - -void -Laser::AdvanceSliceMG (const Fields& fields, const amrex::Geometry& geom, amrex::Real dt, int step) -{ - - HIPACE_PROFILE("Laser::AdvanceSliceMG()"); - - using namespace amrex::literals; - using Complex = amrex::GpuComplex; - constexpr Complex I(0.,1.); - - const amrex::Real dx = geom.CellSize(0); - const amrex::Real dy = geom.CellSize(1); - const amrex::Real dz = geom.CellSize(2); - - const PhysConst phc = get_phys_const(); - const amrex::Real c = phc.c; - const amrex::Real k0 = 2.*MathConst::pi/m_lambda0; - const bool do_avg_rhs = m_MG_average_rhs; - - amrex::MultiFab& nm1j00 = m_slices[WhichLaserSlice::nm1j00]; - amrex::MultiFab& nm1jp1 = m_slices[WhichLaserSlice::nm1jp1]; - amrex::MultiFab& nm1jp2 = m_slices[WhichLaserSlice::nm1jp2]; - amrex::MultiFab& n00j00 = m_slices[WhichLaserSlice::n00j00]; - amrex::MultiFab& n00jp1 = m_slices[WhichLaserSlice::n00jp1]; - amrex::MultiFab& n00jp2 = m_slices[WhichLaserSlice::n00jp2]; - amrex::MultiFab& np1j00 = m_slices[WhichLaserSlice::np1j00]; - amrex::MultiFab& np1jp1 = m_slices[WhichLaserSlice::np1jp1]; - amrex::MultiFab& np1jp2 = m_slices[WhichLaserSlice::np1jp2]; - - amrex::FArrayBox rhs_mg; - amrex::FArrayBox acoeff_real; - amrex::Real acoeff_real_scalar = 0._rt; - amrex::Real acoeff_imag_scalar = 0._rt; - - amrex::Real djn {0.}; - - for ( amrex::MFIter mfi(n00j00, DfltMfi); mfi.isValid(); ++mfi ){ - const amrex::Box& bx = mfi.tilebox(); - const int imin = bx.smallEnd(0); - const int imax = bx.bigEnd (0); - const int jmin = bx.smallEnd(1); - const int jmax = bx.bigEnd (1); - - acoeff_real.resize(bx, 1, amrex::The_Arena()); - rhs_mg.resize(bx, 2, amrex::The_Arena()); - Array3 nm1j00_arr = nm1j00.array(mfi); - Array3 nm1jp1_arr = nm1jp1.array(mfi); - Array3 nm1jp2_arr = nm1jp2.array(mfi); - Array3 n00j00_arr = n00j00.array(mfi); - Array3 n00jp1_arr = n00jp1.array(mfi); - Array3 n00jp2_arr = n00jp2.array(mfi); - Array3 np1jp1_arr = np1jp1.array(mfi); - Array3 np1jp2_arr = np1jp2.array(mfi); - Array3 rhs_mg_arr = rhs_mg.array(); - Array3 acoeff_real_arr = acoeff_real.array(); - Array3 rhs_arr = m_rhs.array(); - - constexpr int lev = 0; - const amrex::FArrayBox& isl_fab = fields.getSlices(lev, WhichSlice::This)[mfi]; - Array3 const isl_arr = isl_fab.array(); - const int chi = Comps[WhichSlice::This]["chi"]; - - // Calculate phase terms. 0 if !m_use_phase - amrex::Real tj00 = 0.; - amrex::Real tjp1 = 0.; - amrex::Real tjp2 = 0.; - - if (m_use_phase) { - int const Nx = bx.length(0); - int const Ny = bx.length(1); - - // Get the central point. - int const imid = (Nx+1)/2; - int const jmid = (Ny+1)/2; - - // Calculate complex arguments (theta) needed - // Just once, on axis, as done in Wake-T - // This is done with a reduce operation, returning the sum of the four elements nearest - // the axis (both real and imag parts, and for the 3 arrays relevant) ... - amrex::ReduceOps< - amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum, - amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum> reduce_op; - amrex::ReduceData< - amrex::Real, amrex::Real, amrex::Real, - amrex::Real, amrex::Real, amrex::Real> reduce_data(reduce_op); - using ReduceTuple = typename decltype(reduce_data)::Type; - reduce_op.eval(bx, reduce_data, - [=] AMREX_GPU_DEVICE (int i, int j, int) -> ReduceTuple - { - if ( ( i == imid-1 || i == imid ) && ( j == jmid-1 || j == jmid ) ) { - return { - n00j00_arr(i,j,0), n00j00_arr(i,j,1), - n00jp1_arr(i,j,0), n00jp1_arr(i,j,1), - n00jp2_arr(i,j,0), n00jp2_arr(i,j,1) - }; - } else { - return {0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt}; - } - }); - // ... and taking the argument of the resulting complex number. - ReduceTuple hv = reduce_data.value(reduce_op); - tj00 = std::atan2(amrex::get<1>(hv), amrex::get<0>(hv)); - tjp1 = std::atan2(amrex::get<3>(hv), amrex::get<2>(hv)); - tjp2 = std::atan2(amrex::get<5>(hv), amrex::get<4>(hv)); - } - - amrex::Real dt1 = tj00 - tjp1; - amrex::Real dt2 = tjp1 - tjp2; - if (dt1 <-1.5_rt*MathConst::pi) dt1 += 2._rt*MathConst::pi; - if (dt1 > 1.5_rt*MathConst::pi) dt1 -= 2._rt*MathConst::pi; - if (dt2 <-1.5_rt*MathConst::pi) dt2 += 2._rt*MathConst::pi; - if (dt2 > 1.5_rt*MathConst::pi) dt2 -= 2._rt*MathConst::pi; - Complex exp1 = amrex::exp(I*(tj00-tjp1)); - Complex exp2 = amrex::exp(I*(tj00-tjp2)); - - // D_j^n as defined in Benedetti's 2017 paper - djn = ( -3._rt*dt1 + dt2 ) / (2._rt*dz); - acoeff_real_scalar = step == 0 ? 6._rt/(c*dt*dz) - : 3._rt/(c*dt*dz) + 2._rt/(c*c*dt*dt); - acoeff_imag_scalar = step == 0 ? -4._rt * ( k0 + djn ) / (c*dt) - : -2._rt * ( k0 + djn ) / (c*dt); - - amrex::ParallelFor( - bx, 1, - [=] AMREX_GPU_DEVICE(int i, int j, int, int) noexcept - { - // Transverse Laplacian of real and imaginary parts of A_j^n-1 - amrex::Real lapR, lapI; - if (step == 0) { - lapR = i>imin && ijmin && jimin && ijmin && jimin && ijmin && jimin && ijmin && j(slice_geom); - } - - const int max_iters = 200; - m_mg->solve2(np1j00[0], rhs_mg, acoeff_real, acoeff_imag_scalar, - m_MG_tolerance_rel, m_MG_tolerance_abs, max_iters, m_MG_verbose); -} - -void -Laser::AdvanceSliceFFT (const Fields& fields, const amrex::Geometry& geom, const amrex::Real dt, int step) -{ - - HIPACE_PROFILE("Laser::AdvanceSliceFFT()"); - - using namespace amrex::literals; - using Complex = amrex::GpuComplex; - constexpr Complex I(0.,1.); - - const amrex::Real dx = geom.CellSize(0); - const amrex::Real dy = geom.CellSize(1); - const amrex::Real dz = geom.CellSize(2); - - const PhysConst phc = get_phys_const(); - const amrex::Real c = phc.c; - const amrex::Real k0 = 2.*MathConst::pi/m_lambda0; - - amrex::MultiFab& nm1j00 = m_slices[WhichLaserSlice::nm1j00]; - amrex::MultiFab& nm1jp1 = m_slices[WhichLaserSlice::nm1jp1]; - amrex::MultiFab& nm1jp2 = m_slices[WhichLaserSlice::nm1jp2]; - amrex::MultiFab& n00j00 = m_slices[WhichLaserSlice::n00j00]; - amrex::MultiFab& n00jp1 = m_slices[WhichLaserSlice::n00jp1]; - amrex::MultiFab& n00jp2 = m_slices[WhichLaserSlice::n00jp2]; - amrex::MultiFab& np1j00 = m_slices[WhichLaserSlice::np1j00]; - amrex::MultiFab& np1jp1 = m_slices[WhichLaserSlice::np1jp1]; - amrex::MultiFab& np1jp2 = m_slices[WhichLaserSlice::np1jp2]; - - for ( amrex::MFIter mfi(n00j00, DfltMfi); mfi.isValid(); ++mfi ){ - const amrex::Box& bx = mfi.tilebox(); - const int imin = bx.smallEnd(0); - const int imax = bx.bigEnd (0); - const int jmin = bx.smallEnd(1); - const int jmax = bx.bigEnd (1); - - // solution: complex array - // The right-hand side is computed and stored in rhs - // Then rhs is Fourier-transformed into rhs_fourier, then multiplied by -1/(k**2+a) - // rhs_fourier is FFT-back-transformed to sol, and sol is normalized and copied into np1j00. - Array3 sol_arr = m_sol.array(); - Array3 rhs_arr = m_rhs.array(); - amrex::Array4 rhs_fourier_arr = m_rhs_fourier.array(); - - Array3 nm1j00_arr = nm1j00.array(mfi); - Array3 nm1jp1_arr = nm1jp1.array(mfi); - Array3 nm1jp2_arr = nm1jp2.array(mfi); - Array3 n00j00_arr = n00j00.array(mfi); - Array3 n00jp1_arr = n00jp1.array(mfi); - Array3 n00jp2_arr = n00jp2.array(mfi); - Array3 np1j00_arr = np1j00.array(mfi); - Array3 np1jp1_arr = np1jp1.array(mfi); - Array3 np1jp2_arr = np1jp2.array(mfi); - - constexpr int lev = 0; - const amrex::FArrayBox& isl_fab = fields.getSlices(lev, WhichSlice::This)[mfi]; - Array3 const isl_arr = isl_fab.array(); - const int chi = Comps[WhichSlice::This]["chi"]; - - int const Nx = bx.length(0); - int const Ny = bx.length(1); - - // Get the central point. Useful to get the on-axis phase and calculate kx and ky. - int const imid = (Nx+1)/2; - int const jmid = (Ny+1)/2; - - // Calculate phase terms. 0 if !m_use_phase - amrex::Real tj00 = 0.; - amrex::Real tjp1 = 0.; - amrex::Real tjp2 = 0.; - - if (m_use_phase) { - // Calculate complex arguments (theta) needed - // Just once, on axis, as done in Wake-T - // This is done with a reduce operation, returning the sum of the four elements nearest - // the axis (both real and imag parts, and for the 3 arrays relevant) ... - amrex::ReduceOps< - amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum, - amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum> reduce_op; - amrex::ReduceData< - amrex::Real, amrex::Real, amrex::Real, - amrex::Real, amrex::Real, amrex::Real> reduce_data(reduce_op); - using ReduceTuple = typename decltype(reduce_data)::Type; - reduce_op.eval(bx, reduce_data, - [=] AMREX_GPU_DEVICE (int i, int j, int) -> ReduceTuple - { - if ( ( i == imid-1 || i == imid ) && ( j == jmid-1 || j == jmid ) ) { - return { - n00j00_arr(i,j,0), n00j00_arr(i,j,1), - n00jp1_arr(i,j,0), n00jp1_arr(i,j,1), - n00jp2_arr(i,j,0), n00jp2_arr(i,j,1) - }; - } else { - return {0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt}; - } - }); - // ... and taking the argument of the resulting complex number. - ReduceTuple hv = reduce_data.value(reduce_op); - tj00 = std::atan2(amrex::get<1>(hv), amrex::get<0>(hv)); - tjp1 = std::atan2(amrex::get<3>(hv), amrex::get<2>(hv)); - tjp2 = std::atan2(amrex::get<5>(hv), amrex::get<4>(hv)); - } - - amrex::Real dt1 = tj00 - tjp1; - amrex::Real dt2 = tjp1 - tjp2; - if (dt1 <-1.5_rt*MathConst::pi) dt1 += 2._rt*MathConst::pi; - if (dt1 > 1.5_rt*MathConst::pi) dt1 -= 2._rt*MathConst::pi; - if (dt2 <-1.5_rt*MathConst::pi) dt2 += 2._rt*MathConst::pi; - if (dt2 > 1.5_rt*MathConst::pi) dt2 -= 2._rt*MathConst::pi; - Complex exp1 = amrex::exp(I*(tj00-tjp1)); - Complex exp2 = amrex::exp(I*(tj00-tjp2)); - - // D_j^n as defined in Benedetti's 2017 paper - amrex::Real djn = ( -3._rt*dt1 + dt2 ) / (2._rt*dz); - amrex::ParallelFor( - bx, 1, - [=] AMREX_GPU_DEVICE(int i, int j, int, int) noexcept - { - // Transverse Laplacian of real and imaginary parts of A_j^n-1 - amrex::Real lapR, lapI; - if (step == 0) { - lapR = i>imin && ijmin && jimin && ijmin && jimin && ijmin && jimin && ijmin && j( m_rhs.dataPtr() ), - reinterpret_cast( m_rhs_fourier.dataPtr() ), - CUFFT_FORWARD); - if ( result != CUFFT_SUCCESS ) { - amrex::Print() << " forward transform using cufftExec failed ! Error: " << - CuFFTUtils::cufftErrorToString(result) << "\n"; - } -#elif defined(AMREX_USE_HIP) - amrex::Abort("Not implemented"); -#else - LaserFFT::VendorExecute( m_plan_fwd ); -#endif - - // Multiply by appropriate factors in Fourier space - amrex::Real dkx = 2.*MathConst::pi/geom.ProbLength(0); - amrex::Real dky = 2.*MathConst::pi/geom.ProbLength(1); - // acoeff_imag is supposed to be a nx*ny array. - // For the sake of simplicity, we evaluate it on-axis only. - const Complex acoeff = - step == 0 ? 6._rt/(c*dt*dz) - I * 4._rt * ( k0 + djn ) / (c*dt) : - 3._rt/(c*dt*dz) + 2._rt/(c*c*dt*dt) - I * 2._rt * ( k0 + djn ) / (c*dt); - amrex::ParallelFor( - bx, - [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - // divide rhs_fourier by -(k^2+a) - amrex::Real kx = (i 0. ? - 1._rt/(kx*kx + ky*ky + acoeff) : 0.; - rhs_fourier_arr(i,j,k,0) *= -inv_k2a; - }); - - // Transform rhs to Fourier space to get solution in sol -#ifdef AMREX_USE_CUDA - cudaStream_t stream_bkw = amrex::Gpu::Device::cudaStream(); - cufftSetStream ( m_plan_bkw, stream_bkw); - result = LaserFFT::VendorExecute( - m_plan_bkw, - reinterpret_cast( m_rhs_fourier.dataPtr() ), - reinterpret_cast( m_sol.dataPtr() ), - CUFFT_INVERSE); - if ( result != CUFFT_SUCCESS ) { - amrex::Print() << " forward transform using cufftExec failed ! Error: " << - CuFFTUtils::cufftErrorToString(result) << "\n"; - } -#elif defined(AMREX_USE_HIP) - amrex::Abort("Not implemented"); -#else - LaserFFT::VendorExecute( m_plan_bkw ); -#endif - - // Normalize and store solution in np1j00[0]. Guard cells are filled with 0s. - amrex::Box grown_bx = bx; - grown_bx.grow(m_slices_nguards); - const amrex::Real inv_numPts = 1./bx.numPts(); - amrex::ParallelFor( - grown_bx, - [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept { - if (i>=imin && i<=imax && j>=jmin && j<=jmax) { - np1j00_arr(i,j,0) = sol_arr(i,j,0).real() * inv_numPts; - np1j00_arr(i,j,1) = sol_arr(i,j,0).imag() * inv_numPts; - } else { - np1j00_arr(i,j,0) = 0._rt; - np1j00_arr(i,j,1) = 0._rt; - } - }); - } -} - -void -Laser::InitLaserSlice (const amrex::Geometry& geom, const int islice) -{ - if (!m_use_laser) return; - - HIPACE_PROFILE("Laser::InitLaserSlice()"); - - using namespace amrex::literals; - using Complex = amrex::GpuComplex; - - // Basic laser parameters and constants - Complex I(0,1); - constexpr int dcomp = 0; - const amrex::Real a0 = m_a0; - const amrex::Real k0 = 2._rt*MathConst::pi/m_lambda0; - const amrex::Real w0 = m_w0[0]; - const amrex::Real x0 = m_position_mean[0]; - const amrex::Real y0 = m_position_mean[1]; - const amrex::Real z0 = m_position_mean[2]; - const amrex::Real L0 = m_L0; - const amrex::Real zfoc = m_focal_distance; - - AMREX_ALWAYS_ASSERT(m_w0[0] == m_w0[1]); - AMREX_ALWAYS_ASSERT(x0 == 0._rt); - AMREX_ALWAYS_ASSERT(y0 == 0._rt); - - // Get grid properties - const auto plo = geom.ProbLoArray(); - amrex::Real const * const dx = geom.CellSize(); - const amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; - - // Envelope quantities common for all this slice - amrex::MultiFab& np1j00 = getSlices(WhichLaserSlice::np1j00); - -#ifdef AMREX_USE_OMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( amrex::MFIter mfi(np1j00, DfltMfiTlng); mfi.isValid(); ++mfi ){ - const amrex::Box& bx = mfi.tilebox(); - amrex::Array4 const & np1j00_arr = np1j00.array(mfi); - - // Initialize a Gaussian laser envelope on slice islice - amrex::ParallelFor( - bx, - [=] AMREX_GPU_DEVICE(int i, int j, int k) - { - amrex::Real z = plo[2] + (islice+0.5_rt)*dx_arr[2] - zfoc; - const amrex::Real x = (i+0.5_rt)*dx_arr[0]+plo[0]; - const amrex::Real y = (j+0.5_rt)*dx_arr[1]+plo[1]; - - // Compute envelope for time step 0 - Complex diffract_factor = 1._rt + I * z * 2._rt/( k0 * w0 * w0 ); - Complex inv_complex_waist_2 = 1._rt /( w0 * w0 * diffract_factor ); - Complex prefactor = a0/diffract_factor; - Complex time_exponent = (z-z0+zfoc)*(z-z0+zfoc)/(L0*L0); - Complex stcfactor = prefactor * amrex::exp( - time_exponent ); - Complex exp_argument = - ( x*x + y*y ) * inv_complex_waist_2; - Complex envelope = stcfactor * amrex::exp( exp_argument ); - np1j00_arr(i,j,k,dcomp ) = envelope.real(); - np1j00_arr(i,j,k,dcomp+1) = envelope.imag(); - } - ); - } + AMREX_ALWAYS_ASSERT(m_position_mean[0] == 0.); + AMREX_ALWAYS_ASSERT(m_position_mean[1] == 0.); } diff --git a/src/laser/MultiLaser.H b/src/laser/MultiLaser.H new file mode 100644 index 0000000000..12e108e361 --- /dev/null +++ b/src/laser/MultiLaser.H @@ -0,0 +1,244 @@ +/* Copyright 2022 + * + * This file is part of HiPACE++. + * + * Authors: MaxThevenet, AlexanderSinn + * Severin Diederichs, atmyers, Angel Ferran Pousa + * License: BSD-3-Clause-LBNL + */ + +#ifndef MULTILASER_H_ +#define MULTILASER_H_ + +#include "Laser.H" +#include "fields/Fields.H" +#include "mg_solver/HpMultiGrid.H" + +#include +#include +#include +#include + +#ifdef AMREX_USE_CUDA +# include +#elif defined(AMREX_USE_HIP) +# include +# include +#else +# include +#endif + +namespace LaserFFT { +#ifdef AMREX_USE_CUDA + using VendorFFT = cufftHandle; + const auto VendorCreate = cufftPlan2d; + const auto VendorDestroy = cufftDestroy; +# ifdef AMREX_USE_FLOAT + const auto VendorExecute = cufftExecC2C; + const auto cufft_type = CUFFT_C2C; + using cufftComplex = cuComplex; +# else + const auto VendorExecute = cufftExecZ2Z; + const auto cufft_type = CUFFT_Z2Z; + using cufftComplex = cuDoubleComplex; +# endif +#elif defined(AMREX_USE_HIP) + using VendorFFT = rocfft_plan; + const auto VendorDestroy = rocfft_plan_destroy; +// TODO +#else +# ifdef AMREX_USE_FLOAT + using VendorFFT = fftwf_plan; + using FFTWComplex = fftwf_complex; + const auto VendorCreate = fftwf_plan_dft_2d; + const auto VendorExecute = fftwf_execute; + const auto VendorDestroy = fftwf_destroy_plan; +# else + using VendorFFT = fftw_plan; + using FFTWComplex = fftw_complex; + const auto VendorCreate = fftw_plan_dft_2d; + const auto VendorExecute = fftw_execute; + const auto VendorDestroy = fftw_destroy_plan; +# endif +#endif +} + +/** \brief describes which slice with respect to the currently calculated is used */ +struct WhichLaserSlice { + // n stands for the time step, j for the longitudinal slice. + // n00 is time step n, nm1 is n-1 and np1 is n+1. Similar notation for slice j. + enum slice { nm1j00, nm1jp1, nm1jp2, n00j00, n00jp1, n00jp2, np1j00, np1jp1, np1jp2, N }; +}; + +class Fields; + +class MultiLaser +{ + + using SpectralFieldLoc = amrex::BaseFab >; + +private: + + static constexpr int m_nslices = WhichLaserSlice::N; + +public: + + /** Constructor */ + explicit MultiLaser () + { + ReadParameters(); + } + + ~MultiLaser () + { + if (!m_use_laser) return; + if (m_solver_type == "fft") { + LaserFFT::VendorDestroy( m_plan_fwd ); + LaserFFT::VendorDestroy( m_plan_bkw ); + } + } + + void ReadParameters (); + + /** get function for the 2D slices + * \param[in] islice slice index + */ + amrex::MultiFab& getSlices (int islice) {return m_slices[islice]; } + + /** get function for the 2D slices (const version) + * \param[in] islice slice index + */ + const amrex::MultiFab& getSlices (int islice) const {return m_slices[islice]; } + + /** Return the 3D FArrayBox containing the laser envelope, const version */ + const amrex::FArrayBox& getFAB () const {return m_F;} + + /** Retun the 3D FArrayBox containing the laser envelope, non-const version */ + amrex::FArrayBox& getFAB () {return m_F;} + + /** getter function, whether 3D envelope is store in host memory */ + int is3dOnHost () const {return m_3d_on_host;} + + /** \brief Allocate laser multifab + * \param[in] slice_ba box array of the slice + * \param[in] slice_dm corresponding distribution mapping + */ + void InitData (const amrex::BoxArray& slice_ba, + const amrex::DistributionMapping& slice_dm); + + /** \brief Initialize 3D laser field on current box. + * + * \param[in] step time step of the simulation + * \param[in] bx Box on which the laser field is initialized + * \param[in] gm Geometry of the problem + */ + void Init3DEnvelope (int step, amrex::Box bx, const amrex::Geometry& gm); + + /** \brief Copy from 2D slice on device to 3D array on host, and vice-versa + * + * \param[in] isl slice index, referring to the 3D slice + * \param[in] to3d if true, copy from 2D slice to 3D array. Otherwise, the other way. + */ + void Copy (int isl, bool to3d); + + /** Wrapper function to advance a laser slice by 1 time step. + * \param[in] fields Field object + * \param[in] geom Geometry object + * \param[in] dt time step of the simulation + * \param[in] step current iteration. Needed because step 0 needs a specific treatment. + */ + void AdvanceSlice (const Fields& fields, const amrex::Geometry& geom, amrex::Real dt, int step); + + /** Advance a laser slice by 1 time step using a multigrid solver. + * The complex phase of the envelope is evaluated on-axis only, but can be generalized to everywhere. + * + * \param[in] fields Field object + * \param[in] geom Geometry object + * \param[in] dt time step of the simulation + * \param[in] step current iteration. Needed because step 0 needs a specific treatment. + */ + void AdvanceSliceMG (const Fields& fields, const amrex::Geometry& geom, amrex::Real dt, int step); + + /** Advance a laser slice by 1 time step using a FFT solver. + * The complex phase of the envelope is evaluated on-axis only. + * + * \param[in] fields Field object + * \param[in] geom Geometry object + * \param[in] dt time step of the simulation + * \param[in] step current iteration. Needed because step 0 needs a specific treatment. + */ + void AdvanceSliceFFT (const Fields& fields, const amrex::Geometry& geom, amrex::Real dt, int step); + + /** Initialize 1 longitudinal slice of the laser, and store it in n00j00 (current time step) + * and nm1j00 (previous time step). + * + * \param[in] geom Geometry object for the slice + * \param[in] islice slice index + */ + void InitLaserSlice (const amrex::Geometry& geom, const int islice); + + bool m_use_laser {false}; /**< whether a laser is used or not */ + +private: + +/* + amrex::Gpu::DeviceVector getA0s (); + amrex::Gpu::DeviceVector getK0s (); + amrex::Gpu::DeviceVector getW0s (); + amrex::Gpu::DeviceVector getX0s (); + amrex::Gpu::DeviceVector getY0s (); + amrex::Gpu::DeviceVector getZ0s (); + amrex::Gpu::DeviceVector getL0s (); + amrex::Gpu::DeviceVector getZfocs (); +*/ + /** Laser central wavelength. + * he central wavelength influences the solver. As long as all the lasers are on the same grid + * (part of MultiLaser), this must be a property of MultiLaser. */ + amrex::Real m_lambda0 {0.}; + amrex::Vector m_names; /**< name of the laser */ + int m_nlasers; /**< Number of laser pulses */ + amrex::Vector m_all_lasers; /**< Each is a laser pulse */ + int m_3d_on_host {0};/** Whether the 3D laser envelope is stored in host or device memory */ + /** Number of guard cells for slices MultiFab */ + amrex::IntVect m_slices_nguards = {-1, -1, -1}; + std::string m_solver_type = "multigrid"; + bool m_use_phase {true}; + amrex::Box m_slice_box; + + /** Nb fields in 3D array: new_real, new_imag, old_real, old_imag */ + int m_nfields_3d {4}; + /** 3D laser data. Could be a vector over levels if MR needed for laser */ + amrex::FArrayBox m_F; + /** Array of N slices required to compute current slice */ + std::array m_slices; + amrex::Real m_MG_tolerance_rel = 1.e-4; + amrex::Real m_MG_tolerance_abs = 0.; + int m_MG_verbose = 0; + /** Whether to use time-averaged RHS in envelope solver. */ + bool m_MG_average_rhs = true; + /** hpmg solver for the envelope solver */ + std::unique_ptr m_mg; + + // Elements for the FFT-based laser envelope solver + // This could belong to AnyFFT etc., with 2 caveats: + // - This solver operates on a FArrayBox instead of a MultiFab, which is more adequate + // - The array in position space must be Complex rather than real, which takes up quite some + // rewriting, see https://github.com/MaxThevenet/hipace/tree/laser_solve, + // not sure what the best way to proceed is. + /** FFTW plan for forward C2C transform to solve Complex Poisson equation */ + LaserFFT::VendorFFT m_plan_fwd; + /** FFTW plan for backward C2C transform to solve Complex Poisson equation */ + LaserFFT::VendorFFT m_plan_bkw; + /** Complex FAB to store the solution (e.g. laser envelope on current slice) */ + SpectralFieldLoc m_sol; + /** Complex FAB to store the RHS in position space */ + SpectralFieldLoc m_rhs; + /** Complex FAB to store the RHS in Fourier space */ + SpectralFieldLoc m_rhs_fourier; +#ifdef AMREX_USE_CUDA + cufftResult m_result_fwd; + cufftResult m_result_bkw; +#endif +}; + +#endif // MULTILASER_H_ diff --git a/src/laser/MultiLaser.cpp b/src/laser/MultiLaser.cpp new file mode 100644 index 0000000000..8afde13647 --- /dev/null +++ b/src/laser/MultiLaser.cpp @@ -0,0 +1,830 @@ +/* Copyright 2022 + * + * This file is part of HiPACE++. + * + * Authors: MaxThevenet, AlexanderSinn + * Severin Diederichs, atmyers, Angel Ferran Pousa + * License: BSD-3-Clause-LBNL + */ + +#include "MultiLaser.H" +#include "utils/Constants.H" +#include "Hipace.H" +#include "utils/HipaceProfilerWrapper.H" +#ifdef AMREX_USE_CUDA +# include "fields/fft_poisson_solver/fft/CuFFTUtils.H" +#endif +#include + +#ifdef AMREX_USE_CUDA +# include +#elif defined(AMREX_USE_HIP) +# include +# include +#else +# include +#endif + +void +MultiLaser::ReadParameters () +{ + amrex::ParmParse pp("lasers"); + getWithParser(pp, "names", m_names); + + m_use_laser = m_names[0] != "no_laser"; + + if (!m_use_laser) return; +#if defined(AMREX_USE_HIP) + amrex::Abort("Laser solver not implemented with HIP"); +#endif + + m_nlasers = m_names.size(); + for (int i = 0; i < m_nlasers; ++i) { + m_all_lasers.emplace_back(Laser(m_names[i])); + } + + getWithParser(pp, "lambda0", m_lambda0); + queryWithParser(pp, "3d_on_host", m_3d_on_host); + queryWithParser(pp, "use_phase", m_use_phase); + queryWithParser(pp, "solver_type", m_solver_type); + AMREX_ALWAYS_ASSERT(m_solver_type == "multigrid" || m_solver_type == "fft"); + + bool mg_param_given = queryWithParser(pp, "MG_tolerance_rel", m_MG_tolerance_rel); + mg_param_given += queryWithParser(pp, "MG_tolerance_abs", m_MG_tolerance_abs); + mg_param_given += queryWithParser(pp, "MG_verbose", m_MG_verbose); + mg_param_given += queryWithParser(pp, "MG_average_rhs", m_MG_average_rhs); + + // Raise warning if user specifies MG parameters without using the MG solver + if (mg_param_given && (m_solver_type != "multigrid")) { + amrex::Print()<<"WARNING: parameters laser.MG_... only active if laser.solver_type = multigrid\n"; + } +} + + +void +MultiLaser::InitData (const amrex::BoxArray& slice_ba, + const amrex::DistributionMapping& slice_dm) +{ + if (!m_use_laser) return; + + HIPACE_PROFILE("MultiLaser::InitData()"); + + // Alloc 2D slices + // 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}; + AMREX_ALWAYS_ASSERT(WhichLaserSlice::N == m_nslices); + for (int islice=0; islice(m_rhs.dataPtr()), + reinterpret_cast(m_rhs_fourier.dataPtr()), + FFTW_FORWARD, FFTW_ESTIMATE); + // Backward FFT plan + m_plan_bkw = LaserFFT::VendorCreate( + fft_size[1], fft_size[0], + reinterpret_cast(m_rhs_fourier.dataPtr()), + reinterpret_cast(m_sol.dataPtr()), + FFTW_BACKWARD, FFTW_ESTIMATE); +#endif + } +} + +void +MultiLaser::Init3DEnvelope (int step, amrex::Box bx, const amrex::Geometry& gm) +{ + + if (!m_use_laser) return; + + HIPACE_PROFILE("MultiLaser::Init3DEnvelope()"); + // Allocate the 3D field on this box + // Note: box has no guard cells + m_F.resize(bx, m_nfields_3d, m_3d_on_host ? amrex::The_Pinned_Arena() : amrex::The_Arena()); + + if (step > 0) return; + + // In order to use the normal Copy function, we use slice np1j00 as a tmp array + // to initialize the laser in the loop over slices below. + // We need to keep the value in np1j00, as it is shifted to np1jp1 and used to compute + // the following slices. This is relevant for the first slices at step 0 of every box + // (except for the head box). + amrex::FArrayBox store_np1j00; + store_np1j00.resize(m_slice_box, 2, amrex::The_Arena()); + store_np1j00.copy(m_slices[WhichLaserSlice::np1j00][0]); + + // Loop over slices + for (int isl = bx.bigEnd(Direction::z); isl >= bx.smallEnd(Direction::z); --isl){ + // Compute initial field on the current (device) slice np1j00 + InitLaserSlice(gm, isl); + // Copy: (device) slice np1j00 to the right location + // in the (host) 3D array of the current time. + Copy(isl, true); + } + + // Reset np1j00 to its original value. + m_slices[WhichLaserSlice::np1j00][0].copy(store_np1j00); +} + +void +MultiLaser::Copy (int isl, bool to3d) +{ + if (!m_use_laser) return; + + using namespace amrex::literals; + + HIPACE_PROFILE("MultiLaser::Copy()"); + + amrex::MultiFab& nm1j00 = m_slices[WhichLaserSlice::nm1j00]; + amrex::MultiFab& nm1jp1 = m_slices[WhichLaserSlice::nm1jp1]; + amrex::MultiFab& nm1jp2 = m_slices[WhichLaserSlice::nm1jp2]; + amrex::MultiFab& n00j00 = m_slices[WhichLaserSlice::n00j00]; + amrex::MultiFab& n00jp1 = m_slices[WhichLaserSlice::n00jp1]; + amrex::MultiFab& n00jp2 = m_slices[WhichLaserSlice::n00jp2]; + amrex::MultiFab& np1j00 = m_slices[WhichLaserSlice::np1j00]; + amrex::MultiFab& np1jp1 = m_slices[WhichLaserSlice::np1jp1]; + amrex::MultiFab& np1jp2 = m_slices[WhichLaserSlice::np1jp2]; + + for ( amrex::MFIter mfi(n00j00, DfltMfi); mfi.isValid(); ++mfi ){ + const amrex::Box& bx = mfi.tilebox(); + Array3 nm1j00_arr = nm1j00.array(mfi); + Array3 nm1jp1_arr = nm1jp1.array(mfi); + Array3 nm1jp2_arr = nm1jp2.array(mfi); + Array3 n00j00_arr = n00j00.array(mfi); + Array3 n00jp1_arr = n00jp1.array(mfi); + Array3 n00jp2_arr = n00jp2.array(mfi); + Array3 np1j00_arr = np1j00.array(mfi); + Array3 np1jp1_arr = np1jp1.array(mfi); + Array3 np1jp2_arr = np1jp2.array(mfi); + amrex::Array4 host_arr = m_F.array(); + amrex::ParallelFor( + bx, 2, + [=] AMREX_GPU_DEVICE(int i, int j, int, int n) noexcept + { + // +2 in 3D array means old. + // 2 components for complex numbers. + if (to3d){ + // this slice into old host + host_arr(i,j,isl,n+2) = n00j00_arr(i,j,n); + // next time slice into new host + host_arr(i,j,isl,n ) = np1j00_arr(i,j,n); + } else { + // Shift slices of step n-1, and get current slice from 3D array + nm1jp2_arr(i,j,n) = nm1jp1_arr(i,j,n); + nm1jp1_arr(i,j,n) = nm1j00_arr(i,j,n); + nm1j00_arr(i,j,n) = host_arr(i,j,isl,n+2); + // Shift slices of step n, and get current slice from 3D array + n00jp2_arr(i,j,n) = n00jp1_arr(i,j,n); + n00jp1_arr(i,j,n) = n00j00_arr(i,j,n); + n00j00_arr(i,j,n) = host_arr(i,j,isl,n); + // Shift slices of step n+1. Current slice will be computed + np1jp2_arr(i,j,n) = np1jp1_arr(i,j,n); + np1jp1_arr(i,j,n) = np1j00_arr(i,j,n); + } + }); + } +} + +void +MultiLaser::AdvanceSlice (const Fields& fields, const amrex::Geometry& geom, amrex::Real dt, int step) +{ + + if (!m_use_laser) return; + + if (m_solver_type == "multigrid") { + AdvanceSliceMG(fields, geom, dt, step); + } else if (m_solver_type == "fft") { + AdvanceSliceFFT(fields, geom, dt, step); + } else { + amrex::Abort("laser.solver_type must be fft or multigrid"); + } +} + +void +MultiLaser::AdvanceSliceMG (const Fields& fields, const amrex::Geometry& geom, amrex::Real dt, int step) +{ + + HIPACE_PROFILE("MultiLaser::AdvanceSliceMG()"); + + using namespace amrex::literals; + using Complex = amrex::GpuComplex; + constexpr Complex I(0.,1.); + + const amrex::Real dx = geom.CellSize(0); + const amrex::Real dy = geom.CellSize(1); + const amrex::Real dz = geom.CellSize(2); + + const PhysConst phc = get_phys_const(); + const amrex::Real c = phc.c; + const amrex::Real k0 = 2.*MathConst::pi/m_lambda0; + const bool do_avg_rhs = m_MG_average_rhs; + + amrex::MultiFab& nm1j00 = m_slices[WhichLaserSlice::nm1j00]; + amrex::MultiFab& nm1jp1 = m_slices[WhichLaserSlice::nm1jp1]; + amrex::MultiFab& nm1jp2 = m_slices[WhichLaserSlice::nm1jp2]; + amrex::MultiFab& n00j00 = m_slices[WhichLaserSlice::n00j00]; + amrex::MultiFab& n00jp1 = m_slices[WhichLaserSlice::n00jp1]; + amrex::MultiFab& n00jp2 = m_slices[WhichLaserSlice::n00jp2]; + amrex::MultiFab& np1j00 = m_slices[WhichLaserSlice::np1j00]; + amrex::MultiFab& np1jp1 = m_slices[WhichLaserSlice::np1jp1]; + amrex::MultiFab& np1jp2 = m_slices[WhichLaserSlice::np1jp2]; + + amrex::FArrayBox rhs_mg; + amrex::FArrayBox acoeff_real; + amrex::Real acoeff_real_scalar = 0._rt; + amrex::Real acoeff_imag_scalar = 0._rt; + + amrex::Real djn {0.}; + + for ( amrex::MFIter mfi(n00j00, DfltMfi); mfi.isValid(); ++mfi ){ + const amrex::Box& bx = mfi.tilebox(); + const int imin = bx.smallEnd(0); + const int imax = bx.bigEnd (0); + const int jmin = bx.smallEnd(1); + const int jmax = bx.bigEnd (1); + + acoeff_real.resize(bx, 1, amrex::The_Arena()); + rhs_mg.resize(bx, 2, amrex::The_Arena()); + Array3 nm1j00_arr = nm1j00.array(mfi); + Array3 nm1jp1_arr = nm1jp1.array(mfi); + Array3 nm1jp2_arr = nm1jp2.array(mfi); + Array3 n00j00_arr = n00j00.array(mfi); + Array3 n00jp1_arr = n00jp1.array(mfi); + Array3 n00jp2_arr = n00jp2.array(mfi); + Array3 np1jp1_arr = np1jp1.array(mfi); + Array3 np1jp2_arr = np1jp2.array(mfi); + Array3 rhs_mg_arr = rhs_mg.array(); + Array3 acoeff_real_arr = acoeff_real.array(); + Array3 rhs_arr = m_rhs.array(); + + constexpr int lev = 0; + const amrex::FArrayBox& isl_fab = fields.getSlices(lev, WhichSlice::This)[mfi]; + Array3 const isl_arr = isl_fab.array(); + const int chi = Comps[WhichSlice::This]["chi"]; + + // Calculate phase terms. 0 if !m_use_phase + amrex::Real tj00 = 0.; + amrex::Real tjp1 = 0.; + amrex::Real tjp2 = 0.; + + if (m_use_phase) { + int const Nx = bx.length(0); + int const Ny = bx.length(1); + + // Get the central point. + int const imid = (Nx+1)/2; + int const jmid = (Ny+1)/2; + + // Calculate complex arguments (theta) needed + // Just once, on axis, as done in Wake-T + // This is done with a reduce operation, returning the sum of the four elements nearest + // the axis (both real and imag parts, and for the 3 arrays relevant) ... + amrex::ReduceOps< + amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum, + amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum> reduce_op; + amrex::ReduceData< + amrex::Real, amrex::Real, amrex::Real, + amrex::Real, amrex::Real, amrex::Real> reduce_data(reduce_op); + using ReduceTuple = typename decltype(reduce_data)::Type; + reduce_op.eval(bx, reduce_data, + [=] AMREX_GPU_DEVICE (int i, int j, int) -> ReduceTuple + { + if ( ( i == imid-1 || i == imid ) && ( j == jmid-1 || j == jmid ) ) { + return { + n00j00_arr(i,j,0), n00j00_arr(i,j,1), + n00jp1_arr(i,j,0), n00jp1_arr(i,j,1), + n00jp2_arr(i,j,0), n00jp2_arr(i,j,1) + }; + } else { + return {0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt}; + } + }); + // ... and taking the argument of the resulting complex number. + ReduceTuple hv = reduce_data.value(reduce_op); + tj00 = std::atan2(amrex::get<1>(hv), amrex::get<0>(hv)); + tjp1 = std::atan2(amrex::get<3>(hv), amrex::get<2>(hv)); + tjp2 = std::atan2(amrex::get<5>(hv), amrex::get<4>(hv)); + } + + amrex::Real dt1 = tj00 - tjp1; + amrex::Real dt2 = tjp1 - tjp2; + if (dt1 <-1.5_rt*MathConst::pi) dt1 += 2._rt*MathConst::pi; + if (dt1 > 1.5_rt*MathConst::pi) dt1 -= 2._rt*MathConst::pi; + if (dt2 <-1.5_rt*MathConst::pi) dt2 += 2._rt*MathConst::pi; + if (dt2 > 1.5_rt*MathConst::pi) dt2 -= 2._rt*MathConst::pi; + Complex exp1 = amrex::exp(I*(tj00-tjp1)); + Complex exp2 = amrex::exp(I*(tj00-tjp2)); + + // D_j^n as defined in Benedetti's 2017 paper + djn = ( -3._rt*dt1 + dt2 ) / (2._rt*dz); + acoeff_real_scalar = step == 0 ? 6._rt/(c*dt*dz) + : 3._rt/(c*dt*dz) + 2._rt/(c*c*dt*dt); + acoeff_imag_scalar = step == 0 ? -4._rt * ( k0 + djn ) / (c*dt) + : -2._rt * ( k0 + djn ) / (c*dt); + + amrex::ParallelFor( + bx, 1, + [=] AMREX_GPU_DEVICE(int i, int j, int, int) noexcept + { + // Transverse Laplacian of real and imaginary parts of A_j^n-1 + amrex::Real lapR, lapI; + if (step == 0) { + lapR = i>imin && ijmin && jimin && ijmin && jimin && ijmin && jimin && ijmin && j(slice_geom); + } + + const int max_iters = 200; + m_mg->solve2(np1j00[0], rhs_mg, acoeff_real, acoeff_imag_scalar, + m_MG_tolerance_rel, m_MG_tolerance_abs, max_iters, m_MG_verbose); +} + +void +MultiLaser::AdvanceSliceFFT (const Fields& fields, const amrex::Geometry& geom, const amrex::Real dt, int step) +{ + + HIPACE_PROFILE("MultiLaser::AdvanceSliceFFT()"); + + using namespace amrex::literals; + using Complex = amrex::GpuComplex; + constexpr Complex I(0.,1.); + + const amrex::Real dx = geom.CellSize(0); + const amrex::Real dy = geom.CellSize(1); + const amrex::Real dz = geom.CellSize(2); + + const PhysConst phc = get_phys_const(); + const amrex::Real c = phc.c; + const amrex::Real k0 = 2.*MathConst::pi/m_lambda0; + + amrex::MultiFab& nm1j00 = m_slices[WhichLaserSlice::nm1j00]; + amrex::MultiFab& nm1jp1 = m_slices[WhichLaserSlice::nm1jp1]; + amrex::MultiFab& nm1jp2 = m_slices[WhichLaserSlice::nm1jp2]; + amrex::MultiFab& n00j00 = m_slices[WhichLaserSlice::n00j00]; + amrex::MultiFab& n00jp1 = m_slices[WhichLaserSlice::n00jp1]; + amrex::MultiFab& n00jp2 = m_slices[WhichLaserSlice::n00jp2]; + amrex::MultiFab& np1j00 = m_slices[WhichLaserSlice::np1j00]; + amrex::MultiFab& np1jp1 = m_slices[WhichLaserSlice::np1jp1]; + amrex::MultiFab& np1jp2 = m_slices[WhichLaserSlice::np1jp2]; + + for ( amrex::MFIter mfi(n00j00, DfltMfi); mfi.isValid(); ++mfi ){ + const amrex::Box& bx = mfi.tilebox(); + const int imin = bx.smallEnd(0); + const int imax = bx.bigEnd (0); + const int jmin = bx.smallEnd(1); + const int jmax = bx.bigEnd (1); + + // solution: complex array + // The right-hand side is computed and stored in rhs + // Then rhs is Fourier-transformed into rhs_fourier, then multiplied by -1/(k**2+a) + // rhs_fourier is FFT-back-transformed to sol, and sol is normalized and copied into np1j00. + Array3 sol_arr = m_sol.array(); + Array3 rhs_arr = m_rhs.array(); + amrex::Array4 rhs_fourier_arr = m_rhs_fourier.array(); + + Array3 nm1j00_arr = nm1j00.array(mfi); + Array3 nm1jp1_arr = nm1jp1.array(mfi); + Array3 nm1jp2_arr = nm1jp2.array(mfi); + Array3 n00j00_arr = n00j00.array(mfi); + Array3 n00jp1_arr = n00jp1.array(mfi); + Array3 n00jp2_arr = n00jp2.array(mfi); + Array3 np1j00_arr = np1j00.array(mfi); + Array3 np1jp1_arr = np1jp1.array(mfi); + Array3 np1jp2_arr = np1jp2.array(mfi); + + constexpr int lev = 0; + const amrex::FArrayBox& isl_fab = fields.getSlices(lev, WhichSlice::This)[mfi]; + Array3 const isl_arr = isl_fab.array(); + const int chi = Comps[WhichSlice::This]["chi"]; + + int const Nx = bx.length(0); + int const Ny = bx.length(1); + + // Get the central point. Useful to get the on-axis phase and calculate kx and ky. + int const imid = (Nx+1)/2; + int const jmid = (Ny+1)/2; + + // Calculate phase terms. 0 if !m_use_phase + amrex::Real tj00 = 0.; + amrex::Real tjp1 = 0.; + amrex::Real tjp2 = 0.; + + if (m_use_phase) { + // Calculate complex arguments (theta) needed + // Just once, on axis, as done in Wake-T + // This is done with a reduce operation, returning the sum of the four elements nearest + // the axis (both real and imag parts, and for the 3 arrays relevant) ... + amrex::ReduceOps< + amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum, + amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum> reduce_op; + amrex::ReduceData< + amrex::Real, amrex::Real, amrex::Real, + amrex::Real, amrex::Real, amrex::Real> reduce_data(reduce_op); + using ReduceTuple = typename decltype(reduce_data)::Type; + reduce_op.eval(bx, reduce_data, + [=] AMREX_GPU_DEVICE (int i, int j, int) -> ReduceTuple + { + if ( ( i == imid-1 || i == imid ) && ( j == jmid-1 || j == jmid ) ) { + return { + n00j00_arr(i,j,0), n00j00_arr(i,j,1), + n00jp1_arr(i,j,0), n00jp1_arr(i,j,1), + n00jp2_arr(i,j,0), n00jp2_arr(i,j,1) + }; + } else { + return {0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt}; + } + }); + // ... and taking the argument of the resulting complex number. + ReduceTuple hv = reduce_data.value(reduce_op); + tj00 = std::atan2(amrex::get<1>(hv), amrex::get<0>(hv)); + tjp1 = std::atan2(amrex::get<3>(hv), amrex::get<2>(hv)); + tjp2 = std::atan2(amrex::get<5>(hv), amrex::get<4>(hv)); + } + + amrex::Real dt1 = tj00 - tjp1; + amrex::Real dt2 = tjp1 - tjp2; + if (dt1 <-1.5_rt*MathConst::pi) dt1 += 2._rt*MathConst::pi; + if (dt1 > 1.5_rt*MathConst::pi) dt1 -= 2._rt*MathConst::pi; + if (dt2 <-1.5_rt*MathConst::pi) dt2 += 2._rt*MathConst::pi; + if (dt2 > 1.5_rt*MathConst::pi) dt2 -= 2._rt*MathConst::pi; + Complex exp1 = amrex::exp(I*(tj00-tjp1)); + Complex exp2 = amrex::exp(I*(tj00-tjp2)); + + // D_j^n as defined in Benedetti's 2017 paper + amrex::Real djn = ( -3._rt*dt1 + dt2 ) / (2._rt*dz); + amrex::ParallelFor( + bx, 1, + [=] AMREX_GPU_DEVICE(int i, int j, int, int) noexcept + { + // Transverse Laplacian of real and imaginary parts of A_j^n-1 + amrex::Real lapR, lapI; + if (step == 0) { + lapR = i>imin && ijmin && jimin && ijmin && jimin && ijmin && jimin && ijmin && j( m_rhs.dataPtr() ), + reinterpret_cast( m_rhs_fourier.dataPtr() ), + CUFFT_FORWARD); + if ( result != CUFFT_SUCCESS ) { + amrex::Print() << " forward transform using cufftExec failed ! Error: " << + CuFFTUtils::cufftErrorToString(result) << "\n"; + } +#elif defined(AMREX_USE_HIP) + amrex::Abort("Not implemented"); +#else + LaserFFT::VendorExecute( m_plan_fwd ); +#endif + + // Multiply by appropriate factors in Fourier space + amrex::Real dkx = 2.*MathConst::pi/geom.ProbLength(0); + amrex::Real dky = 2.*MathConst::pi/geom.ProbLength(1); + // acoeff_imag is supposed to be a nx*ny array. + // For the sake of simplicity, we evaluate it on-axis only. + const Complex acoeff = + step == 0 ? 6._rt/(c*dt*dz) - I * 4._rt * ( k0 + djn ) / (c*dt) : + 3._rt/(c*dt*dz) + 2._rt/(c*c*dt*dt) - I * 2._rt * ( k0 + djn ) / (c*dt); + amrex::ParallelFor( + bx, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + // divide rhs_fourier by -(k^2+a) + amrex::Real kx = (i 0. ? + 1._rt/(kx*kx + ky*ky + acoeff) : 0.; + rhs_fourier_arr(i,j,k,0) *= -inv_k2a; + }); + + // Transform rhs to Fourier space to get solution in sol +#ifdef AMREX_USE_CUDA + cudaStream_t stream_bkw = amrex::Gpu::Device::cudaStream(); + cufftSetStream ( m_plan_bkw, stream_bkw); + result = LaserFFT::VendorExecute( + m_plan_bkw, + reinterpret_cast( m_rhs_fourier.dataPtr() ), + reinterpret_cast( m_sol.dataPtr() ), + CUFFT_INVERSE); + if ( result != CUFFT_SUCCESS ) { + amrex::Print() << " forward transform using cufftExec failed ! Error: " << + CuFFTUtils::cufftErrorToString(result) << "\n"; + } +#elif defined(AMREX_USE_HIP) + amrex::Abort("Not implemented"); +#else + LaserFFT::VendorExecute( m_plan_bkw ); +#endif + + // Normalize and store solution in np1j00[0]. Guard cells are filled with 0s. + amrex::Box grown_bx = bx; + grown_bx.grow(m_slices_nguards); + const amrex::Real inv_numPts = 1./bx.numPts(); + amrex::ParallelFor( + grown_bx, + [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept { + if (i>=imin && i<=imax && j>=jmin && j<=jmax) { + np1j00_arr(i,j,0) = sol_arr(i,j,0).real() * inv_numPts; + np1j00_arr(i,j,1) = sol_arr(i,j,0).imag() * inv_numPts; + } else { + np1j00_arr(i,j,0) = 0._rt; + np1j00_arr(i,j,1) = 0._rt; + } + }); + } +} + +void +MultiLaser::InitLaserSlice (const amrex::Geometry& geom, const int islice) +{ + if (!m_use_laser) return; + + HIPACE_PROFILE("MultiLaser::InitLaserSlice()"); + + using namespace amrex::literals; + using Complex = amrex::GpuComplex; + + // Basic laser parameters and constants + Complex I(0,1); + constexpr int dcomp = 0; + const amrex::Real k0 = 2._rt*MathConst::pi/m_lambda0; + + // Get grid properties + const auto plo = geom.ProbLoArray(); + amrex::Real const * const dx = geom.CellSize(); + const amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; + + // Envelope quantities common for all this slice + amrex::MultiFab& np1j00 = getSlices(WhichLaserSlice::np1j00); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( amrex::MFIter mfi(np1j00, DfltMfiTlng); mfi.isValid(); ++mfi ){ + const amrex::Box& bx = mfi.tilebox(); + amrex::Array4 const & np1j00_arr = np1j00.array(mfi); + + // Initialize a Gaussian laser envelope on slice islice + for (int ilaser=0; ilaser +MultiLaser::getA0s () +{ + amrex::Gpu::DeviceVector a0s(m_nlasers); + for (int i=0; i +MultiLaser::getK0s () +{ + amrex::Gpu::DeviceVector k0s(m_nlasers); + for (int i=0; i +MultiLaser::getW0s () +{ + amrex::Gpu::DeviceVector w0s(m_nlasers); + for (int i=0; i +MultiLaser::getX0s () +{ + amrex::Gpu::DeviceVector x0s(m_nlasers); + for (int i=0; i +MultiLaser::getY0s () +{ + amrex::Gpu::DeviceVector y0s(m_nlasers); + for (int i=0; i +MultiLaser::getX0s () +{ + amrex::Gpu::DeviceVector x0s(m_nlasers); + for (int i=0; i getX0s (); + amrex::Gpu::DeviceVector getY0s (); + amrex::Gpu::DeviceVector getZ0s (); + amrex::Gpu::DeviceVector getL0s (); + amrex::Gpu::DeviceVector getZfocs (); +*/ diff --git a/src/particles/deposition/ExplicitDeposition.H b/src/particles/deposition/ExplicitDeposition.H index 1e0275bd7b..3ab6723d4c 100644 --- a/src/particles/deposition/ExplicitDeposition.H +++ b/src/particles/deposition/ExplicitDeposition.H @@ -14,12 +14,12 @@ /** Depose Sx and Sy of particles in species plasma into the current 2D slice in fields * \param[in] plasma species of which the current is deposited * \param[in,out] fields the general field class, modified by this function - * \param[in] laser that affects the plasma during the deposition + * \param[in] multi_laser Lasers that affects the plasma during the deposition * \param[in] gm Geometry of the simulation, to get the cell size etc. * \param[in] lev MR level */ void -ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, const Laser& laser, +ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, const MultiLaser& multi_laser, amrex::Geometry const& gm, const int lev); #endif // EXPLICITDEPOSITION_H_ diff --git a/src/particles/deposition/ExplicitDeposition.cpp b/src/particles/deposition/ExplicitDeposition.cpp index f22deead7e..67f7c1fae0 100644 --- a/src/particles/deposition/ExplicitDeposition.cpp +++ b/src/particles/deposition/ExplicitDeposition.cpp @@ -16,7 +16,7 @@ #include "AMReX_GpuLaunch.H" void -ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, const Laser& laser, +ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, const MultiLaser& multi_laser, amrex::Geometry const& gm, const int lev) { HIPACE_PROFILE("ExplicitDeposition()"); using namespace amrex::literals; @@ -48,8 +48,8 @@ ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, const Laser plasma.m_can_ionize ? soa.GetIntData(PlasmaIdx::ion_lev).data() : nullptr; // Construct empty Array4 with one z slice so that Array3 constructor works for no laser - const Array3 a_laser_arr = laser.m_use_laser ? - laser.getSlices(WhichLaserSlice::n00j00).const_array(pti) : + const Array3 a_laser_arr = multi_laser.m_use_laser ? + multi_laser.getSlices(WhichLaserSlice::n00j00).const_array(pti) : amrex::Array4(nullptr, {0,0,0}, {0,0,1}, 0); const amrex::Real x_pos_offset = GetPosOffset(0, gm, isl_fab.box()); @@ -71,7 +71,7 @@ ExplicitDeposition (PlasmaParticleContainer& plasma, Fields& fields, const Laser amrex::TypeList, // depos_order amrex::CompileTimeOptions, // can_ionize amrex::CompileTimeOptions>{}, // use_laser - {Hipace::m_depos_order_xy, plasma.m_can_ionize, laser.m_use_laser}, + {Hipace::m_depos_order_xy, plasma.m_can_ionize, multi_laser.m_use_laser}, pti.numParticles(), [=] AMREX_GPU_DEVICE (int ip, auto a_depos_order, auto can_ionize, auto use_laser) { constexpr int depos_order = a_depos_order.value; diff --git a/src/particles/deposition/PlasmaDepositCurrent.H b/src/particles/deposition/PlasmaDepositCurrent.H index 404717fa79..a68da95e8f 100644 --- a/src/particles/deposition/PlasmaDepositCurrent.H +++ b/src/particles/deposition/PlasmaDepositCurrent.H @@ -12,13 +12,13 @@ #include "particles/sorting/TileSort.H" #include "fields/Fields.H" #include "utils/Constants.H" -#include "laser/Laser.H" +#include "laser/MultiLaser.H" #include "Hipace.H" /** Depose current of particles in species plasma into the current 2D slice in fields * \param[in] plasma species of which the current is deposited * \param[in,out] fields the general field class, modified by this function - * \param[in] laser that affects the plasma during the deposition + * \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 @@ -31,7 +31,7 @@ * \param[in] bin_size tile size (square) */ void -DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const Laser& laser, +DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const MultiLaser& multi_laser, const int which_slice, const bool temp_slice, const bool deposit_jx_jy, const bool deposit_jz, const bool deposit_rho, bool deposit_chi, amrex::Geometry const& gm, int const lev, diff --git a/src/particles/deposition/PlasmaDepositCurrent.cpp b/src/particles/deposition/PlasmaDepositCurrent.cpp index 15b785bf42..2ccc6f44fb 100644 --- a/src/particles/deposition/PlasmaDepositCurrent.cpp +++ b/src/particles/deposition/PlasmaDepositCurrent.cpp @@ -20,7 +20,7 @@ void -DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const Laser& laser, +DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const MultiLaser& multi_laser, const int which_slice, const bool temp_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, @@ -62,7 +62,7 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const Laser& l amrex::Vector& tmp_dens = fields.getTmpDensities(); // extract the laser Fields - const amrex::MultiFab& a_mf = laser.getSlices(WhichLaserSlice::n00j00); + const amrex::MultiFab& a_mf = multi_laser.getSlices(WhichLaserSlice::n00j00); // Offset for converting positions to indexes const amrex::Real x_pos_offset = GetPosOffset(0, gm, isl_fab.box()); @@ -86,7 +86,7 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const Laser& l // Extract laser array from MultiFab const Array3 a_laser_arr = - laser.m_use_laser ? a_mf[pti].const_array() : amrex::Array4(); + multi_laser.m_use_laser ? a_mf[pti].const_array() : amrex::Array4(); // Extract box properties const amrex::Real invvol = Hipace::m_normalized_units ? 1._rt : 1._rt/(dx[0]*dx[1]*dx[2]); @@ -174,7 +174,7 @@ DepositCurrent (PlasmaParticleContainer& plasma, Fields & fields, const Laser& l Hipace::m_outer_depos_loop, plasma.m_can_ionize, do_tiling, - laser.m_use_laser + multi_laser.m_use_laser }, num_particles, [=] AMREX_GPU_DEVICE (int idx, auto depos_order, auto outer_depos_loop, diff --git a/src/particles/plasma/MultiPlasma.H b/src/particles/plasma/MultiPlasma.H index c83949fe8f..17d29fbee0 100644 --- a/src/particles/plasma/MultiPlasma.H +++ b/src/particles/plasma/MultiPlasma.H @@ -11,7 +11,7 @@ #include "PlasmaParticleContainer.H" #include "particles/sorting/TileSort.H" #include "fields/Fields.H" -#include "laser/Laser.H" +#include "laser/MultiLaser.H" #include "particles/collisions/CoulombCollision.H" class MultiPlasma @@ -42,7 +42,7 @@ public: /** Loop over plasma species and depose their currents into the current 2D slice in fields * * \param[in,out] fields the general field class, modified by this function - * \param[in] laser that affects the plasma during the deposition + * \param[in] multi_laser Lasers 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 @@ -53,18 +53,18 @@ public: * \param[in] lev MR level */ void DepositCurrent ( - Fields & fields, const Laser & laser, int which_slice, bool temp_slice, bool deposit_jx_jy, + Fields & fields, const MultiLaser & multi_laser, int which_slice, bool temp_slice, bool deposit_jx_jy, bool deposit_jz, bool deposit_rho, bool deposit_chi, amrex::Geometry const& gm, int const lev); /** Loop over plasma species and depose Sx and Sy into the current 2D slice in fields * * \param[in,out] fields the general field class, modified by this function - * \param[in] laser that affects the plasma during the deposition + * \param[in] multi_laser Lasers that affects the plasma during the deposition * \param[in] gm Geometry of the simulation, to get the cell size etc. * \param[in] lev MR level */ - void ExplicitDeposition (Fields& fields, const Laser& laser, + void ExplicitDeposition (Fields& fields, const MultiLaser& multi_laser, amrex::Geometry const& gm, int const lev); /** \brief Return max density, to compute the adaptive time step. @@ -77,7 +77,7 @@ public: /** \brief Loop over plasma species and Gather fields, update forces and push particles * * \param[in,out] fields the general field class, modified by this function - * \param[in] laser that affects the plasma during the deposition + * \param[in] multi_laser Lasers that affects the plasma during the deposition * \param[in] gm Geometry of the simulation, to get the cell size etc. * \param[in] temp_slice if true, the temporary data (x_temp, ...) will be used * \param[in] do_push boolean to define if plasma particles are pushed @@ -86,7 +86,7 @@ public: * \param[in] lev MR level */ void AdvanceParticles ( - const Fields & fields, const Laser & laser, amrex::Geometry const& gm, bool temp_slice, + const Fields & fields, const MultiLaser & multi_laser, amrex::Geometry const& gm, bool temp_slice, bool do_push, bool do_update, bool do_shift, int lev); /** \brief Resets the particle position x, y, to x_prev, y_prev @@ -99,13 +99,13 @@ public: /** \brief Loop over plasma species and deposit their neutralizing background, if needed * * \param[in,out] fields the general field class, modified by this function - * \param[in] laser that affects the plasma during the deposition + * \param[in] multi_laser that affects the plasma during the deposition * \param[in] which_slice slice in which the densities are deposited * \param[in] gm Geometry of the simulation, to get the cell size etc. * \param[in] nlev number of MR levels */ void DepositNeutralizingBackground ( - Fields & fields, const Laser & laser, int which_slice, amrex::Geometry const& gm, + Fields & fields, const MultiLaser & multi_laser, int which_slice, amrex::Geometry const& gm, int const nlev); /** Calculates Ionization Probability and makes new Plasma Particles diff --git a/src/particles/plasma/MultiPlasma.cpp b/src/particles/plasma/MultiPlasma.cpp index edb6f66404..b519069dab 100644 --- a/src/particles/plasma/MultiPlasma.cpp +++ b/src/particles/plasma/MultiPlasma.cpp @@ -82,34 +82,34 @@ MultiPlasma::maxDensity () const void MultiPlasma::DepositCurrent ( - Fields & fields, const Laser & laser, int which_slice, bool temp_slice, bool deposit_jx_jy, - bool deposit_jz, bool deposit_rho, bool deposit_chi, amrex::Geometry const& gm, - int const lev) + Fields & fields, const MultiLaser & multi_laser, int which_slice, bool temp_slice, + bool deposit_jx_jy, bool deposit_jz, bool deposit_rho, bool deposit_chi, + amrex::Geometry const& gm, int const lev) { for (int i=0; i const& a_arr = use_laser ? diff --git a/src/salame/Salame.cpp b/src/salame/Salame.cpp index 9f0d4e8262..f4362bb4b5 100644 --- a/src/salame/Salame.cpp +++ b/src/salame/Salame.cpp @@ -30,14 +30,14 @@ SalameModule (Hipace* hipace, const int n_iter, const bool do_advance, int& last // STEP 1: Calculate what Ez would be with the initial SALAME beam weight // advance plasma to the temp slice, only shift once - hipace->m_multi_plasma.AdvanceParticles(hipace->m_fields, hipace->m_laser, hipace->Geom(lev), + hipace->m_multi_plasma.AdvanceParticles(hipace->m_fields, hipace->m_multi_laser, hipace->Geom(lev), true, true, true, iter==0, lev); hipace->m_fields.duplicate(lev, WhichSlice::Salame, {"jx", "jy"}, WhichSlice::Next, {"jx_beam", "jy_beam"}); // deposit plasma jx and jy on the next temp slice, to the SALANE slice - hipace->m_multi_plasma.DepositCurrent(hipace->m_fields, hipace->m_laser, + hipace->m_multi_plasma.DepositCurrent(hipace->m_fields, hipace->m_multi_laser, WhichSlice::Salame, true, true, false, false, false, hipace->Geom(lev), lev); // use an initial guess of zero for Bx and By in MG solver to reduce relative error @@ -64,7 +64,7 @@ SalameModule (Hipace* hipace, const int n_iter, const bool do_advance, int& last if (do_advance) { SalameOnlyAdvancePlasma(hipace, lev); - hipace->m_multi_plasma.DepositCurrent(hipace->m_fields, hipace->m_laser, + hipace->m_multi_plasma.DepositCurrent(hipace->m_fields, hipace->m_multi_laser, WhichSlice::Salame, true, true, false, false, false, hipace->Geom(lev), lev); } else { SalameGetJxJyFromBxBy(hipace, lev); @@ -110,7 +110,7 @@ SalameModule (Hipace* hipace, const int n_iter, const bool do_advance, int& last hipace->InitializeSxSyWithBeam(lev); - hipace->m_multi_plasma.ExplicitDeposition(hipace->m_fields, hipace->m_laser, + hipace->m_multi_plasma.ExplicitDeposition(hipace->m_fields, hipace->m_multi_laser, hipace->Geom(lev), lev); hipace->ExplicitMGSolveBxBy(lev, WhichSlice::This); diff --git a/tests/laser_blowout_wake_explicit.1Rank.sh b/tests/laser_blowout_wake_explicit.1Rank.sh index ea20714db8..495dc4719a 100755 --- a/tests/laser_blowout_wake_explicit.1Rank.sh +++ b/tests/laser_blowout_wake_explicit.1Rank.sh @@ -29,11 +29,11 @@ mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_normalized \ beams.names = no_beam \ geometry.prob_lo = -20. -20. -8 \ geometry.prob_hi = 20. 20. 6 \ - laser.use_laser = 1 \ + lasers.names = laser \ + lasers.lambda0 = .8e-6 \ laser.a0 = 4.5 \ - laser.lambda0 = .8e-6 \ laser.position_mean = 0. 0. 0 \ - laser.w0 = 4 4 \ + laser.w0 = 4 \ laser.L0 = 2 \ amr.n_cell = 128 128 100 \ diff --git a/tests/laser_blowout_wake_explicit.SI.1Rank.sh b/tests/laser_blowout_wake_explicit.SI.1Rank.sh index 485ac58b8c..f7211abde1 100755 --- a/tests/laser_blowout_wake_explicit.SI.1Rank.sh +++ b/tests/laser_blowout_wake_explicit.SI.1Rank.sh @@ -28,11 +28,11 @@ mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_SI \ beams.names = no_beam \ geometry.prob_lo = -20.*kp_inv -20.*kp_inv -8.*kp_inv \ geometry.prob_hi = 20.*kp_inv 20.*kp_inv 6.*kp_inv \ - laser.use_laser = 1 \ + lasers.names = laser \ + lasers.lambda0 = .8e-6 \ laser.a0 = 4.5 \ - laser.lambda0 = .8e-6 \ laser.position_mean = 0. 0. 0 \ - laser.w0 = 4.*kp_inv 4.*kp_inv \ + laser.w0 = 4.*kp_inv \ laser.L0 = 2.*kp_inv \ amr.n_cell = 128 128 100 \ diff --git a/tests/laser_evolution.SI.2Rank.sh b/tests/laser_evolution.SI.2Rank.sh index 4dec0b8734..7c816e7df0 100755 --- a/tests/laser_evolution.SI.2Rank.sh +++ b/tests/laser_evolution.SI.2Rank.sh @@ -27,7 +27,7 @@ TEST_NAME="${FILE_NAME%.*}" # Run the simulation with multigrid Poisson solver mpiexec -n 2 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_SI \ - laser.solver_type = multigrid \ + lasers.solver_type = multigrid \ hipace.file_prefix = $TEST_NAME # Compare the result with theory $HIPACE_EXAMPLE_DIR/analysis_laser_vacuum.py --output-dir=$TEST_NAME @@ -36,7 +36,7 @@ rm -rf $TEST_NAME # Run the simulation with FFT Poisson solver mpiexec -n 2 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_SI \ - laser.solver_type = fft \ + lasers.solver_type = fft \ hipace.file_prefix = $TEST_NAME # Compare the result with theory $HIPACE_EXAMPLE_DIR/analysis_laser_vacuum.py --output-dir=$TEST_NAME