diff --git a/docs/source/run/parameters.rst b/docs/source/run/parameters.rst index 167da805b0..9bd1d99ce2 100644 --- a/docs/source/run/parameters.rst +++ b/docs/source/run/parameters.rst @@ -159,10 +159,10 @@ Geometry Maximum level of mesh refinement. Currently, mesh refinement is supported up to level `2`. Note, that the mesh refinement algorithm is still in active development and should be used with care. -* ``geometry.patch_lo`` (3 `float`) +* ``geometry.prob_lo`` (3 `float`) Lower end of the simulation box in x, y and z. -* ``geometry.patch_hi`` (3 `float`) +* ``geometry.prob_hi`` (3 `float`) Higher end of the simulation box in x, y and z. * ``geometry.is_periodic`` (3 `bool`) @@ -188,6 +188,16 @@ Geometry * ``mr_lev2.patch_hi`` (3 `float`) Upper end of the refined grid in x, y and z. +* ``lasers.n_cell`` (2 `integer`) + Number of cells in x and y for the laser grid. + The number of cells in the zeta direction is calculated from ``patch_lo`` and ``patch_hi``. + +* ``lasers.patch_lo`` (3 `float`) + Lower end of the laser grid in x, y and z. + +* ``lasers.patch_hi`` (3 `float`) + Upper end of the laser grid in x, y and z. + Time step --------- @@ -805,6 +815,10 @@ Parameters starting with ``lasers.`` apply to all laser pulses, parameters start * ``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. +* ``lasers.interp_order`` (`int`) optional (default `1`) + Transverse shape order for the laser to field interpolation of aabs and + the field to laser interpolation of chi. Currently, `0,1,2,3` are implemented. + * ``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. @@ -898,11 +912,15 @@ Field diagnostics The names of all field diagnostics, separated by a space. Multiple diagnostics can be used to limit the output to only a few relevant regions to save on file size. To run without field diagnostics, choose the name ``no_field_diag``. - If mesh refinement is used, the default becomes ``lev0 lev1`` or ``lev0 lev1 lev2``. + Depending on whether mesh refinement or a laser is used, the default becomes + a subset of ``lev0 lev1 lev2 laser_diag``. -* `` or diagnostic.level`` (`integer`) optional (default `0`) - From which mesh refinement level the diagnostics should be collected. - If ```` is equal to ``lev1``, the default for this parameter becomes 1 etc. +* `` or diagnostic.base_geometry`` (`string`) optional (default `level_0`) + Which geometry the diagnostics should be based on. + Available geometries are `level_0`, `level_1`, `level_2` and `laser`, + depending on if MR or a laser is used. + If ```` is equal to ``lev0 lev1 lev2 laser_diag``, the default for this parameter + becomes ``level_0 level_1 level_2 laser``respectively. * ``.output_period`` (`integer`) optional (default `0`) Output period for fields. No output is given for ``.output_period = 0``. @@ -935,8 +953,12 @@ Field diagnostics If ``rho`` is explicitly mentioned as ``field_data``, it is deposited by the plasma to be available as a diagnostic. Similarly if ``rho_`` is explicitly mentioned, the charge density of that plasma species will be separately available as a diagnostic. - When a laser pulse is used, the laser complex envelope ``laserEnvelope`` is available. + When a laser pulse is used, the laser complex envelope ``laserEnvelope`` is available + in the ``laser`` base geometry. The plasma proper density (n/gamma) is then also accessible via ``chi``. + A field can be removed from the list, for example, after it has been included through ``all``, + by adding ``remove_`` after it has been added. If a field is added and removed + multiple times, the last occurrence takes precedence. * `` or diagnostic.patch_lo`` (3 `float`) optional (default `-infinity -infinity -infinity`) Lower limit for the diagnostic grid. diff --git a/src/Hipace.H b/src/Hipace.H index 03e60a69ef..a11cf6d88c 100644 --- a/src/Hipace.H +++ b/src/Hipace.H @@ -282,7 +282,7 @@ private: amrex::Vector< CoulombCollision > m_all_collisions; void InitDiagnostics (const int step); - void FillFieldDiagnostics (const int lev, int islice); + void FillFieldDiagnostics (const int current_N_level, int islice); void FillBeamDiagnostics (const int step); void WriteDiagnostics (const int step); void FlushDiagnostics (); diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 5a2ce8e164..b25b7d44e3 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -75,7 +75,7 @@ Hipace::Hipace () : m_multi_plasma(), m_adaptive_time_step(m_multi_beam.get_nbeams()), m_multi_laser(), - m_diags(m_N_level) + m_diags(m_N_level, m_multi_laser.UseLaser()) { amrex::ParmParse pp;// Traditionally, max_step and stop_time do not have prefix. queryWithParser(pp, "max_step", m_max_step); @@ -162,7 +162,10 @@ Hipace::Hipace () : MakeGeometry(); - m_use_laser = m_multi_laser.m_use_laser; + // use level 0 as default for laser geometry + m_multi_laser.MakeLaserGeometry(m_3D_geom[0]); + + m_use_laser = m_multi_laser.UseLaser(); queryWithParser(pph, "collisions", m_collision_names); /** Initialize the collision objects */ @@ -196,15 +199,14 @@ Hipace::InitData () amrex::Print() << "using the leapfrog plasma particle pusher\n"; #endif + m_multi_laser.InitData(); + for (int lev=0; lev= m_3D_geom[0].Domain().smallEnd(2)) { m_multi_buffer.get_data(islice-1, m_multi_beam, m_multi_laser, WhichBeamSlice::Next); @@ -576,12 +577,10 @@ Hipace::SolveOneSlice (int islice, int step) m_fields.InSituComputeDiags(step, m_physical_time, islice, m_3D_geom[0], m_max_step, m_max_time); // get laser insitu diagnostics - m_multi_laser.InSituComputeDiags(step, m_physical_time, islice, m_3D_geom[0], m_max_step, m_max_time); + m_multi_laser.InSituComputeDiags(step, m_physical_time, islice, m_max_step, m_max_time); // copy fields (and laser) to diagnostic array - for (int lev=0; lev m_comps_output; /**< Component names to Write to output file */ /** Component indexes to Write to output file */ amrex::Gpu::DeviceVector m_comps_output_idx; /** Vector over levels, all fields */ amrex::FArrayBox m_F; - using complex_type = amrex::GpuComplex; /** FAB for laser */ - amrex::BaseFab m_F_laser; + amrex::BaseFab> m_F_laser; amrex::Geometry m_geom_io; /**< Diagnostics geometry */ bool m_has_field; /**< if there is field output to write */ /** Number of iterations between consecutive output dumps. * Default is 0, meaning no output */ int m_output_period = 0; - /** Name of the laser in input and output files */ - std::string m_laser_io_name = "laserEnvelope"; }; @@ -57,10 +57,10 @@ class Diagnostic public: /** \brief Constructor */ - explicit Diagnostic (int nlev); + explicit Diagnostic (int nlev, bool use_laser); /** \brief Determine which data to output */ - void Initialize (const int lev, bool do_laser); + void Initialize (int nlev, bool use_laser); /** \brief return names of the beams to output */ amrex::Vector& getBeamNames () { return m_output_beam_names; } @@ -80,7 +80,7 @@ public: int output_step, int max_step, amrex::Real output_time, amrex::Real max_time) { - return (fd.m_nfields > 0 || fd.m_do_laser) && + return (fd.m_nfields > 0) && utils::doDiagnostics(fd.m_output_period, output_step, max_step, output_time, max_time); } @@ -148,16 +148,15 @@ public: /** \brief resizes the FArrayBox of the diagnostics to the currently calculated box * - * \param[in] a_domain box to which the Geometry of the diagnostics will be resized to - * \param[in] lev MR level - * \param[in] geom geometry of the full simulation domain + * \param[in] field_geom geometry of the full simulation domain for all field MR levels + * \param[in] laser_geom geometry of the full simulation domain for the laser * \param[in] output_step current step index * \param[in] max_step last step index * \param[in] output_time physical time of current step * \param[in] max_time physical time of last step */ - void ResizeFDiagFAB (const amrex::Box a_domain, const int lev, - amrex::Geometry const& geom, int output_step, int max_step, + void ResizeFDiagFAB (amrex::Vector& field_geom, + amrex::Geometry const& laser_geom, int output_step, int max_step, amrex::Real output_time, amrex::Real max_time); private: diff --git a/src/diagnostics/Diagnostic.cpp b/src/diagnostics/Diagnostic.cpp index bb4a90bd7f..963bb1b381 100644 --- a/src/diagnostics/Diagnostic.cpp +++ b/src/diagnostics/Diagnostic.cpp @@ -10,15 +10,26 @@ #include "utils/HipaceProfilerWrapper.H" #include -Diagnostic::Diagnostic (int nlev) +#include +#include +#include +#include +#include + +Diagnostic::Diagnostic (int nlev, bool use_laser) { amrex::ParmParse ppd("diagnostic"); amrex::ParmParse pph("hipace"); + // Make the default diagnostic objects, subset of: lev0, lev1, lev2, laser_diag amrex::Vector field_diag_names{}; - - for(int ilev = 0; ilev diag_name_to_default_geometry{}; + // for each geometry name, is it based on fields or laser + std::map geometry_name_to_geom_type{}; + // for each geometry name, if its for fields what MR level is it on + std::map geometry_name_to_level{}; + // for each geometry, what output components are available + std::map> geometry_name_to_output_comps{}; + // in case there is an error, generate a string with all available geometries and components + std::stringstream all_comps_error_str{}; + for (int lev = 0; lev'.\n"; + + // keep track of all components from the input and later assert that they were all used + std::map is_global_comp_used{}; + for (auto& fd : m_field_data) { amrex::ParmParse pp(fd.m_diag_name); - queryWithParserAlt(pp, "field_data", fd.m_comps_output, ppd); + std::string base_geom_name = "level_0"; - amrex::Vector all_field_comps{}; - all_field_comps.reserve(Comps[WhichSlice::This].size() + (do_laser ? 1 : 0)); - for (const auto& [comp, idx] : Comps[WhichSlice::This]) { - all_field_comps.push_back(comp); - } - if (do_laser) { - all_field_comps.push_back(fd.m_laser_io_name); + if (diag_name_to_default_geometry.count(fd.m_diag_name) > 0) { + base_geom_name = diag_name_to_default_geometry.at(fd.m_diag_name); } - if(fd.m_comps_output.empty()) { - fd.m_comps_output = all_field_comps; + queryWithParserAlt(pp, "base_geometry", base_geom_name, ppd); + + if (geometry_name_to_geom_type.count(base_geom_name) > 0) { + fd.m_base_geom_type = geometry_name_to_geom_type.at(base_geom_name); + fd.m_level = geometry_name_to_level.at(base_geom_name); } else { - for (const std::string& comp_name : fd.m_comps_output) { - if (comp_name == "all" || comp_name == "All") { - fd.m_comps_output = all_field_comps; - break; - } - if (comp_name == "none" || comp_name == "None") { - fd.m_comps_output.clear(); - break; - } - if (std::find(all_field_comps.begin(), all_field_comps.end(), comp_name) == all_field_comps.end()) { - std::stringstream error_str{}; - error_str << "Unknown field diagnostics component: " << comp_name <<"\nmust be " - << "'all', 'none' or a subset of:"; - for (auto& comp : all_field_comps) { - error_str << " " << comp; - } - amrex::Abort(error_str.str()); - } - } + amrex::Abort("Unknown diagnostics base_geometry: '" + base_geom_name + "'!\n" + + all_comps_error_str.str()); } - // remove laser from fd.m_comps_output because it is output separately - for (auto it = fd.m_comps_output.begin(); it != fd.m_comps_output.end();) { - if (*it == fd.m_laser_io_name) { - it = fd.m_comps_output.erase(it); - fd.m_do_laser = true; + amrex::Vector use_comps{}; + const bool use_local_comps = queryWithParser(pp, "field_data", use_comps); + if (!use_local_comps) { + queryWithParser(ppd, "field_data", use_comps); + } + + // set to store all used components to avoid duplicates + std::set comps_set{}; + + if (use_comps.empty()) { + // by default output all components + use_comps.push_back("all"); + } + + // iterate through the user-provided components from left to right + for (const std::string& comp_name : use_comps) { + if (comp_name == "all" || comp_name == "All") { + is_global_comp_used[comp_name] = true; + // insert all available components + comps_set.insert(geometry_name_to_output_comps[base_geom_name].begin(), + geometry_name_to_output_comps[base_geom_name].end()); + } else if (comp_name == "none" || comp_name == "None") { + is_global_comp_used[comp_name] = true; + // remove all components + comps_set.clear(); + } else if (geometry_name_to_output_comps[base_geom_name].count(comp_name) > 0) { + is_global_comp_used[comp_name] = true; + // insert requested component + comps_set.insert(comp_name); + } else if (comp_name.find("remove_") == 0 && + geometry_name_to_output_comps[base_geom_name].count( + comp_name.substr(std::string("remove_").size(), comp_name.size())) > 0) { + is_global_comp_used[comp_name] = true; + // remove requested component + comps_set.erase(comp_name.substr(std::string("remove_").size(), comp_name.size())); + } else if (use_local_comps) { + // if field_data was specified through , + // assert that all components exist in the geometry + amrex::Abort("Unknown diagnostics field_data '" + comp_name + + "' in base_geometry '" + base_geom_name + "'!\n" + + all_comps_error_str.str()); } else { - ++it; + // if field_data was specified through diagnostic, + // check later that all components are at least used by one of the diagnostics + is_global_comp_used.try_emplace(comp_name, false); } } + fd.m_comps_output.assign(comps_set.begin(), comps_set.end()); fd.m_nfields = fd.m_comps_output.size(); - amrex::Gpu::PinnedVector local_comps_output_idx(fd.m_nfields); - for(int i = 0; i < fd.m_nfields; ++i) { - local_comps_output_idx[i] = Comps[WhichSlice::This][fd.m_comps_output[i]]; + // copy the indexes of m_comps_output to the GPU + if (fd.m_base_geom_type == FieldDiagnosticData::geom_type::field) { + amrex::Gpu::PinnedVector local_comps_output_idx(fd.m_nfields); + for(int i = 0; i < fd.m_nfields; ++i) { + local_comps_output_idx[i] = Comps[WhichSlice::This][fd.m_comps_output[i]]; + } + fd.m_comps_output_idx.assign(local_comps_output_idx.begin(), local_comps_output_idx.end()); } - fd.m_comps_output_idx.assign(local_comps_output_idx.begin(), local_comps_output_idx.end()); + } + + // check that all components are at least used by one of the diagnostics + for (auto& [key, val] : is_global_comp_used) { + if (!val) { + amrex::Abort("Unknown or unused component in diagnostic.field_data.\n'" + + key + "' does not belong to any diagnostic.names!\n" + + all_comps_error_str.str()); + } + } - if (m_field_data.size() != 1) { + // if there are multiple diagnostic objects with the same m_base_geom_type (colliding component + // names), append the name of the diagnostic object to the component name in the output + for (auto& fd : m_field_data) { + if (1 < std::count_if(m_field_data.begin(), m_field_data.end(), [&] (auto& fd2) { + return fd.m_base_geom_type == fd2.m_base_geom_type; + })) { for (auto& comp_name : fd.m_comps_output) { comp_name += "_" + fd.m_diag_name; } - if (fd.m_do_laser) { - fd.m_laser_io_name += "_" + fd.m_diag_name; - } } } @@ -219,17 +296,27 @@ Diagnostic::Initialize (const int lev, bool do_laser) { } void -Diagnostic::ResizeFDiagFAB (const amrex::Box a_domain, const int lev, - amrex::Geometry const& geom, int output_step, int max_step, +Diagnostic::ResizeFDiagFAB (amrex::Vector& field_geom, + amrex::Geometry const& laser_geom, int output_step, int max_step, amrex::Real output_time, amrex::Real max_time) { AMREX_ALWAYS_ASSERT(m_initialized); for (auto& fd : m_field_data) { - if (fd.m_level != lev) continue; + amrex::Geometry geom; - amrex::Box domain = a_domain; + // choose the geometry of the diagnostic + switch (fd.m_base_geom_type) { + case FieldDiagnosticData::geom_type::field: + geom = field_geom[fd.m_level]; + break; + case FieldDiagnosticData::geom_type::laser: + geom = laser_geom; + break; + } + + amrex::Box domain = geom.Domain(); if (fd.m_include_ghost_cells) { domain.grow(Fields::m_slices_nguards); @@ -279,12 +366,15 @@ Diagnostic::ResizeFDiagFAB (const amrex::Box a_domain, const int lev, if(fd.m_has_field) { HIPACE_PROFILE("Diagnostic::ResizeFDiagFAB()"); - fd.m_F.resize(domain, fd.m_nfields, amrex::The_Pinned_Arena()); - fd.m_F.setVal(0); - - if (fd.m_do_laser) { - fd.m_F_laser.resize(domain, 1, amrex::The_Pinned_Arena()); - fd.m_F_laser.setVal({0,0}); + switch (fd.m_base_geom_type) { + case FieldDiagnosticData::geom_type::field: + fd.m_F.resize(domain, fd.m_nfields, amrex::The_Pinned_Arena()); + fd.m_F.setVal(0); + break; + case FieldDiagnosticData::geom_type::laser: + fd.m_F_laser.resize(domain, fd.m_nfields, amrex::The_Pinned_Arena()); + fd.m_F_laser.setVal({0,0}); + break; } } } diff --git a/src/diagnostics/OpenPMDWriter.cpp b/src/diagnostics/OpenPMDWriter.cpp index 4f91f039d8..c66ce30933 100644 --- a/src/diagnostics/OpenPMDWriter.cpp +++ b/src/diagnostics/OpenPMDWriter.cpp @@ -90,21 +90,14 @@ OpenPMDWriter::WriteFieldData ( // todo: periodicity/boundary, field solver, particle pusher, etc. auto meshes = iteration.meshes; - amrex::Vector varnames = fd.m_comps_output; - - if (fd.m_do_laser) { - // laser must be at the end of varnames - varnames.push_back(fd.m_laser_io_name); - } - // loop over field components - for ( int icomp = 0; icomp < varnames.size(); ++icomp ) + for ( int icomp = 0; icomp < fd.m_nfields; ++icomp ) { - const bool is_laser_comp = varnames[icomp].find("laserEnvelope") == 0; + const bool is_laser_comp = fd.m_base_geom_type == FieldDiagnosticData::geom_type::laser; // "B" "x" (todo) // "Bx" "" (just for now) - openPMD::Mesh field = meshes[varnames[icomp]]; + openPMD::Mesh field = meshes[fd.m_comps_output[icomp]]; openPMD::MeshRecordComponent field_comp = field[openPMD::MeshRecordComponent::SCALAR]; // meta-data diff --git a/src/fields/Fields.H b/src/fields/Fields.H index 96e6b485a4..8095ffbac8 100644 --- a/src/fields/Fields.H +++ b/src/fields/Fields.H @@ -133,14 +133,14 @@ public: } /** \brief Copy between the full FArrayBox and slice MultiFab. * - * \param[in] lev MR level + * \param[in] current_N_level number of MR levels active on the current slice * \param[in] i_slice z slice in which to write the data * \param[in,out] fd data for field diagnostics - * \param[in] calc_geom main geometry + * \param[in] field_geom main field geometry * \param[in] multi_laser MultiLaser object */ - void Copy (const int lev, const int i_slice, FieldDiagnosticData& fd, - const amrex::Geometry& calc_geom, MultiLaser& multi_laser); + void Copy (const int current_N_level, const int i_slice, FieldDiagnosticData& fd, + const amrex::Vector& field_geom, MultiLaser& multi_laser); /** \brief Initialize all required fields to zero and interpolate from lev-1 to lev if needed * diff --git a/src/fields/Fields.cpp b/src/fields/Fields.cpp index efcf9d8db3..5f6654cdbc 100644 --- a/src/fields/Fields.cpp +++ b/src/fields/Fields.cpp @@ -437,15 +437,15 @@ Multiply (const amrex::IntVect box_grow, amrex::MultiFab dst, } void -Fields::Copy (const int lev, const int i_slice, FieldDiagnosticData& fd, - const amrex::Geometry& calc_geom, MultiLaser& multi_laser) +Fields::Copy (const int current_N_level, const int i_slice, FieldDiagnosticData& fd, + const amrex::Vector& field_geom, MultiLaser& multi_laser) { HIPACE_PROFILE("Fields::Copy()"); constexpr int depos_order_xy = 1; constexpr int depos_order_z = 1; constexpr int depos_order_offset = depos_order_z / 2 + 1; - const amrex::Real poff_calc_z = GetPosOffset(2, calc_geom, calc_geom.Domain()); + const amrex::Real poff_calc_z = GetPosOffset(2, field_geom[0], field_geom[0].Domain()); const amrex::Real poff_diag_x = GetPosOffset(0, fd.m_geom_io, fd.m_geom_io.Domain()); const amrex::Real poff_diag_y = GetPosOffset(1, fd.m_geom_io, fd.m_geom_io.Domain()); const amrex::Real poff_diag_z = GetPosOffset(2, fd.m_geom_io, fd.m_geom_io.Domain()); @@ -454,21 +454,21 @@ Fields::Copy (const int lev, const int i_slice, FieldDiagnosticData& fd, // Calculate to which diag_fab slices this slice could contribute const int i_slice_min = i_slice - depos_order_offset; const int i_slice_max = i_slice + depos_order_offset; - const amrex::Real pos_slice_min = i_slice_min * calc_geom.CellSize(2) + poff_calc_z; - const amrex::Real pos_slice_max = i_slice_max * calc_geom.CellSize(2) + poff_calc_z; + const amrex::Real pos_slice_min = i_slice_min * field_geom[0].CellSize(2) + poff_calc_z; + const amrex::Real pos_slice_max = i_slice_max * field_geom[0].CellSize(2) + poff_calc_z; int k_min = static_cast(amrex::Math::round((pos_slice_min - poff_diag_z) * fd.m_geom_io.InvCellSize(2))); const int k_max = static_cast(amrex::Math::round((pos_slice_max - poff_diag_z) * fd.m_geom_io.InvCellSize(2))); - amrex::Box diag_box = fd.m_F.box(); + amrex::Box diag_box = fd.m_geom_io.Domain(); if (fd.m_slice_dir != 2) { // Put contributions from i_slice to different diag_fab slices in GPU vector m_rel_z_vec.resize(k_max+1-k_min); m_rel_z_vec_cpu.resize(k_max+1-k_min); for (int k=k_min; k<=k_max; ++k) { const amrex::Real pos = k * fd.m_geom_io.CellSize(2) + poff_diag_z; - const amrex::Real mid_i_slice = (pos - poff_calc_z)*calc_geom.InvCellSize(2); + const amrex::Real mid_i_slice = (pos - poff_calc_z)*field_geom[0].InvCellSize(2); amrex::Real sz_cell[depos_order_z + 1]; const int k_cell = compute_shape_factor(sz_cell, mid_i_slice); m_rel_z_vec_cpu[k-k_min] = 0; @@ -495,17 +495,18 @@ Fields::Copy (const int lev, const int i_slice, FieldDiagnosticData& fd, } else { m_rel_z_vec.resize(1); m_rel_z_vec_cpu.resize(1); - const amrex::Real pos_z = i_slice * calc_geom.CellSize(2) + poff_calc_z; + const amrex::Real pos_z = i_slice * field_geom[0].CellSize(2) + poff_calc_z; if (fd.m_geom_io.ProbLo(2) <= pos_z && pos_z <= fd.m_geom_io.ProbHi(2)) { - m_rel_z_vec_cpu[0] = calc_geom.CellSize(2); + m_rel_z_vec_cpu[0] = field_geom[0].CellSize(2); k_min = 0; } else { return; } } if (diag_box.isEmpty()) return; - auto& slice_mf = m_slices[lev]; - auto slice_func = interpolated_field_xy{{slice_mf}, calc_geom}; + auto& slice_mf = m_slices[fd.m_level]; + auto slice_func = interpolated_field_xy{{slice_mf}, field_geom[fd.m_level]}; auto& laser_mf = multi_laser.getSlices(); auto laser_func = interpolated_field_xy{{laser_mf}, multi_laser.GetLaserGeom()}; @@ -527,7 +528,8 @@ Fields::Copy (const int lev, const int i_slice, FieldDiagnosticData& fd, const amrex::Real dx = fd.m_geom_io.CellSize(0); const amrex::Real dy = fd.m_geom_io.CellSize(1); - if (fd.m_nfields > 0) { + if (fd.m_base_geom_type == FieldDiagnosticData::geom_type::field && + current_N_level > fd.m_level) { auto slice_array = slice_func.array(mfi); amrex::Array4 diag_array = fd.m_F.array(); amrex::ParallelFor(diag_box, fd.m_nfields, @@ -538,17 +540,16 @@ Fields::Copy (const int lev, const int i_slice, FieldDiagnosticData& fd, const int m = n[diag_comps]; diag_array(i,j,k,n) += rel_z_data[k-k_min] * slice_array(x,y,m); }); - } - - if (fd.m_do_laser) { + } else if (fd.m_base_geom_type == FieldDiagnosticData::geom_type::laser && + multi_laser.UseLaser(i_slice)) { auto laser_array = laser_func.array(mfi); - amrex::Array4 diag_array_laser = fd.m_F_laser.array(); + amrex::Array4> diag_array_laser = fd.m_F_laser.array(); amrex::ParallelFor(diag_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const amrex::Real x = i * dx + poff_diag_x; const amrex::Real y = j * dy + poff_diag_y; - diag_array_laser(i,j,k) += FieldDiagnosticData::complex_type { + diag_array_laser(i,j,k) += amrex::GpuComplex { rel_z_data[k-k_min] * laser_array(x,y,WhichLaserSlice::n00j00_r), rel_z_data[k-k_min] * laser_array(x,y,WhichLaserSlice::n00j00_i) }; diff --git a/src/laser/MultiLaser.H b/src/laser/MultiLaser.H index e8d516ec1b..4dc44f7f1e 100644 --- a/src/laser/MultiLaser.H +++ b/src/laser/MultiLaser.H @@ -11,7 +11,6 @@ #define MULTILASER_H_ #include "Laser.H" -#include "fields/Fields.H" #include "mg_solver/HpMultiGrid.H" #include "fields/fft_poisson_solver/fft/AnyFFT.H" @@ -43,12 +42,16 @@ namespace WhichLaserSlice { np1jp1_i, np1jp2_r, np1jp2_i, + chi, + chi_initial, N }; } class Fields; +class MultiPlasma; + class MultiLaser { @@ -70,14 +73,13 @@ public: /** get function for the 2D slices (const version) */ const amrex::MultiFab& getSlices () const {return m_slices; } - /** \brief Allocate laser multifab - * \param[in] slice_ba box array of the slice - * \param[in] slice_dm corresponding distribution mapping - * \param[in] geom_3D 3D Geometry for level 0 + /** \brief Make Laser geometry + * \param[in] field_geom_3D 3D Geometry for level 0 */ - void InitData (const amrex::BoxArray& slice_ba, - const amrex::DistributionMapping& slice_dm, - const amrex::Geometry& geom_3D); + void MakeLaserGeometry (const amrex::Geometry& field_geom_3D); + + /** \brief Allocate laser multifab */ + void InitData (); /** \brief Initialize on slice of the 3D laser field. * @@ -86,85 +88,89 @@ public: */ void InitSliceEnvelope (const int islice, const int comp); - /** \brief Read in a laser from an openPMD file - * - * \param[in] gm Geometry for level 0 - */ - void GetEnvelopeFromFileHelper (const amrex::Geometry& gm); + /** \brief Read in a laser from an openPMD file */ + void GetEnvelopeFromFileHelper (); - /** \brief Read in a laser from an openPMD file - * - * \param[in] gm Geometry for level 0 - */ + /** \brief Read in a laser from an openPMD file */ template - void GetEnvelopeFromFile (const amrex::Geometry& gm); + void GetEnvelopeFromFile (); /** \brief Shift 2D slices in zeta + * \param[in] islice slice index */ - void ShiftLaserSlices (); + void ShiftLaserSlices (const int islice); /** Write Aabs into Fields MultiFab + * \param[in] islice slice index * \param[in] current_N_level number of MR levels active on the current slice * \param[in] fields Field object - * \param[in] geom Geometry of the problem + * \param[in] field_geom Geometry of the problem + */ + void UpdateLaserAabs (const int islice, const int current_N_level, Fields& fields, + amrex::Vector const& field_geom); + + /** Put Chi from the fields and initial chi into the chi component of the laser + * \param[in] fields Field object + * \param[in] geom_field_lev0 Geometry of the fields on MR level 0 */ - void UpdateLaserAabs (const int current_N_level, Fields& fields, - amrex::Vector const& geom); + void InterpolateChi (const Fields& fields, amrex::Geometry const& geom_field_lev0); + + /** Fill the chi_initial component of the laser using the Plasma density function + * \param[in] multi_plasma multi plasma to get the density function, charge and mass + */ + void SetInitialChi (const MultiPlasma& multi_plasma); /** Wrapper function to advance a laser slice by 1 time step. + * \param[in] islice slice index * \param[in] fields Field object * \param[in] dt time step of the simulation * \param[in] step current iteration. Needed because step 0 needs a specific treatment. + * \param[in] geom_field_lev0 Geometry of the fields on MR level 0 */ - void AdvanceSlice (const Fields& fields, amrex::Real dt, int step); + void AdvanceSlice (const int islice, const Fields& fields, amrex::Real dt, int step, + amrex::Geometry const& geom_field_lev0); /** 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] dt time step of the simulation * \param[in] step current iteration. Needed because step 0 needs a specific treatment. */ - void AdvanceSliceMG (const Fields& fields, amrex::Real dt, int step); + void AdvanceSliceMG (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] dt time step of the simulation * \param[in] step current iteration. Needed because step 0 needs a specific treatment. */ - void AdvanceSliceFFT (const Fields& fields, amrex::Real dt, int step); + void AdvanceSliceFFT (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 * \param[in] comp laser component to initialize */ - void InitLaserSlice (const amrex::Geometry& geom, const int islice, const int comp); + void InitLaserSlice (const int islice, const int comp); /** Compute in-situ laser diagnostics of current slice, store in member variable * \param[in] step current time step * \param[in] time physical time * \param[in] islice current slice, on which diags are computed. - * \param[in] geom3D Geometry of the problem * \param[in] max_step maximum time step of simulation * \param[in] max_time maximum time of simulation */ - void InSituComputeDiags (int step, amrex::Real time, int islice, const amrex::Geometry& geom3D, + void InSituComputeDiags (int step, amrex::Real time, int islice, int max_step, amrex::Real max_time); /** Dump in-situ reduced diagnostics to file. * \param[in] step current time step * \param[in] time physical time - * \param[in] geom3D Geometry object for the whole domain * \param[in] max_step maximum time step of simulation * \param[in] max_time maximum time of simulation */ - void InSituWriteToFile (int step, amrex::Real time, const amrex::Geometry& geom3D, - int max_step, amrex::Real max_time); + void InSituWriteToFile (int step, amrex::Real time, int max_step, amrex::Real max_time); /** Get the central wavelength */ amrex::Real GetLambda0 () const { return m_lambda0; } @@ -172,10 +178,25 @@ public: /** Get the geometry of the Laser Box */ const amrex::Geometry& GetLaserGeom () const { return m_laser_geom_3D; } - bool m_use_laser {false}; /**< whether a laser is used or not */ + /** If the laser geometry includes this slice + * \param[in] islice slice index + */ + bool HasSlice (const int islice) const { + return GetLaserGeom().Domain().smallEnd(2) <= islice && + islice <= GetLaserGeom().Domain().bigEnd(2); + } + + /** If the laser is used */ + bool UseLaser () const { return m_use_laser; } + + /** If the laser is used and the laser geometry includes this slice + * \param[in] islice slice index + */ + bool UseLaser (const int islice) const { return m_use_laser && HasSlice(islice); } private: + bool m_use_laser {false}; /**< whether a laser is used or not */ /** 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. */ @@ -187,8 +208,16 @@ private: amrex::IntVect m_slices_nguards = {-1, -1, -1}; std::string m_solver_type = "multigrid"; bool m_use_phase {true}; - amrex::Box m_slice_box; + /** 3D Laser Geometry */ amrex::Geometry m_laser_geom_3D; + /** xy slice BoxArray. Contains only one box */ + amrex::BoxArray m_laser_slice_ba; + /** xy slice DistributionMapping */ + amrex::DistributionMapping m_laser_slice_dm; + /** slice box of laser */ + amrex::Box m_slice_box; + /** interpolation order for laser to field and field to laser operations */ + int m_interp_order = 1; /** if the lasers are initialized from openPMD file */ bool m_laser_from_file = false; diff --git a/src/laser/MultiLaser.cpp b/src/laser/MultiLaser.cpp index ba77cbcdc5..ac84172959 100644 --- a/src/laser/MultiLaser.cpp +++ b/src/laser/MultiLaser.cpp @@ -9,7 +9,10 @@ #include "MultiLaser.H" #include "utils/Constants.H" +#include "fields/Fields.H" #include "Hipace.H" +#include "particles/plasma/MultiPlasma.H" +#include "particles/particles_utils/ShapeFactors.H" #include "utils/HipaceProfilerWrapper.H" #include "utils/DeprecatedInput.H" #include "utils/InsituUtil.H" @@ -45,6 +48,8 @@ MultiLaser::ReadParameters () 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"); + queryWithParser(pp, "interp_order", m_interp_order); + AMREX_ALWAYS_ASSERT(m_interp_order <= 3 && m_interp_order >= 0); 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); @@ -67,9 +72,63 @@ MultiLaser::ReadParameters () void -MultiLaser::InitData (const amrex::BoxArray& slice_ba, - const amrex::DistributionMapping& slice_dm, - const amrex::Geometry& geom_3D) +MultiLaser::MakeLaserGeometry (const amrex::Geometry& field_geom_3D) +{ + if (!m_use_laser) return; + amrex::ParmParse pp("lasers"); + + // use field_geom_3D as the default + std::array n_cells_laser {field_geom_3D.Domain().length(0), + field_geom_3D.Domain().length(1)}; + std::array patch_lo_laser { + field_geom_3D.ProbDomain().lo(0), + field_geom_3D.ProbDomain().lo(1), + field_geom_3D.ProbDomain().lo(2)}; + std::array patch_hi_laser { + field_geom_3D.ProbDomain().hi(0), + field_geom_3D.ProbDomain().hi(1), + field_geom_3D.ProbDomain().hi(2)}; + + // get parameters from user input + queryWithParser(pp, "n_cell", n_cells_laser); + queryWithParser(pp, "patch_lo", patch_lo_laser); + queryWithParser(pp, "patch_hi", patch_hi_laser); + + // round zeta lo and hi to full cells + const amrex::Real pos_offset_z = GetPosOffset(2, field_geom_3D, field_geom_3D.Domain()); + + const int zeta_lo = std::max( field_geom_3D.Domain().smallEnd(2), + int(amrex::Math::round((patch_lo_laser[2] - pos_offset_z) * field_geom_3D.InvCellSize(2))) + ); + + const int zeta_hi = std::min( field_geom_3D.Domain().bigEnd(2), + int(amrex::Math::round((patch_hi_laser[2] - pos_offset_z) * field_geom_3D.InvCellSize(2))) + ); + + patch_lo_laser[2] = (zeta_lo-0.5)*field_geom_3D.CellSize(2) + pos_offset_z; + patch_hi_laser[2] = (zeta_hi+0.5)*field_geom_3D.CellSize(2) + pos_offset_z; + + // make the boxes + const amrex::Box domain_3D_laser{amrex::IntVect(0, 0, zeta_lo), + amrex::IntVect(n_cells_laser[0]-1, n_cells_laser[1]-1, zeta_hi)}; + + const amrex::RealBox real_box(patch_lo_laser, patch_hi_laser); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(real_box.volume() > 0., "Laser box must have positive volume"); + + // make the geometry, slice box and ba and dm + m_laser_geom_3D.define(domain_3D_laser, real_box, amrex::CoordSys::cartesian, {0, 0, 0}); + + m_slice_box = domain_3D_laser; + m_slice_box.setSmall(2, 0); + m_slice_box.setBig(2, 0); + + m_laser_slice_ba.define(m_slice_box); + m_laser_slice_dm.define(amrex::Vector({amrex::ParallelDescriptor::MyProc()})); +} + +void +MultiLaser::InitData () { if (!m_use_laser) return; @@ -80,13 +139,10 @@ MultiLaser::InitData (const amrex::BoxArray& slice_ba, int nguards_xy = (Hipace::m_depos_order_xy + 1) / 2 + 1; m_slices_nguards = {nguards_xy, nguards_xy, 0}; m_slices.define( - slice_ba, slice_dm, WhichLaserSlice::N, m_slices_nguards, + m_laser_slice_ba, m_laser_slice_dm, WhichLaserSlice::N, m_slices_nguards, amrex::MFInfo().SetArena(amrex::The_Arena())); m_slices.setVal(0.0); - m_laser_geom_3D = geom_3D; - m_slice_box = slice_ba[0]; - if (m_solver_type == "fft") { m_sol.resize(m_slice_box, 1, amrex::The_Arena()); m_rhs.resize(m_slice_box, 1, amrex::The_Arena()); @@ -108,7 +164,7 @@ MultiLaser::InitData (const amrex::BoxArray& slice_ba, if (m_laser_from_file) { if (Hipace::HeadRank()) { m_F_input_file.resize(m_laser_geom_3D.Domain(), 2, amrex::The_Pinned_Arena()); - GetEnvelopeFromFileHelper(m_laser_geom_3D); + GetEnvelopeFromFileHelper(); } #ifdef AMREX_USE_MPI // need to communicate m_lambda0 as it is read in from the input file only by the head rank @@ -136,7 +192,7 @@ MultiLaser::InitData (const amrex::BoxArray& slice_ba, void MultiLaser::InitSliceEnvelope (const int islice, const int comp) { - if (!m_use_laser) return; + if (!UseLaser(islice)) return; HIPACE_PROFILE("MultiLaser::InitSliceEnvelope()"); @@ -147,13 +203,13 @@ MultiLaser::InitSliceEnvelope (const int islice, const int comp) m_slices[0].copy(m_F_input_file, src_box, 0, m_slice_box, comp, 2); } else { // Compute initial field on the current (device) slice comp and comp + 1 - InitLaserSlice(m_laser_geom_3D, islice, comp); + InitLaserSlice(islice, comp); } } void -MultiLaser::GetEnvelopeFromFileHelper (const amrex::Geometry& gm) { +MultiLaser::GetEnvelopeFromFileHelper () { HIPACE_PROFILE("MultiLaser::GetEnvelopeFromFileHelper()"); #ifdef HIPACE_USE_OPENPMD @@ -195,14 +251,13 @@ MultiLaser::GetEnvelopeFromFileHelper (const amrex::Geometry& gm) { } if (input_type == openPMD::Datatype::CFLOAT) { - GetEnvelopeFromFile>(gm); + GetEnvelopeFromFile>(); } else if (input_type == openPMD::Datatype::CDOUBLE) { - GetEnvelopeFromFile>(gm); + GetEnvelopeFromFile>(); } else { amrex::Abort("Unknown Datatype used in Laser input file. Must use CDOUBLE or CFLOAT\n"); } #else - amrex::ignore_unused(gm); amrex::Abort("loading a laser envelope from an external file requires openPMD support: " "Add HiPACE_OPENPMD=ON when compiling HiPACE++.\n"); #endif // HIPACE_USE_OPENPMD @@ -210,7 +265,7 @@ MultiLaser::GetEnvelopeFromFileHelper (const amrex::Geometry& gm) { template void -MultiLaser::GetEnvelopeFromFile (const amrex::Geometry& gm) { +MultiLaser::GetEnvelopeFromFile () { using namespace amrex::literals; @@ -219,7 +274,7 @@ MultiLaser::GetEnvelopeFromFile (const amrex::Geometry& gm) { const PhysConst phc = get_phys_const(); const amrex::Real clight = phc.c; - const amrex::Box& domain = gm.Domain(); + const amrex::Box& domain = m_laser_geom_3D.Domain(); auto series = openPMD::Series( m_input_file_path , openPMD::Access::READ_ONLY ); auto laser = series.iterations[m_file_num_iteration].meshes[m_file_envelope_name]; @@ -257,13 +312,13 @@ MultiLaser::GetEnvelopeFromFile (const amrex::Geometry& gm) { series.flush(); constexpr int interp_order_xy = 1; - const amrex::Real dx = gm.CellSize(Direction::x); - const amrex::Real dy = gm.CellSize(Direction::y); - const amrex::Real dz = gm.CellSize(Direction::z); - const amrex::Real xmin = gm.ProbLo(Direction::x)+dx/2; - const amrex::Real ymin = gm.ProbLo(Direction::y)+dy/2; - const amrex::Real zmin = gm.ProbLo(Direction::z)+dz/2; - const amrex::Real zmax = gm.ProbHi(Direction::z)-dz/2; + const amrex::Real dx = m_laser_geom_3D.CellSize(Direction::x); + const amrex::Real dy = m_laser_geom_3D.CellSize(Direction::y); + const amrex::Real dz = m_laser_geom_3D.CellSize(Direction::z); + const amrex::Real xmin = m_laser_geom_3D.ProbLo(Direction::x)+dx/2; + const amrex::Real ymin = m_laser_geom_3D.ProbLo(Direction::y)+dy/2; + const amrex::Real zmin = m_laser_geom_3D.ProbLo(Direction::z)+dz/2; + const amrex::Real zmax = m_laser_geom_3D.ProbHi(Direction::z)-dz/2; const int imin = domain.smallEnd(0); const int jmin = domain.smallEnd(1); const int kmin = domain.smallEnd(2); @@ -436,9 +491,9 @@ MultiLaser::GetEnvelopeFromFile (const amrex::Geometry& gm) { } void -MultiLaser::ShiftLaserSlices () +MultiLaser::ShiftLaserSlices (const int islice) { - if (!m_use_laser) return; + if (!UseLaser(islice)) return; HIPACE_PROFILE("MultiLaser::ShiftLaserSlices()"); @@ -470,50 +525,218 @@ MultiLaser::ShiftLaserSlices () } void -MultiLaser::UpdateLaserAabs (const int current_N_level, Fields& fields, - amrex::Vector const& geom) +MultiLaser::UpdateLaserAabs (const int islice, const int current_N_level, Fields& fields, + amrex::Vector const& field_geom) { if (!m_use_laser) return; + if (!HasSlice(islice) && !HasSlice(islice + 1)) return; HIPACE_PROFILE("MultiLaser::UpdateLaserAabs()"); + if (!HasSlice(islice)) { + // set aabs to zero if there is no laser on this slice + // we only need to do this if the previous slice (slice + 1) had a laser + for (int lev=0; lev laser_arr = m_slices.const_array(mfi); const Array2 field_arr = fields.getSlices(0).array(mfi, Comps[WhichSlice::This]["aabs"]); - amrex::ParallelFor(mfi.growntilebox(), - [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept { + const amrex::Real poff_field_x = GetPosOffset(0, field_geom[0], field_geom[0].Domain()); + const amrex::Real poff_field_y = GetPosOffset(1, field_geom[0], field_geom[0].Domain()); + const amrex::Real poff_laser_x = GetPosOffset(0, m_laser_geom_3D, m_laser_geom_3D.Domain()); + const amrex::Real poff_laser_y = GetPosOffset(1, m_laser_geom_3D, m_laser_geom_3D.Domain()); + + const amrex::Real dx_field = field_geom[0].CellSize(0); + const amrex::Real dy_field = field_geom[0].CellSize(1); + const amrex::Real dx_laser_inv = m_laser_geom_3D.InvCellSize(0); + const amrex::Real dy_laser_inv = m_laser_geom_3D.InvCellSize(1); + + const int x_lo = m_slice_box.smallEnd(0); + const int x_hi = m_slice_box.bigEnd(0); + const int y_lo = m_slice_box.smallEnd(1); + const int y_hi = m_slice_box.bigEnd(1); + + amrex::ParallelFor( + amrex::TypeList>{}, + {m_interp_order}, + mfi.growntilebox(), + [=] AMREX_GPU_DEVICE(int i, int j, int, auto interp_order) noexcept { using namespace WhichLaserSlice; - field_arr(i,j) = abssq(laser_arr(i,j,n00j00_r), laser_arr(i,j,n00j00_i)); + const amrex::Real x = i * dx_field + poff_field_x; + const amrex::Real y = j * dy_field + poff_field_y; + + const amrex::Real xmid = (x - poff_laser_x) * dx_laser_inv; + const amrex::Real ymid = (y - poff_laser_y) * dy_laser_inv; + + amrex::Real aabs = 0; + + // interpolate from laser grid to fields grid + for (int iy=0; iy<=interp_order; ++iy) { + for (int ix=0; ix<=interp_order; ++ix) { + auto [shape_x, cell_x] = + compute_single_shape_factor(xmid, ix); + auto [shape_y, cell_y] = + compute_single_shape_factor(ymid, iy); + + if (x_lo <= cell_x && cell_x <= x_hi && y_lo <= cell_y && cell_y <= y_hi) { + aabs += shape_x*shape_y*abssq(laser_arr(cell_x, cell_y, n00j00_r), + laser_arr(cell_x, cell_y, n00j00_i)); + } + } + } + + field_arr(i,j) = aabs; }); } // interpolate aabs to higher MR levels for (int lev=1; lev laser_arr_chi = m_slices.array(mfi, WhichLaserSlice::chi_initial); + + // put chi from the plasma density function on the laser grid as if it were deposited there, + // this works even outside the field grid + // note that the effect of temperature / non-zero u is ignored here + for (auto& plasma : multi_plasma.m_all_plasmas) { + + const PhysConst pc = get_phys_const(); + const amrex::Real c_t = pc.c * Hipace::m_physical_time; + amrex::Real chi_factor = plasma.GetCharge() * plasma.GetCharge() * pc.mu0 / plasma.GetMass(); + if (plasma.m_can_ionize) { + chi_factor *= plasma.m_init_ion_lev * plasma.m_init_ion_lev; + } + + auto density_func = plasma.m_density_func; + + const amrex::Real poff_laser_x = GetPosOffset(0, m_laser_geom_3D, m_laser_geom_3D.Domain()); + const amrex::Real poff_laser_y = GetPosOffset(1, m_laser_geom_3D, m_laser_geom_3D.Domain()); + + const amrex::Real dx_laser = m_laser_geom_3D.CellSize(0); + const amrex::Real dy_laser = m_laser_geom_3D.CellSize(1); + + amrex::ParallelFor(mfi.growntilebox(), + [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept { + const amrex::Real x = i * dx_laser + poff_laser_x; + const amrex::Real y = j * dy_laser + poff_laser_y; + + laser_arr_chi(i, j) += density_func(x, y, c_t) * chi_factor; + }); + } + } +} + +void +MultiLaser::InterpolateChi (const Fields& fields, amrex::Geometry const& geom_field_lev0) +{ + HIPACE_PROFILE("MultiLaser::InterpolateChi()"); + + for ( amrex::MFIter mfi(m_slices, DfltMfi); mfi.isValid(); ++mfi ){ + Array3 laser_arr = m_slices.array(mfi); + Array2 field_arr_chi = + fields.getSlices(0).array(mfi, Comps[WhichSlice::This]["chi"]); + + const amrex::Real poff_laser_x = GetPosOffset(0, m_laser_geom_3D, m_laser_geom_3D.Domain()); + const amrex::Real poff_laser_y = GetPosOffset(1, m_laser_geom_3D, m_laser_geom_3D.Domain()); + const amrex::Real poff_field_x = GetPosOffset(0, geom_field_lev0, geom_field_lev0.Domain()); + const amrex::Real poff_field_y = GetPosOffset(1, geom_field_lev0, geom_field_lev0.Domain()); + + const amrex::Real dx_laser = m_laser_geom_3D.CellSize(0); + const amrex::Real dy_laser = m_laser_geom_3D.CellSize(1); + const amrex::Real dx_laser_inv = m_laser_geom_3D.InvCellSize(0); + const amrex::Real dy_laser_inv = m_laser_geom_3D.InvCellSize(1); + const amrex::Real dx_field = geom_field_lev0.CellSize(0); + const amrex::Real dy_field = geom_field_lev0.CellSize(1); + const amrex::Real dx_field_inv = geom_field_lev0.InvCellSize(0); + const amrex::Real dy_field_inv = geom_field_lev0.InvCellSize(1); + + amrex::Box field_box = fields.getSlices(0)[mfi].box(); + // Even in the valid domain, + // chi near the boundaries is incorrect due to >0 deposition order. + field_box.grow(-2*Fields::m_slices_nguards); + + const amrex::Real pos_x_lo = field_box.smallEnd(0) * dx_field + poff_field_x; + const amrex::Real pos_x_hi = field_box.bigEnd(0) * dx_field + poff_field_x; + const amrex::Real pos_y_lo = field_box.smallEnd(1) * dy_field + poff_field_y; + const amrex::Real pos_y_hi = field_box.bigEnd(1) * dy_field + poff_field_y; + + // the indexes of the laser box where the fields box ends + const int x_lo = amrex::Math::ceil((pos_x_lo - poff_laser_x) * dx_laser_inv); + const int x_hi = amrex::Math::floor((pos_x_hi - poff_laser_x) * dx_laser_inv); + const int y_lo = amrex::Math::ceil((pos_y_lo - poff_laser_y) * dy_laser_inv); + const int y_hi = amrex::Math::floor((pos_y_hi - poff_laser_y) * dy_laser_inv); + + amrex::ParallelFor( + amrex::TypeList>{}, + {m_interp_order}, + mfi.growntilebox(), + [=] AMREX_GPU_DEVICE(int i, int j, int, auto interp_order) noexcept { + const amrex::Real x = i * dx_laser + poff_laser_x; + const amrex::Real y = j * dy_laser + poff_laser_y; + + const amrex::Real xmid = (x - poff_field_x) * dx_field_inv; + const amrex::Real ymid = (y - poff_field_y) * dy_field_inv; + + amrex::Real chi = 0; + + if (x_lo <= i && i <= x_hi && y_lo <= j && j <= y_hi) { + // interpolate chi from fields to laser + for (int iy=0; iy<=interp_order; ++iy) { + for (int ix=0; ix<=interp_order; ++ix) { + auto [shape_x, cell_x] = + compute_single_shape_factor(xmid, ix); + auto [shape_y, cell_y] = + compute_single_shape_factor(ymid, iy); + + chi += shape_x*shape_y*field_arr_chi(cell_x, cell_y); + } + } + } else { + // get initial chi outside the fields box + chi = laser_arr(i, j, WhichLaserSlice::chi_initial); + } + + laser_arr(i, j, WhichLaserSlice::chi) = chi; + }); + } +} + +void +MultiLaser::AdvanceSlice (const int islice, const Fields& fields, amrex::Real dt, int step, + amrex::Geometry const& geom_field_lev0) +{ + + if (!UseLaser(islice)) return; + + InterpolateChi(fields, geom_field_lev0); if (m_solver_type == "multigrid") { - AdvanceSliceMG(fields, dt, step); + AdvanceSliceMG(dt, step); } else if (m_solver_type == "fft") { - AdvanceSliceFFT(fields, dt, step); + AdvanceSliceFFT(dt, step); } else { amrex::Abort("laser.solver_type must be fft or multigrid"); } } void -MultiLaser::AdvanceSliceMG (const Fields& fields, amrex::Real dt, int step) +MultiLaser::AdvanceSliceMG (amrex::Real dt, int step) { HIPACE_PROFILE("MultiLaser::AdvanceSliceMG()"); @@ -552,11 +775,6 @@ MultiLaser::AdvanceSliceMG (const Fields& fields, amrex::Real dt, int step) Array3 rhs_mg_arr = rhs_mg.array(); Array3 acoeff_real_arr = acoeff_real.array(); - constexpr int lev = 0; - const amrex::FArrayBox& isl_fab = fields.getSlices(lev)[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.; @@ -651,7 +869,7 @@ MultiLaser::AdvanceSliceMG (const Fields& fields, amrex::Real dt, int step) const Complex anp1jp1 = arr(i, j, np1jp1_r) + I * arr(i, j, np1jp1_i); const Complex anp1jp2 = arr(i, j, np1jp2_r) + I * arr(i, j, np1jp2_i); acoeff_real_arr(i,j,0) = do_avg_rhs ? - acoeff_real_scalar + isl_arr(i,j,chi) : acoeff_real_scalar; + acoeff_real_scalar + arr(i, j, chi) : acoeff_real_scalar; Complex rhs; if (step == 0) { @@ -665,9 +883,9 @@ MultiLaser::AdvanceSliceMG (const Fields& fields, amrex::Real dt, int step) - lapA + ( -6._rt/(c*dt*dz) + 4._rt*I*djn/(c*dt) + I*4._rt*k0/(c*dt) ) * an00j00; if (do_avg_rhs) { - rhs += isl_arr(i,j,chi) * an00j00; + rhs += arr(i, j, chi) * an00j00; } else { - rhs += isl_arr(i,j,chi) * an00j00 * 2._rt; + rhs += arr(i, j, chi) * an00j00 * 2._rt; } } else { const Complex anm1jp1 = arr(i, j, nm1jp1_r) + I * arr(i, j, nm1jp1_i); @@ -680,9 +898,9 @@ MultiLaser::AdvanceSliceMG (const Fields& fields, amrex::Real dt, int step) - lapA + ( -3._rt/(c*dt*dz) + 2._rt*I*djn/(c*dt) + 2._rt/(c*c*dt*dt) + I*2._rt*k0/(c*dt) ) * anm1j00; if (do_avg_rhs) { - rhs += isl_arr(i,j,chi) * anm1j00; + rhs += arr(i, j, chi) * anm1j00; } else { - rhs += isl_arr(i,j,chi) * an00j00 * 2._rt; + rhs += arr(i, j, chi) * an00j00 * 2._rt; } } rhs_mg_arr(i,j,0) = rhs.real(); @@ -703,7 +921,7 @@ MultiLaser::AdvanceSliceMG (const Fields& fields, amrex::Real dt, int step) } void -MultiLaser::AdvanceSliceFFT (const Fields& fields, const amrex::Real dt, int step) +MultiLaser::AdvanceSliceFFT (const amrex::Real dt, int step) { HIPACE_PROFILE("MultiLaser::AdvanceSliceFFT()"); @@ -737,11 +955,6 @@ MultiLaser::AdvanceSliceFFT (const Fields& fields, const amrex::Real dt, int ste Array3 arr = m_slices.array(mfi); - constexpr int lev = 0; - const amrex::FArrayBox& isl_fab = fields.getSlices(lev)[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); @@ -839,7 +1052,7 @@ MultiLaser::AdvanceSliceFFT (const Fields& fields, const amrex::Real dt, int ste rhs = + 8._rt/(c*dt*dz)*(-anp1jp1+an00jp1)*exp1 + 2._rt/(c*dt*dz)*(+anp1jp2-an00jp2)*exp2 - + 2._rt * isl_arr(i,j,chi) * an00j00 + + 2._rt * arr(i, j, chi) * an00j00 - lapA + ( -6._rt/(c*dt*dz) + 4._rt*I*djn/(c*dt) + I*4._rt*k0/(c*dt) ) * an00j00; } else { @@ -850,7 +1063,7 @@ MultiLaser::AdvanceSliceFFT (const Fields& fields, const amrex::Real dt, int ste + 4._rt/(c*dt*dz)*(-anp1jp1+anm1jp1)*exp1 + 1._rt/(c*dt*dz)*(+anp1jp2-anm1jp2)*exp2 - 4._rt/(c*c*dt*dt)*an00j00 - + 2._rt * isl_arr(i,j,chi) * an00j00 + + 2._rt * arr(i, j, chi) * an00j00 - lapA + ( -3._rt/(c*dt*dz) + 2._rt*I*djn/(c*dt) + 2._rt/(c*c*dt*dt) + I*2._rt*k0/(c*dt) ) * anm1j00; } @@ -902,10 +1115,8 @@ MultiLaser::AdvanceSliceFFT (const Fields& fields, const amrex::Real dt, int ste } void -MultiLaser::InitLaserSlice (const amrex::Geometry& geom, const int islice, const int comp) +MultiLaser::InitLaserSlice (const int islice, const int comp) { - if (!m_use_laser) return; - HIPACE_PROFILE("MultiLaser::InitLaserSlice()"); using namespace amrex::literals; @@ -916,9 +1127,10 @@ MultiLaser::InitLaserSlice (const amrex::Geometry& geom, const int islice, const 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]}; + const amrex::Real poff_x = GetPosOffset(0, m_laser_geom_3D, m_laser_geom_3D.Domain()); + const amrex::Real poff_y = GetPosOffset(1, m_laser_geom_3D, m_laser_geom_3D.Domain()); + const amrex::Real poff_z = GetPosOffset(2, m_laser_geom_3D, m_laser_geom_3D.Domain()); + const amrex::GpuArray dx_arr = m_laser_geom_3D.CellSizeArray(); #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) @@ -943,9 +1155,9 @@ MultiLaser::InitLaserSlice (const amrex::Geometry& geom, const int islice, const bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { - const amrex::Real z = plo[2] + (islice+0.5_rt)*dx_arr[2] - z0; - const amrex::Real x = (i+0.5_rt)*dx_arr[0]+plo[0]-x0; - const amrex::Real y = (j+0.5_rt)*dx_arr[1]+plo[1]-y0; + const amrex::Real x = i * dx_arr[0] + poff_x - x0; + const amrex::Real y = j * dx_arr[1] + poff_y - y0; + const amrex::Real z = islice * dx_arr[2] + poff_z - z0; // Coordinate rotation in yz plane for a laser propagating at an angle. const amrex::Real yp=std::cos(propagation_angle_yz+PFT_yz)*y-std::sin(propagation_angle_yz+PFT_yz)*z; const amrex::Real zp=std::sin(propagation_angle_yz+PFT_yz)*y+std::cos(propagation_angle_yz+PFT_yz)*z; @@ -973,9 +1185,10 @@ MultiLaser::InitLaserSlice (const amrex::Geometry& geom, const int islice, const } void -MultiLaser::InSituComputeDiags (int step, amrex::Real time, int islice, const amrex::Geometry& geom3D, +MultiLaser::InSituComputeDiags (int step, amrex::Real time, int islice, int max_step, amrex::Real max_time) { + if (!UseLaser(islice)) return; if (!utils::doDiagnostics(m_insitu_period, step, max_step, time, max_time)) return; HIPACE_PROFILE("MultiLaser::InSituComputeDiags()"); @@ -985,17 +1198,18 @@ MultiLaser::InSituComputeDiags (int step, amrex::Real time, int islice, const am AMREX_ALWAYS_ASSERT(m_insitu_rdata.size()>0 && m_insitu_sum_rdata.size()>0 && m_insitu_cdata.size()>0); - const int nslices = geom3D.Domain().length(2); - const amrex::Real poff_x = GetPosOffset(0, geom3D, geom3D.Domain()); - const amrex::Real poff_y = GetPosOffset(1, geom3D, geom3D.Domain()); - const amrex::Real dx = geom3D.CellSize(0); - const amrex::Real dy = geom3D.CellSize(1); - const amrex::Real dxdydz = dx * dy * geom3D.CellSize(2); - - const int xmid_lo = geom3D.Domain().smallEnd(0) + (geom3D.Domain().length(0) - 1) / 2; - const int xmid_hi = geom3D.Domain().smallEnd(0) + (geom3D.Domain().length(0)) / 2; - const int ymid_lo = geom3D.Domain().smallEnd(1) + (geom3D.Domain().length(1) - 1) / 2; - const int ymid_hi = geom3D.Domain().smallEnd(1) + (geom3D.Domain().length(1)) / 2; + const int nslices = m_laser_geom_3D.Domain().length(2); + const int laser_slice = islice - m_laser_geom_3D.Domain().smallEnd(2); + const amrex::Real poff_x = GetPosOffset(0, m_laser_geom_3D, m_laser_geom_3D.Domain()); + const amrex::Real poff_y = GetPosOffset(1, m_laser_geom_3D, m_laser_geom_3D.Domain()); + const amrex::Real dx = m_laser_geom_3D.CellSize(0); + const amrex::Real dy = m_laser_geom_3D.CellSize(1); + const amrex::Real dxdydz = dx * dy * m_laser_geom_3D.CellSize(2); + + const int xmid_lo = m_laser_geom_3D.Domain().smallEnd(0) + (m_laser_geom_3D.Domain().length(0) - 1) / 2; + const int xmid_hi = m_laser_geom_3D.Domain().smallEnd(0) + (m_laser_geom_3D.Domain().length(0)) / 2; + const int ymid_lo = m_laser_geom_3D.Domain().smallEnd(1) + (m_laser_geom_3D.Domain().length(1) - 1) / 2; + const int ymid_hi = m_laser_geom_3D.Domain().smallEnd(1) + (m_laser_geom_3D.Domain().length(1)) / 2; const amrex::Real mid_factor = (xmid_lo == xmid_hi ? 1._rt : 0.5_rt) * (ymid_lo == ymid_hi ? 1._rt : 0.5_rt); @@ -1037,10 +1251,10 @@ MultiLaser::InSituComputeDiags (int step, amrex::Real time, int islice, const am amrex::constexpr_for<0, m_insitu_nrp>( [&] (auto idx) { if (idx == 0) { - m_insitu_rdata[islice + idx * nslices] = amrex::get(a); + m_insitu_rdata[laser_slice + idx * nslices] = amrex::get(a); m_insitu_sum_rdata[idx] = std::max(m_insitu_sum_rdata[idx], amrex::get(a)); } else { - m_insitu_rdata[islice + idx * nslices] = amrex::get(a)*dxdydz; + m_insitu_rdata[laser_slice + idx * nslices] = amrex::get(a)*dxdydz; m_insitu_sum_rdata[idx] += amrex::get(a)*dxdydz; } } @@ -1048,15 +1262,15 @@ MultiLaser::InSituComputeDiags (int step, amrex::Real time, int islice, const am amrex::constexpr_for<0, m_insitu_ncp>( [&] (auto idx) { - m_insitu_cdata[islice + idx * nslices] = amrex::get(a) * mid_factor; + m_insitu_cdata[laser_slice + idx * nslices] = amrex::get(a) * mid_factor; } ); } void -MultiLaser::InSituWriteToFile (int step, amrex::Real time, const amrex::Geometry& geom3D, - int max_step, amrex::Real max_time) +MultiLaser::InSituWriteToFile (int step, amrex::Real time, int max_step, amrex::Real max_time) { + if (!m_use_laser) return; if (!utils::doDiagnostics(m_insitu_period, step, max_step, time, max_time)) return; HIPACE_PROFILE("MultiLaser::InSituWriteToFile()"); @@ -1074,7 +1288,7 @@ MultiLaser::InSituWriteToFile (int step, amrex::Real time, const amrex::Geometry std::ofstream ofs{m_insitu_file_prefix + "/reduced_laser." + pad_rank_num + ".txt", std::ofstream::out | std::ofstream::app | std::ofstream::binary}; - const int nslices_int = geom3D.Domain().length(2); + const int nslices_int = m_laser_geom_3D.Domain().length(2); const std::size_t nslices = static_cast(nslices_int); const int is_normalized_units = Hipace::m_normalized_units; @@ -1084,8 +1298,8 @@ MultiLaser::InSituWriteToFile (int step, amrex::Real time, const amrex::Geometry {"time" , &time}, {"step" , &step}, {"n_slices" , &nslices_int}, - {"z_lo" , &geom3D.ProbLo()[2]}, - {"z_hi" , &geom3D.ProbHi()[2]}, + {"z_lo" , &m_laser_geom_3D.ProbLo()[2]}, + {"z_hi" , &m_laser_geom_3D.ProbHi()[2]}, {"is_normalized_units", &is_normalized_units}, {"max(|a|^2)" , &m_insitu_rdata[0], nslices}, {"[|a|^2]" , &m_insitu_rdata[1*nslices], nslices}, diff --git a/src/utils/MultiBuffer.H b/src/utils/MultiBuffer.H index 97b19be127..dfb872fdd5 100644 --- a/src/utils/MultiBuffer.H +++ b/src/utils/MultiBuffer.H @@ -18,7 +18,7 @@ class MultiBuffer public: // initialize MultiBuffer and open initial receive requests - void initialize (int nslices, MultiBeam& beams, bool use_laser, amrex::Box laser_box); + void initialize (int nslices, MultiBeam& beams, MultiLaser& laser); // receive data from previous rank and unpack it into MultiBeam and MultiLaser void get_data (int slice, MultiBeam& beams, MultiLaser& laser, int beam_slice); @@ -113,9 +113,7 @@ private: bool m_async_memcpy = true; int m_nslices = 0; int m_nbeams = 0; - bool m_use_laser = false; int m_laser_ncomp = 4; - amrex::Box m_laser_slice_box {}; /** How many slices of beam particles can be received in advance */ int m_max_leading_slices = std::numeric_limits::max(); /** How many slices of beam particles can be stored before being sent */ @@ -144,10 +142,10 @@ private: void make_progress (int slice, bool is_blocking, int current_slice); // write MultiBeam sizes into the metadata array - void write_metadata (int slice, MultiBeam& beams, int beam_slice); + void write_metadata (int slice, MultiBeam& beams, MultiLaser& laser, int beam_slice); // helper function to get location of individual arrays inside a buffer - std::size_t get_buffer_offset (int slice, offset_type type, MultiBeam& beams, + std::size_t get_buffer_offset (int slice, offset_type type, MultiBeam& beams, MultiLaser& laser, int ibeam, int comp); // copy gpu array into buffer at buffer_offset, either dtoh or dtod diff --git a/src/utils/MultiBuffer.cpp b/src/utils/MultiBuffer.cpp index 6be36b9443..0da539264c 100644 --- a/src/utils/MultiBuffer.cpp +++ b/src/utils/MultiBuffer.cpp @@ -52,7 +52,7 @@ void MultiBuffer::free_buffer (int slice) { m_datanodes[slice].m_buffer_size = 0; } -void MultiBuffer::initialize (int nslices, MultiBeam& beams, bool use_laser, amrex::Box laser_box) { +void MultiBuffer::initialize (int nslices, MultiBeam& beams, MultiLaser& laser) { amrex::ParmParse pp("comms_buffer"); @@ -62,8 +62,6 @@ void MultiBuffer::initialize (int nslices, MultiBeam& beams, bool use_laser, amr m_nslices = nslices; m_nbeams = beams.get_nbeams(); - m_use_laser = use_laser; - m_laser_slice_box = laser_box; m_rank_send_to = (rank_id + 1) % n_ranks; m_rank_receive_from = (rank_id - 1 + n_ranks) % n_ranks; @@ -119,8 +117,8 @@ void MultiBuffer::initialize (int nslices, MultiBeam& beams, bool use_laser, amr } } - if (m_use_laser) { - size_estimate += m_laser_slice_box.d_numPts() * nslices + if (laser.UseLaser()) { + size_estimate += laser.GetLaserGeom().Domain().numPts() * m_laser_ncomp * sizeof(amrex::Real); } @@ -395,7 +393,7 @@ void MultiBuffer::get_data (int slice, MultiBeam& beams, MultiLaser& laser, int for (int b = 0; b < m_nbeams; ++b) { beams.getBeam(b).intializeSlice(slice, beam_slice); } - if (m_use_laser) { + if (laser.UseLaser(slice)) { using namespace WhichLaserSlice; const int laser_comp = (beam_slice == WhichBeamSlice::Next) ? n00jp2_r : n00j00_r; laser.InitSliceEnvelope(slice, laser_comp); @@ -445,7 +443,7 @@ void MultiBuffer::put_data (int slice, MultiBeam& beams, MultiLaser& laser, int m_datanodes[slice].m_metadata_progress = comm_progress::sim_completed; } else { // pack and asynchronously send buffer - write_metadata(slice, beams, beam_slice); + write_metadata(slice, beams, laser, beam_slice); m_datanodes[slice].m_metadata_progress = comm_progress::ready_to_send; if (m_async_memcpy) { if (slice < m_nslices - 1) { @@ -600,12 +598,12 @@ void MultiBuffer::put_time (amrex::Real time) { #endif } -void MultiBuffer::write_metadata (int slice, MultiBeam& beams, int beam_slice) { +void MultiBuffer::write_metadata (int slice, MultiBeam& beams, MultiLaser& laser, int beam_slice) { for (int b = 0; b < m_nbeams; ++b) { // write number of beam particles (per beam) get_metadata_location(slice)[b + 1] = beams.getBeam(b).getNumParticles(beam_slice); } - std::size_t offset = get_buffer_offset(slice, offset_type::total, beams, 0, 0); + std::size_t offset = get_buffer_offset(slice, offset_type::total, beams, laser, 0, 0); // write total buffer size get_metadata_location(slice)[0] = (offset+sizeof(storage_type)-1) / sizeof(storage_type); m_datanodes[slice].m_buffer_size = get_metadata_location(slice)[0]; @@ -615,7 +613,7 @@ void MultiBuffer::write_metadata (int slice, MultiBeam& beams, int beam_slice) { } std::size_t MultiBuffer::get_buffer_offset (int slice, offset_type type, MultiBeam& beams, - int ibeam, int comp) { + MultiLaser& laser, int ibeam, int comp) { // calculate offset for each chunk of data in one place // to ensure consistency between packing and unpacking std::size_t offset = 0; @@ -656,12 +654,12 @@ std::size_t MultiBuffer::get_buffer_offset (int slice, offset_type type, MultiBe } // add offset for laser, if used - if (m_use_laser) { + if (laser.UseLaser(slice)) { for (int lcomp = 0; lcomp < m_laser_ncomp; ++lcomp) { if (type == offset_type::laser && lcomp == comp) { return offset; } - offset += m_laser_slice_box.numPts() * sizeof(amrex::Real); + offset += laser.getSlices()[0].box().numPts() * sizeof(amrex::Real); } } @@ -758,7 +756,8 @@ void MultiBuffer::pack_data (int slice, MultiBeam& beams, MultiLaser& laser, int if (beam.communicateIdCpuComponent()) { // only pack idcpu component if it should be communicated - memcpy_to_buffer(slice, get_buffer_offset(slice, offset_type::beam_idcpu, beams, b, 0), + memcpy_to_buffer(slice, get_buffer_offset(slice, offset_type::beam_idcpu, + beams, laser, b, 0), soa.GetIdCPUData().dataPtr(), num_particles * sizeof(std::uint64_t)); } @@ -766,7 +765,8 @@ void MultiBuffer::pack_data (int slice, MultiBeam& beams, MultiLaser& laser, int for (int rcomp = 0; rcomp < beam.numRealComponents(); ++rcomp) { // only pack real component if it should be communicated if (beam.communicateRealComponent(rcomp)) { - memcpy_to_buffer(slice, get_buffer_offset(slice, offset_type::beam_real, beams, b, rcomp), + memcpy_to_buffer(slice, get_buffer_offset(slice, offset_type::beam_real, + beams, laser, b, rcomp), soa.GetRealData(rcomp).dataPtr(), num_particles * sizeof(amrex::Real)); } @@ -775,23 +775,24 @@ void MultiBuffer::pack_data (int slice, MultiBeam& beams, MultiLaser& laser, int for (int icomp = 0; icomp < beam.numIntComponents(); ++icomp) { // only pack int component if it should be communicated if (beam.communicateIntComponent(icomp)) { - memcpy_to_buffer(slice, get_buffer_offset(slice, offset_type::beam_int, beams, b, icomp), + memcpy_to_buffer(slice, get_buffer_offset(slice, offset_type::beam_int, + beams, laser, b, icomp), soa.GetIntData(icomp).dataPtr(), num_particles * sizeof(int)); } } } - if (m_use_laser) { + if (laser.UseLaser(slice)) { using namespace WhichLaserSlice; const int laser_comp_0_1 = (beam_slice == WhichBeamSlice::Next) ? np1jp2_r : np1j00_r; const int laser_comp_2_3 = (beam_slice == WhichBeamSlice::Next) ? n00jp2_r : n00j00_r; // copy real and imag components in one operation - memcpy_to_buffer(slice, get_buffer_offset(slice, offset_type::laser, beams, 0, 0), + memcpy_to_buffer(slice, get_buffer_offset(slice, offset_type::laser, beams, laser, 0, 0), laser.getSlices()[0].dataPtr(laser_comp_0_1), - 2 * m_laser_slice_box.numPts() * sizeof(amrex::Real)); - memcpy_to_buffer(slice, get_buffer_offset(slice, offset_type::laser, beams, 0, 2), + 2 * laser.getSlices()[0].box().numPts() * sizeof(amrex::Real)); + memcpy_to_buffer(slice, get_buffer_offset(slice, offset_type::laser, beams, laser, 0, 2), laser.getSlices()[0].dataPtr(laser_comp_2_3), - 2 * m_laser_slice_box.numPts() * sizeof(amrex::Real)); + 2 * laser.getSlices()[0].box().numPts() * sizeof(amrex::Real)); } amrex::Gpu::streamSynchronize(); for (int b = 0; b < m_nbeams; ++b) { @@ -809,7 +810,8 @@ void MultiBuffer::unpack_data (int slice, MultiBeam& beams, MultiLaser& laser, i if (beam.communicateIdCpuComponent()) { // only undpack idcpu component if it should be communicated - memcpy_from_buffer(slice, get_buffer_offset(slice, offset_type::beam_idcpu, beams, b, 0), + memcpy_from_buffer(slice, get_buffer_offset(slice, offset_type::beam_idcpu, + beams, laser, b, 0), soa.GetIdCPUData().dataPtr(), num_particles * sizeof(std::uint64_t)); } else { @@ -824,7 +826,8 @@ void MultiBuffer::unpack_data (int slice, MultiBeam& beams, MultiLaser& laser, i for (int rcomp = 0; rcomp < beam.numRealComponents(); ++rcomp) { if (beam.communicateRealComponent(rcomp)) { // only unpack real component if it should be communicated - memcpy_from_buffer(slice, get_buffer_offset(slice, offset_type::beam_real, beams, b, rcomp), + memcpy_from_buffer(slice, get_buffer_offset(slice, offset_type::beam_real, + beams, laser, b, rcomp), soa.GetRealData(rcomp).dataPtr(), num_particles * sizeof(amrex::Real)); } else { @@ -839,7 +842,8 @@ void MultiBuffer::unpack_data (int slice, MultiBeam& beams, MultiLaser& laser, i for (int icomp = 0; icomp < beam.numIntComponents(); ++icomp) { if (beam.communicateIntComponent(icomp)) { // only unpack int component if it should be communicated - memcpy_from_buffer(slice, get_buffer_offset(slice, offset_type::beam_int, beams, b, icomp), + memcpy_from_buffer(slice, get_buffer_offset(slice, offset_type::beam_int, + beams, laser, b, icomp), soa.GetIntData(icomp).dataPtr(), num_particles * sizeof(int)); } else { @@ -851,17 +855,17 @@ void MultiBuffer::unpack_data (int slice, MultiBeam& beams, MultiLaser& laser, i } } } - if (m_use_laser) { + if (laser.UseLaser(slice)) { using namespace WhichLaserSlice; const int laser_comp_0_1 = (beam_slice == WhichBeamSlice::Next) ? n00jp2_r : n00j00_r; const int laser_comp_2_3 = (beam_slice == WhichBeamSlice::Next) ? nm1jp2_r : nm1j00_r; // copy real and imag components in one operation - memcpy_from_buffer(slice, get_buffer_offset(slice, offset_type::laser, beams, 0, 0), + memcpy_from_buffer(slice, get_buffer_offset(slice, offset_type::laser, beams, laser, 0, 0), laser.getSlices()[0].dataPtr(laser_comp_0_1), - 2 * m_laser_slice_box.numPts() * sizeof(amrex::Real)); - memcpy_from_buffer(slice, get_buffer_offset(slice, offset_type::laser, beams, 0, 2), + 2 * laser.getSlices()[0].box().numPts() * sizeof(amrex::Real)); + memcpy_from_buffer(slice, get_buffer_offset(slice, offset_type::laser, beams, laser, 0, 2), laser.getSlices()[0].dataPtr(laser_comp_2_3), - 2 * m_laser_slice_box.numPts() * sizeof(amrex::Real)); + 2 * laser.getSlices()[0].box().numPts() * sizeof(amrex::Real)); } amrex::Gpu::streamSynchronize(); }