From fb7583127b6eb0a384634b34adf07712ca2c4220 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Mon, 22 Apr 2024 22:13:51 +0200 Subject: [PATCH 01/24] Use correct geometry for Laser diagnostics --- src/fields/Fields.cpp | 3 ++- src/laser/MultiLaser.H | 3 +++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/fields/Fields.cpp b/src/fields/Fields.cpp index 3d1ea83b18..67495c6df1 100644 --- a/src/fields/Fields.cpp +++ b/src/fields/Fields.cpp @@ -488,7 +488,8 @@ Fields::Copy (const int lev, const int i_slice, FieldDiagnosticData& fd, auto& slice_mf = m_slices[lev]; auto slice_func = interpolated_field_xy{{slice_mf}, calc_geom}; auto& laser_mf = multi_laser.getSlices(); - auto laser_func = interpolated_field_xy{{laser_mf}, calc_geom}; + auto laser_func = interpolated_field_xy{{laser_mf}, multi_laser.GetLaserGeom()}; #ifdef AMREX_USE_GPU // This async copy happens on the same stream as the ParallelFor below, which uses the copied array. diff --git a/src/laser/MultiLaser.H b/src/laser/MultiLaser.H index 0495fa68fb..e0619d331a 100644 --- a/src/laser/MultiLaser.H +++ b/src/laser/MultiLaser.H @@ -224,6 +224,9 @@ public: /** Get the central wavelength */ amrex::Real GetLambda0 () const { return m_lambda0; } + /** 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 */ private: From 60df9d0c45739dd679b5c129a9142798775b55ea Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Wed, 24 Apr 2024 22:17:45 +0200 Subject: [PATCH 02/24] separate grid for laser --- src/Hipace.cpp | 13 +-- src/laser/MultiLaser.H | 54 ++++++------ src/laser/MultiLaser.cpp | 174 +++++++++++++++++++++++++-------------- 3 files changed, 143 insertions(+), 98 deletions(-) diff --git a/src/Hipace.cpp b/src/Hipace.cpp index faca043427..30a11834cf 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -156,6 +156,9 @@ Hipace::Hipace () : MakeGeometry(); + // use level 0 as default for laser geometry + m_multi_laser.MakeLaserGeometry(m_3D_geom[0]); + m_use_laser = m_multi_laser.m_use_laser; queryWithParser(pph, "collisions", m_collision_names); @@ -190,12 +193,10 @@ Hipace::InitData () amrex::Print() << "using the leapfrog plasma particle pusher\n"; #endif + m_multi_laser.InitData(); + for (int lev=0; lev - void GetEnvelopeFromFile (const amrex::Geometry& gm); + void GetEnvelopeFromFile (); /** \brief Shift 2D slices in zeta */ @@ -166,6 +160,8 @@ public: void UpdateLaserAabs (const int current_N_level, Fields& fields, amrex::Vector const& geom); + void InterpolateChi (const Fields& fields, amrex::Geometry const& geom_field_lev0); + /** Wrapper function to advance a laser slice by 1 time step. * \param[in] fields Field object * \param[in] dt time step of the simulation @@ -176,50 +172,44 @@ public: /** 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; } @@ -239,8 +229,14 @@ 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; /** 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 13be96812c..f70fbde7ff 100644 --- a/src/laser/MultiLaser.cpp +++ b/src/laser/MultiLaser.cpp @@ -83,9 +83,53 @@ 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"); + + 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)}; + + queryWithParser(pp, "n_cell", n_cells_laser); + queryWithParser(pp, "patch_lo", patch_lo_laser); + queryWithParser(pp, "patch_hi", patch_hi_laser); + + const amrex::Real pos_offset_z = GetPosOffset(2, field_geom_3D, field_geom_3D.Domain()); + + const int zeta_lo = + int(amrex::Math::round((patch_lo_laser[2] - pos_offset_z) * field_geom_3D.InvCellSize(2))); + + const int zeta_hi = + 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; + + 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)}; + + m_laser_geom_3D.define(domain_3D_laser, amrex::RealBox(patch_lo_laser, patch_hi_laser), + 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; @@ -96,13 +140,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()); @@ -172,7 +213,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 @@ -211,13 +252,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 @@ -259,14 +300,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 @@ -274,7 +314,7 @@ MultiLaser::GetEnvelopeFromFileHelper (const amrex::Geometry& gm) { template void -MultiLaser::GetEnvelopeFromFile (const amrex::Geometry& gm) { +MultiLaser::GetEnvelopeFromFile () { using namespace amrex::literals; @@ -283,7 +323,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]; @@ -321,13 +361,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); @@ -561,23 +601,42 @@ MultiLaser::UpdateLaserAabs (const int current_N_level, Fields& fields, } } +void +MultiLaser::InterpolateChi (const Fields& fields, amrex::Geometry const&) +{ + HIPACE_PROFILE("MultiLaser::InterpolateChi()"); + + for ( amrex::MFIter mfi(m_slices, DfltMfi); mfi.isValid(); ++mfi ){ + + Array2 laser_arr_chi = m_slices.array(mfi, WhichLaserSlice::chi); + Array2 field_arr_chi = fields.getSlices(0).array(mfi, Comps[WhichSlice::This]["chi"]); + + amrex::ParallelFor(mfi.growntilebox(), + [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept { + laser_arr_chi(i, j) = field_arr_chi(i, j); + }); + } +} + void MultiLaser::AdvanceSlice (const Fields& fields, amrex::Real dt, int step) { if (!m_use_laser) return; + InterpolateChi(fields, m_laser_geom_3D); + 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()"); @@ -616,11 +675,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.; @@ -715,7 +769,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) { @@ -729,9 +783,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); @@ -744,9 +798,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(); @@ -767,7 +821,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()"); @@ -801,11 +855,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); @@ -903,7 +952,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 { @@ -914,7 +963,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; } @@ -1023,7 +1072,7 @@ 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; @@ -1037,8 +1086,8 @@ 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 auto plo = m_laser_geom_3D.ProbLoArray(); + amrex::Real const * const dx = m_laser_geom_3D.CellSize(); const amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; #ifdef AMREX_USE_OMP @@ -1093,7 +1142,7 @@ 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 (!utils::doDiagnostics(m_insitu_period, step, max_step, time, max_time)) return; @@ -1105,17 +1154,17 @@ 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 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); @@ -1174,8 +1223,7 @@ MultiLaser::InSituComputeDiags (int step, amrex::Real time, int islice, const am } 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 (!utils::doDiagnostics(m_insitu_period, step, max_step, time, max_time)) return; HIPACE_PROFILE("MultiLaser::InSituWriteToFile()"); @@ -1194,7 +1242,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; @@ -1204,8 +1252,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}, From 8bb09753b8a16651bfa1ccb5cca3e5273c237020 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Sun, 28 Apr 2024 01:26:24 +0200 Subject: [PATCH 03/24] add min and max for zeta --- src/laser/MultiLaser.cpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/laser/MultiLaser.cpp b/src/laser/MultiLaser.cpp index f70fbde7ff..b280d5300f 100644 --- a/src/laser/MultiLaser.cpp +++ b/src/laser/MultiLaser.cpp @@ -105,11 +105,13 @@ MultiLaser::MakeLaserGeometry (const amrex::Geometry& field_geom_3D) const amrex::Real pos_offset_z = GetPosOffset(2, field_geom_3D, field_geom_3D.Domain()); - const int zeta_lo = - int(amrex::Math::round((patch_lo_laser[2] - pos_offset_z) * field_geom_3D.InvCellSize(2))); + 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 = - int(amrex::Math::round((patch_hi_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; From 5088f4a411256ca7c88a7a9e4c9ae026321786a9 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Sun, 28 Apr 2024 03:22:48 +0200 Subject: [PATCH 04/24] add interpolation --- src/Hipace.cpp | 11 +++--- src/fields/Fields.cpp | 8 +++-- src/laser/MultiLaser.H | 16 ++++++--- src/laser/MultiLaser.cpp | 71 +++++++++++++++++++++++++++++++++------ src/utils/MultiBuffer.H | 2 +- src/utils/MultiBuffer.cpp | 16 ++++----- 6 files changed, 92 insertions(+), 32 deletions(-) diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 30a11834cf..051d35c474 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -212,10 +212,7 @@ Hipace::InitData () m_fields.checkInit(); - m_multi_buffer.initialize(m_3D_geom[0].Domain().length(2), - m_multi_beam, - m_use_laser, - m_use_laser ? m_multi_laser.getSlices()[0].box() : amrex::Box{}); + m_multi_buffer.initialize(m_3D_geom[0].Domain().length(2), m_multi_beam, m_multi_laser); amrex::ParmParse pph("hipace"); bool do_output_input = false; @@ -483,7 +480,7 @@ Hipace::SolveOneSlice (int islice, int step) } // write laser aabs into fields MultiFab - m_multi_laser.UpdateLaserAabs(current_N_level, m_fields, m_3D_geom); + m_multi_laser.UpdateLaserAabs(islice, current_N_level, m_fields, m_3D_geom); // deposit current 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); @@ -612,7 +609,7 @@ Hipace::SolveOneSlice (int islice, int step) m_multi_beam.shiftBeamSlices(); - m_multi_laser.ShiftLaserSlices(); + m_multi_laser.ShiftLaserSlices(islice); } void diff --git a/src/fields/Fields.cpp b/src/fields/Fields.cpp index 3d1ea83b18..fcbf1cfca9 100644 --- a/src/fields/Fields.cpp +++ b/src/fields/Fields.cpp @@ -426,6 +426,8 @@ Fields::Copy (const int lev, const int i_slice, FieldDiagnosticData& fd, constexpr int depos_order_z = 1; constexpr int depos_order_offset = depos_order_z / 2 + 1; + const amrex::Geometry& laser_geom = multi_laser.GetLaserGeom(); + const amrex::Real poff_calc_z = GetPosOffset(2, calc_geom, calc_geom.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()); @@ -507,7 +509,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_nfields > 0 && + calc_geom.smallEnd(2) <= i_slice && i_slice <= calc_geom.bigEnd(2)) { auto slice_array = slice_func.array(mfi); amrex::Array4 diag_array = fd.m_F.array(); amrex::ParallelFor(diag_box, fd.m_nfields, @@ -520,7 +523,8 @@ Fields::Copy (const int lev, const int i_slice, FieldDiagnosticData& fd, }); } - if (fd.m_do_laser) { + if (fd.m_do_laser && + laser_geom.smallEnd(2) <= i_slice && i_slice <= laser_geom.bigEnd(2)) { auto laser_array = laser_func.array(mfi); amrex::Array4 diag_array_laser = fd.m_F_laser.array(); amrex::ParallelFor(diag_box, diff --git a/src/laser/MultiLaser.H b/src/laser/MultiLaser.H index 7289edfde6..2f54344c3d 100644 --- a/src/laser/MultiLaser.H +++ b/src/laser/MultiLaser.H @@ -149,25 +149,28 @@ public: 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 */ - void UpdateLaserAabs (const int current_N_level, Fields& fields, - amrex::Vector const& geom); + void UpdateLaserAabs (const int islice, const int current_N_level, Fields& fields, + amrex::Vector const& field_geom); void InterpolateChi (const Fields& fields, amrex::Geometry const& geom_field_lev0); /** 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. */ - void AdvanceSlice (const Fields& fields, amrex::Real dt, int step); + void AdvanceSlice (const int islice, const Fields& fields, 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. @@ -214,6 +217,11 @@ public: /** Get the central wavelength */ amrex::Real GetLambda0 () const { return m_lambda0; } + bool HasSlice (const int islice) const { + return m_laser_geom_3D.Domain.smallEnd(2) <= islice && + islice <= m_laser_geom_3D.Domain.bigEnd(2); + } + bool m_use_laser {false}; /**< whether a laser is used or not */ private: diff --git a/src/laser/MultiLaser.cpp b/src/laser/MultiLaser.cpp index b280d5300f..51159283cc 100644 --- a/src/laser/MultiLaser.cpp +++ b/src/laser/MultiLaser.cpp @@ -10,6 +10,7 @@ #include "MultiLaser.H" #include "utils/Constants.H" #include "Hipace.H" +#include "particles/particles_utils/ShapeFactors.H" #include "utils/HipaceProfilerWrapper.H" #include "utils/DeprecatedInput.H" #include "utils/InsituUtil.H" @@ -244,6 +245,7 @@ void MultiLaser::InitSliceEnvelope (const int islice, const int comp) { if (!m_use_laser) return; + if (!HasSlice(islice)) return; HIPACE_PROFILE("MultiLaser::InitSliceEnvelope()"); @@ -542,9 +544,10 @@ MultiLaser::GetEnvelopeFromFile () { } void -MultiLaser::ShiftLaserSlices () +MultiLaser::ShiftLaserSlices (const int islice) { if (!m_use_laser) return; + if (!HasSlice(islice)) return; HIPACE_PROFILE("MultiLaser::ShiftLaserSlices()"); @@ -576,30 +579,73 @@ 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; HIPACE_PROFILE("MultiLaser::UpdateLaserAabs()"); + constexpr interp_order = 1; + + if (!HasSlice(islice)) { + 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"]); + 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(mfi.growntilebox(), [=] AMREX_GPU_DEVICE(int i, int j, int) 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; + + 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; lev0); 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); @@ -1208,10 +1258,10 @@ MultiLaser::InSituComputeDiags (int step, amrex::Real time, int islice, 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; } } @@ -1219,7 +1269,7 @@ MultiLaser::InSituComputeDiags (int step, amrex::Real time, int islice, 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; } ); } @@ -1227,6 +1277,7 @@ MultiLaser::InSituComputeDiags (int step, amrex::Real time, int islice, void 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()"); diff --git a/src/utils/MultiBuffer.H b/src/utils/MultiBuffer.H index 97b19be127..ccb8cb5092 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); diff --git a/src/utils/MultiBuffer.cpp b/src/utils/MultiBuffer.cpp index 6be36b9443..670491beed 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,7 +62,7 @@ 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_use_laser = laser.m_use_laser; m_laser_slice_box = laser_box; m_rank_send_to = (rank_id + 1) % n_ranks; @@ -781,17 +781,17 @@ void MultiBuffer::pack_data (int slice, MultiBeam& beams, MultiLaser& laser, int } } } - if (m_use_laser) { + if (m_use_laser && laser.HasSlice(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), laser.getSlices()[0].dataPtr(laser_comp_0_1), - 2 * m_laser_slice_box.numPts() * sizeof(amrex::Real)); + 2 * laser.getSlices()[0].box().numPts() * sizeof(amrex::Real)); memcpy_to_buffer(slice, get_buffer_offset(slice, offset_type::laser, beams, 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) { @@ -851,17 +851,17 @@ void MultiBuffer::unpack_data (int slice, MultiBeam& beams, MultiLaser& laser, i } } } - if (m_use_laser) { + if (m_use_laser && laser.HasSlice(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), laser.getSlices()[0].dataPtr(laser_comp_0_1), - 2 * m_laser_slice_box.numPts() * sizeof(amrex::Real)); + 2 * laser.getSlices()[0].box().numPts() * sizeof(amrex::Real)); memcpy_from_buffer(slice, get_buffer_offset(slice, offset_type::laser, beams, 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(); } From 33bc851b7b2337850f7b1ba5dec8ebe2f0afe512 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Sun, 28 Apr 2024 04:39:37 +0200 Subject: [PATCH 05/24] add more interpolation --- src/Hipace.cpp | 6 ++--- src/laser/MultiLaser.H | 16 ++++++++---- src/laser/MultiLaser.cpp | 53 ++++++++++++++++++++++++++++++--------- src/utils/MultiBuffer.H | 6 ++--- src/utils/MultiBuffer.cpp | 50 +++++++++++++++++++----------------- 5 files changed, 84 insertions(+), 47 deletions(-) diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 051d35c474..c6c2ea8506 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -159,7 +159,7 @@ Hipace::Hipace () : // use level 0 as default for laser geometry m_multi_laser.MakeLaserGeometry(m_3D_geom[0]); - m_use_laser = m_multi_laser.m_use_laser; + m_use_laser = m_multi_laser.UseLaser(); queryWithParser(pph, "collisions", m_collision_names); /** Initialize the collision objects */ @@ -197,7 +197,7 @@ Hipace::InitData () for (int lev=0; lev laser_arr_chi = m_slices.array(mfi, WhichLaserSlice::chi); 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_field_inv = geom_field_lev0.InvCellSize(0); + const amrex::Real dy_field_inv = geom_field_lev0.InvCellSize(1); + + const int x_lo = fields.getSlices(0)[mfi].box().smallEnd(0); + const int x_hi = fields.getSlices(0)[mfi].box().bigEnd(0); + const int y_lo = fields.getSlices(0)[mfi].box().smallEnd(1); + const int y_hi = fields.getSlices(0)[mfi].box().bigEnd(1); + amrex::ParallelFor(mfi.growntilebox(), [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept { - laser_arr_chi(i, j) = field_arr_chi(i, j); + const amrex::Real x = i * dx_laser + m_laser_geom_3D; + const amrex::Real y = j * dy_laser + m_laser_geom_3D; + + 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; + + 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) { + chi += shape_x*shape_y*field_arr_chi(cell_x, cell_y); + } + } + } + + laser_arr_chi(i, j) = chi; }); } } @@ -670,8 +703,7 @@ void MultiLaser::AdvanceSlice (const int islice, const Fields& fields, amrex::Real dt, int step) { - if (!m_use_laser) return; - if (!HasSlice(islice)) return; + if (!UseLaser(islice)) return; InterpolateChi(fields, m_laser_geom_3D); @@ -1123,8 +1155,6 @@ MultiLaser::AdvanceSliceFFT (const amrex::Real dt, int step) void MultiLaser::InitLaserSlice (const int islice, const int comp) { - if (!m_use_laser) return; - HIPACE_PROFILE("MultiLaser::InitLaserSlice()"); using namespace amrex::literals; @@ -1194,8 +1224,7 @@ void MultiLaser::InSituComputeDiags (int step, amrex::Real time, int islice, int max_step, amrex::Real max_time) { - if (!m_use_laser) return; - if (!HasSlice(islice)) return; + if (!UseLaser(islice)) return; if (!utils::doDiagnostics(m_insitu_period, step, max_step, time, max_time)) return; HIPACE_PROFILE("MultiLaser::InSituComputeDiags()"); diff --git a/src/utils/MultiBuffer.H b/src/utils/MultiBuffer.H index ccb8cb5092..dfb872fdd5 100644 --- a/src/utils/MultiBuffer.H +++ b/src/utils/MultiBuffer.H @@ -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 670491beed..0da539264c 100644 --- a/src/utils/MultiBuffer.cpp +++ b/src/utils/MultiBuffer.cpp @@ -62,8 +62,6 @@ void MultiBuffer::initialize (int nslices, MultiBeam& beams, MultiLaser& laser) m_nslices = nslices; m_nbeams = beams.get_nbeams(); - m_use_laser = laser.m_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, MultiLaser& laser) } } - 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,21 +775,22 @@ 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 && laser.HasSlice(slice)) { + 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 * laser.getSlices()[0].box().numPts() * sizeof(amrex::Real)); - memcpy_to_buffer(slice, get_buffer_offset(slice, offset_type::laser, beams, 0, 2), + memcpy_to_buffer(slice, get_buffer_offset(slice, offset_type::laser, beams, laser, 0, 2), laser.getSlices()[0].dataPtr(laser_comp_2_3), 2 * laser.getSlices()[0].box().numPts() * sizeof(amrex::Real)); } @@ -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,15 +855,15 @@ void MultiBuffer::unpack_data (int slice, MultiBeam& beams, MultiLaser& laser, i } } } - if (m_use_laser && laser.HasSlice(slice)) { + 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 * laser.getSlices()[0].box().numPts() * sizeof(amrex::Real)); - memcpy_from_buffer(slice, get_buffer_offset(slice, offset_type::laser, beams, 0, 2), + memcpy_from_buffer(slice, get_buffer_offset(slice, offset_type::laser, beams, laser, 0, 2), laser.getSlices()[0].dataPtr(laser_comp_2_3), 2 * laser.getSlices()[0].box().numPts() * sizeof(amrex::Real)); } From 4a2a92abf074afb7003ae9b286716f528e40cb41 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Sun, 28 Apr 2024 04:42:31 +0200 Subject: [PATCH 06/24] fix box --- src/fields/Fields.cpp | 6 ++++-- src/laser/MultiLaser.H | 4 ++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/fields/Fields.cpp b/src/fields/Fields.cpp index 3be383e7a2..ee588357b3 100644 --- a/src/fields/Fields.cpp +++ b/src/fields/Fields.cpp @@ -510,7 +510,8 @@ Fields::Copy (const int lev, const int i_slice, FieldDiagnosticData& fd, const amrex::Real dy = fd.m_geom_io.CellSize(1); if (fd.m_nfields > 0 && - calc_geom.smallEnd(2) <= i_slice && i_slice <= calc_geom.bigEnd(2)) { + calc_geom.Domain().smallEnd(2) <= i_slice && + i_slice <= calc_geom.Domain().bigEnd(2)) { auto slice_array = slice_func.array(mfi); amrex::Array4 diag_array = fd.m_F.array(); amrex::ParallelFor(diag_box, fd.m_nfields, @@ -524,7 +525,8 @@ Fields::Copy (const int lev, const int i_slice, FieldDiagnosticData& fd, } if (fd.m_do_laser && - laser_geom.smallEnd(2) <= i_slice && i_slice <= laser_geom.bigEnd(2)) { + laser_geom.Domain().smallEnd(2) <= i_slice && + i_slice <= laser_geom.Domain().bigEnd(2)) { auto laser_array = laser_func.array(mfi); amrex::Array4 diag_array_laser = fd.m_F_laser.array(); amrex::ParallelFor(diag_box, diff --git a/src/laser/MultiLaser.H b/src/laser/MultiLaser.H index f7a052b2ef..787f3134d1 100644 --- a/src/laser/MultiLaser.H +++ b/src/laser/MultiLaser.H @@ -224,8 +224,8 @@ public: * \param[in] islice slice index */ bool HasSlice (const int islice) const { - return GetLaserGeom().Domain.smallEnd(2) <= islice && - islice <= GetLaserGeom().Domain.bigEnd(2); + return GetLaserGeom().Domain().smallEnd(2) <= islice && + islice <= GetLaserGeom().Domain().bigEnd(2); } bool UseLaser () const { return m_use_laser; } From 440c2841e4716a1c603e79f90b65876c8131b0a9 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Sun, 28 Apr 2024 04:54:01 +0200 Subject: [PATCH 07/24] fix doc --- src/fields/Fields.cpp | 13 ++++--------- src/laser/MultiLaser.H | 2 +- src/laser/MultiLaser.cpp | 4 ++-- 3 files changed, 7 insertions(+), 12 deletions(-) diff --git a/src/fields/Fields.cpp b/src/fields/Fields.cpp index ee588357b3..24bb1c924f 100644 --- a/src/fields/Fields.cpp +++ b/src/fields/Fields.cpp @@ -426,8 +426,6 @@ Fields::Copy (const int lev, const int i_slice, FieldDiagnosticData& fd, constexpr int depos_order_z = 1; constexpr int depos_order_offset = depos_order_z / 2 + 1; - const amrex::Geometry& laser_geom = multi_laser.GetLaserGeom(); - const amrex::Real poff_calc_z = GetPosOffset(2, calc_geom, calc_geom.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()); @@ -490,7 +488,8 @@ Fields::Copy (const int lev, const int i_slice, FieldDiagnosticData& fd, auto& slice_mf = m_slices[lev]; auto slice_func = interpolated_field_xy{{slice_mf}, calc_geom}; auto& laser_mf = multi_laser.getSlices(); - auto laser_func = interpolated_field_xy{{laser_mf}, laser_geom}; + auto laser_func = interpolated_field_xy{{laser_mf}, multi_laser.GetLaserGeom()}; #ifdef AMREX_USE_GPU // This async copy happens on the same stream as the ParallelFor below, which uses the copied array. @@ -509,9 +508,7 @@ 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 && - calc_geom.Domain().smallEnd(2) <= i_slice && - i_slice <= calc_geom.Domain().bigEnd(2)) { + if (fd.m_nfields > 0) { auto slice_array = slice_func.array(mfi); amrex::Array4 diag_array = fd.m_F.array(); amrex::ParallelFor(diag_box, fd.m_nfields, @@ -524,9 +521,7 @@ Fields::Copy (const int lev, const int i_slice, FieldDiagnosticData& fd, }); } - if (fd.m_do_laser && - laser_geom.Domain().smallEnd(2) <= i_slice && - i_slice <= laser_geom.Domain().bigEnd(2)) { + if (fd.m_do_laser && multi_laser.UseLaser(i_slice)) { auto laser_array = laser_func.array(mfi); amrex::Array4 diag_array_laser = fd.m_F_laser.array(); amrex::ParallelFor(diag_box, diff --git a/src/laser/MultiLaser.H b/src/laser/MultiLaser.H index 787f3134d1..95aa72d972 100644 --- a/src/laser/MultiLaser.H +++ b/src/laser/MultiLaser.H @@ -157,7 +157,7 @@ public: * \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); diff --git a/src/laser/MultiLaser.cpp b/src/laser/MultiLaser.cpp index 34f77bc6b7..eb7f6adcf9 100644 --- a/src/laser/MultiLaser.cpp +++ b/src/laser/MultiLaser.cpp @@ -583,7 +583,7 @@ MultiLaser::UpdateLaserAabs (const int islice, const int current_N_level, Fields if (!m_use_laser) return; HIPACE_PROFILE("MultiLaser::UpdateLaserAabs()"); - constexpr interp_order = 1; + constexpr int interp_order = 1; if (!HasSlice(islice)) { for (int lev=0; lev Date: Sun, 28 Apr 2024 04:57:56 +0200 Subject: [PATCH 08/24] fix shape factor --- src/laser/MultiLaser.cpp | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/laser/MultiLaser.cpp b/src/laser/MultiLaser.cpp index eb7f6adcf9..96b8c39611 100644 --- a/src/laser/MultiLaser.cpp +++ b/src/laser/MultiLaser.cpp @@ -627,8 +627,10 @@ MultiLaser::UpdateLaserAabs (const int islice, const int current_N_level, Fields 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); + 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), @@ -685,8 +687,10 @@ MultiLaser::InterpolateChi (const Fields& fields, amrex::Geometry const& geom_fi 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); + 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) { chi += shape_x*shape_y*field_arr_chi(cell_x, cell_y); From 89d2ba0157a4c1d45b641bfb115127e195163686 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Sun, 28 Apr 2024 04:59:03 +0200 Subject: [PATCH 09/24] fix pos offset --- src/laser/MultiLaser.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/laser/MultiLaser.cpp b/src/laser/MultiLaser.cpp index 96b8c39611..b41a3c618c 100644 --- a/src/laser/MultiLaser.cpp +++ b/src/laser/MultiLaser.cpp @@ -677,8 +677,8 @@ MultiLaser::InterpolateChi (const Fields& fields, amrex::Geometry const& geom_fi amrex::ParallelFor(mfi.growntilebox(), [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept { - const amrex::Real x = i * dx_laser + m_laser_geom_3D; - const amrex::Real y = j * dy_laser + m_laser_geom_3D; + 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; From f35a5f1d6de7261a9652ef640c12d61fd6204c9c Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Sun, 28 Apr 2024 05:14:26 +0200 Subject: [PATCH 10/24] msvc fix --- src/laser/MultiLaser.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/laser/MultiLaser.cpp b/src/laser/MultiLaser.cpp index b41a3c618c..4c75ae5925 100644 --- a/src/laser/MultiLaser.cpp +++ b/src/laser/MultiLaser.cpp @@ -583,7 +583,7 @@ MultiLaser::UpdateLaserAabs (const int islice, const int current_N_level, Fields if (!m_use_laser) return; HIPACE_PROFILE("MultiLaser::UpdateLaserAabs()"); - constexpr int interp_order = 1; + static constexpr int interp_order = 1; if (!HasSlice(islice)) { for (int lev=0; lev Date: Sun, 28 Apr 2024 20:46:09 +0200 Subject: [PATCH 11/24] diagnostic rework for laser --- src/Hipace.H | 2 +- src/Hipace.cpp | 21 ++-- src/diagnostics/Diagnostic.H | 27 ++--- src/diagnostics/Diagnostic.cpp | 189 ++++++++++++++++++++---------- src/diagnostics/OpenPMDWriter.cpp | 13 +- src/fields/Fields.H | 8 +- src/fields/Fields.cpp | 37 +++--- 7 files changed, 174 insertions(+), 123 deletions(-) diff --git a/src/Hipace.H b/src/Hipace.H index 7e46ea4dc4..63a7f6324e 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 c6c2ea8506..1f31a7b63d 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -69,7 +69,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); @@ -197,9 +197,10 @@ Hipace::InitData () 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..f501fc9128 100644 --- a/src/diagnostics/Diagnostic.cpp +++ b/src/diagnostics/Diagnostic.cpp @@ -10,15 +10,24 @@ #include "utils/HipaceProfilerWrapper.H" #include -Diagnostic::Diagnostic (int nlev) +#include +#include +#include +#include + +Diagnostic::Diagnostic (int nlev, bool use_laser) { amrex::ParmParse ppd("diagnostic"); amrex::ParmParse pph("hipace"); amrex::Vector field_diag_names{}; - - for(int ilev = 0; ilev diag_name_to_default_geometry{}; + std::map geometry_name_to_geom_type{}; + std::map geometry_name_to_level{}; + std::map> geometry_name_to_output_comps{}; + std::stringstream all_comps_error_str{}; + for (int lev = 0; lev global_output_comps{}; + std::map is_global_comp_used{}; + queryWithParser(ppd, "field_data", global_output_comps); + for (auto& comp : global_output_comps) { + is_global_comp_used[comp] = false; + } + 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 (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 (do_laser) { - all_field_comps.push_back(fd.m_laser_io_name); + + 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 { + amrex::Abort("Unknown diagnostics base_geometry: '" + base_geom_name + "'!\n" + + all_comps_error_str.str()); } - if(fd.m_comps_output.empty()) { - fd.m_comps_output = all_field_comps; + amrex::Vector local_output_comps{}; + const bool use_local_comps = queryWithParser(pp, "field_data", local_output_comps); + + amrex::Vector use_comps{}; + use_comps = use_local_comps ? local_output_comps : global_output_comps; + + if (use_comps.empty()) { + fd.m_comps_output = geometry_name_to_output_comps[base_geom_name]; } else { - for (const std::string& comp_name : fd.m_comps_output) { + for (const std::string& comp_name : use_comps) { if (comp_name == "all" || comp_name == "All") { - fd.m_comps_output = all_field_comps; + fd.m_comps_output = geometry_name_to_output_comps[base_geom_name]; break; - } - if (comp_name == "none" || comp_name == "None") { + } else 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; + } else if (std::find(geometry_name_to_output_comps[base_geom_name].begin(), + geometry_name_to_output_comps[base_geom_name].end(), comp_name) != + geometry_name_to_output_comps[base_geom_name].end()) { + if (is_global_comp_used.count(comp_name) > 0) { + is_global_comp_used[comp_name] = true; } - amrex::Abort(error_str.str()); + fd.m_comps_output.push_back(comp_name); + } else if (use_local_comps) { + amrex::Abort("Unknown diagnostics field_data '" + comp_name + + "' in 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; - } else { - ++it; + fd.m_nfields = fd.m_comps_output.size(); + + 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_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]]; + 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()); } - fd.m_comps_output_idx.assign(local_comps_output_idx.begin(), local_comps_output_idx.end()); + } - if (m_field_data.size() != 1) { + 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 +268,26 @@ 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; + 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 +337,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 45132492a7..80395dcc18 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 24bb1c924f..7a5b8172cd 100644 --- a/src/fields/Fields.cpp +++ b/src/fields/Fields.cpp @@ -418,15 +418,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()); @@ -435,21 +435,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; @@ -476,17 +476,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}; + if (diag_box.isEmpty() || fd.m_nfields <= 0) return; + 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()}; @@ -508,7 +509,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, @@ -519,17 +521,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 && multi_laser.UseLaser(i_slice)) { + } 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) }; From b62746b79229077e79d3bd8df88152398e0e2927 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Sun, 28 Apr 2024 20:56:08 +0200 Subject: [PATCH 12/24] fix assert for all and none component --- src/diagnostics/Diagnostic.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/diagnostics/Diagnostic.cpp b/src/diagnostics/Diagnostic.cpp index f501fc9128..b7f8299b7d 100644 --- a/src/diagnostics/Diagnostic.cpp +++ b/src/diagnostics/Diagnostic.cpp @@ -189,9 +189,15 @@ Diagnostic::Initialize (int nlev, bool use_laser) { } else { for (const std::string& comp_name : use_comps) { if (comp_name == "all" || comp_name == "All") { + if (is_global_comp_used.count(comp_name) > 0) { + is_global_comp_used[comp_name] = true; + } fd.m_comps_output = geometry_name_to_output_comps[base_geom_name]; break; } else if (comp_name == "none" || comp_name == "None") { + if (is_global_comp_used.count(comp_name) > 0) { + is_global_comp_used[comp_name] = true; + } fd.m_comps_output.clear(); break; } else if (std::find(geometry_name_to_output_comps[base_geom_name].begin(), From ef35062703bb8a555f53cb881537bb2c55528633 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Sun, 28 Apr 2024 21:54:42 +0200 Subject: [PATCH 13/24] add option to remove field --- src/diagnostics/Diagnostic.cpp | 43 +++++++++++++++++----------------- 1 file changed, 22 insertions(+), 21 deletions(-) diff --git a/src/diagnostics/Diagnostic.cpp b/src/diagnostics/Diagnostic.cpp index b7f8299b7d..70d6688b04 100644 --- a/src/diagnostics/Diagnostic.cpp +++ b/src/diagnostics/Diagnostic.cpp @@ -12,6 +12,7 @@ #include #include +#include #include #include @@ -124,7 +125,7 @@ Diagnostic::Initialize (int nlev, bool use_laser) { std::map diag_name_to_default_geometry{}; std::map geometry_name_to_geom_type{}; std::map geometry_name_to_level{}; - std::map> geometry_name_to_output_comps{}; + std::map> geometry_name_to_output_comps{}; std::stringstream all_comps_error_str{}; for (int lev = 0; lev use_comps{}; use_comps = use_local_comps ? local_output_comps : global_output_comps; + std::set comps_set{}; + if (use_comps.empty()) { - fd.m_comps_output = geometry_name_to_output_comps[base_geom_name]; + comps_set.insert(geometry_name_to_output_comps[base_geom_name].begin(), + geometry_name_to_output_comps[base_geom_name].end()); } else { for (const std::string& comp_name : use_comps) { if (comp_name == "all" || comp_name == "All") { - if (is_global_comp_used.count(comp_name) > 0) { - is_global_comp_used[comp_name] = true; - } - fd.m_comps_output = geometry_name_to_output_comps[base_geom_name]; - break; + is_global_comp_used[comp_name] = true; + 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") { - if (is_global_comp_used.count(comp_name) > 0) { - is_global_comp_used[comp_name] = true; - } - fd.m_comps_output.clear(); - break; - } else if (std::find(geometry_name_to_output_comps[base_geom_name].begin(), - geometry_name_to_output_comps[base_geom_name].end(), comp_name) != - geometry_name_to_output_comps[base_geom_name].end()) { - if (is_global_comp_used.count(comp_name) > 0) { - is_global_comp_used[comp_name] = true; - } - fd.m_comps_output.push_back(comp_name); + is_global_comp_used[comp_name] = true; + comps_set.clear(); + } else if (geometry_name_to_output_comps[base_geom_name].count(comp_name) > 0) { + is_global_comp_used[comp_name] = true; + 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; + comps_set.erase(comp_name.substr(std::string("remove_").size(), comp_name.size())); } else if (use_local_comps) { amrex::Abort("Unknown diagnostics field_data '" + comp_name + "' in base_geometry '" + base_geom_name + "'!\n" + @@ -215,6 +215,7 @@ Diagnostic::Initialize (int nlev, bool use_laser) { } } + fd.m_comps_output.assign(comps_set.begin(), comps_set.end()); fd.m_nfields = fd.m_comps_output.size(); if (fd.m_base_geom_type == FieldDiagnosticData::geom_type::field) { From 61b16dc7a41c0f3678ea97bbca02bf924b30ea4b Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Mon, 29 Apr 2024 02:11:03 +0200 Subject: [PATCH 14/24] bug fix --- src/diagnostics/Diagnostic.cpp | 22 ++++++++++------------ src/fields/Fields.cpp | 2 +- src/laser/MultiLaser.cpp | 21 +++++++++++++-------- 3 files changed, 24 insertions(+), 21 deletions(-) diff --git a/src/diagnostics/Diagnostic.cpp b/src/diagnostics/Diagnostic.cpp index 70d6688b04..c0025b1632 100644 --- a/src/diagnostics/Diagnostic.cpp +++ b/src/diagnostics/Diagnostic.cpp @@ -136,7 +136,7 @@ Diagnostic::Initialize (int nlev, bool use_laser) { all_comps_error_str << "Available components in base_geometry '" << geom_name << "':\n "; for (const auto& [comp, idx] : Comps[WhichSlice::This]) { geometry_name_to_output_comps[geom_name].insert(comp); - all_comps_error_str << "'" << comp << "' "; + all_comps_error_str << comp << " "; } all_comps_error_str << "\n"; } @@ -149,16 +149,12 @@ Diagnostic::Initialize (int nlev, bool use_laser) { geometry_name_to_level.emplace(geom_name, 0); all_comps_error_str << "Available components in base_geometry '" << geom_name << "':\n "; geometry_name_to_output_comps[geom_name].insert(laser_io_name); - all_comps_error_str << "'" << laser_io_name << "'\n"; + all_comps_error_str << laser_io_name << "\n"; } - all_comps_error_str << "Additionally, 'all' and 'none' are supported as field_data\n"; + all_comps_error_str << "Additionally, 'all' and 'none' are supported as field_data\n" + << "Components can be removed after 'all' by using 'remove_'.\n"; - amrex::Vector global_output_comps{}; std::map is_global_comp_used{}; - queryWithParser(ppd, "field_data", global_output_comps); - for (auto& comp : global_output_comps) { - is_global_comp_used[comp] = false; - } for (auto& fd : m_field_data) { amrex::ParmParse pp(fd.m_diag_name); @@ -179,11 +175,11 @@ Diagnostic::Initialize (int nlev, bool use_laser) { all_comps_error_str.str()); } - amrex::Vector local_output_comps{}; - const bool use_local_comps = queryWithParser(pp, "field_data", local_output_comps); - amrex::Vector use_comps{}; - use_comps = use_local_comps ? local_output_comps : global_output_comps; + const bool use_local_comps = queryWithParser(pp, "field_data", use_comps); + if (!use_local_comps) { + queryWithParser(ppd, "field_data", use_comps); + } std::set comps_set{}; @@ -211,6 +207,8 @@ Diagnostic::Initialize (int nlev, bool use_laser) { amrex::Abort("Unknown diagnostics field_data '" + comp_name + "' in base_geometry '" + base_geom_name + "'!\n" + all_comps_error_str.str()); + } else { + is_global_comp_used.try_emplace(comp_name, false); } } } diff --git a/src/fields/Fields.cpp b/src/fields/Fields.cpp index 7a5b8172cd..43a2e6ffab 100644 --- a/src/fields/Fields.cpp +++ b/src/fields/Fields.cpp @@ -484,7 +484,7 @@ Fields::Copy (const int current_N_level, const int i_slice, FieldDiagnosticData& return; } } - if (diag_box.isEmpty() || fd.m_nfields <= 0) return; + if (diag_box.isEmpty()) return; auto& slice_mf = m_slices[fd.m_level]; auto slice_func = interpolated_field_xy{{slice_mf}, field_geom[fd.m_level]}; diff --git a/src/laser/MultiLaser.cpp b/src/laser/MultiLaser.cpp index 4c75ae5925..8a99c8cd3f 100644 --- a/src/laser/MultiLaser.cpp +++ b/src/laser/MultiLaser.cpp @@ -120,8 +120,11 @@ MultiLaser::MakeLaserGeometry (const amrex::Geometry& field_geom_3D) 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)}; - m_laser_geom_3D.define(domain_3D_laser, amrex::RealBox(patch_lo_laser, patch_hi_laser), - amrex::CoordSys::cartesian, {0, 0, 0}); + 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"); + + 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); @@ -581,6 +584,7 @@ MultiLaser::UpdateLaserAabs (const int islice, const int current_N_level, Fields amrex::Vector const& field_geom) { if (!m_use_laser) return; + if (!HasSlice(islice) && !HasSlice(islice + 1)) return; HIPACE_PROFILE("MultiLaser::UpdateLaserAabs()"); static constexpr int interp_order = 1; @@ -1169,9 +1173,10 @@ MultiLaser::InitLaserSlice (const int islice, const int comp) const amrex::Real k0 = 2._rt*MathConst::pi/m_lambda0; // Get grid properties - const auto plo = m_laser_geom_3D.ProbLoArray(); - amrex::Real const * const dx = m_laser_geom_3D.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()) @@ -1195,9 +1200,9 @@ MultiLaser::InitLaserSlice (const int islice, const int comp) 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)*y-std::sin(propagation_angle_yz)*z; const amrex::Real zp=std::sin(propagation_angle_yz)*y+std::cos(propagation_angle_yz)*z; From 1b3b435a7fd349622f1249e693480994128c6a6d Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Mon, 29 Apr 2024 04:51:39 +0200 Subject: [PATCH 15/24] add SetInitialChi --- src/Hipace.cpp | 2 + src/laser/MultiLaser.H | 6 +++ src/laser/MultiLaser.cpp | 102 ++++++++++++++++++++++++++++++--------- 3 files changed, 87 insertions(+), 23 deletions(-) diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 1f31a7b63d..29e8c28472 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -376,6 +376,8 @@ Hipace::Evolve () // Only reset plasma after receiving time step, to use proper density m_multi_plasma.InitData(m_slice_ba, m_slice_dm, m_slice_geom, m_3D_geom); + m_multi_laser.SetInitialChi(m_multi_plasma); + // deposit neutralizing background if (m_interpolate_neutralizing_background) { if (m_do_tiling) { diff --git a/src/laser/MultiLaser.H b/src/laser/MultiLaser.H index 95aa72d972..d1ea9d44f8 100644 --- a/src/laser/MultiLaser.H +++ b/src/laser/MultiLaser.H @@ -13,6 +13,7 @@ #include "Laser.H" #include "fields/Fields.H" #include "mg_solver/HpMultiGrid.H" +#include "particles/plasma/MultiPlasma.H" #include #include @@ -90,6 +91,7 @@ namespace WhichLaserSlice { np1jp2_r, np1jp2_i, chi, + chi_initial, N }; } @@ -164,6 +166,8 @@ public: void InterpolateChi (const Fields& fields, amrex::Geometry const& geom_field_lev0); + 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 @@ -254,6 +258,8 @@ private: 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 8a99c8cd3f..580ab3f930 100644 --- a/src/laser/MultiLaser.cpp +++ b/src/laser/MultiLaser.cpp @@ -62,6 +62,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); @@ -587,7 +589,6 @@ MultiLaser::UpdateLaserAabs (const int islice, const int current_N_level, Fields if (!HasSlice(islice) && !HasSlice(islice + 1)) return; HIPACE_PROFILE("MultiLaser::UpdateLaserAabs()"); - static constexpr int interp_order = 1; if (!HasSlice(islice)) { for (int lev=0; lev>{}, + {m_interp_order}, + mfi.growntilebox(), + [=] AMREX_GPU_DEVICE(int i, int j, int, auto interp_order) noexcept { using namespace WhichLaserSlice; const amrex::Real x = i * dx_field + poff_field_x; @@ -637,8 +641,8 @@ MultiLaser::UpdateLaserAabs (const int islice, const int current_N_level, Fields 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)); + aabs += shape_x*shape_y*abssq(laser_arr(cell_x, cell_y, n00j00_r), + laser_arr(cell_x, cell_y, n00j00_i)); } } } @@ -653,16 +657,51 @@ MultiLaser::UpdateLaserAabs (const int islice, const int current_N_level, Fields } } +void +MultiLaser::SetInitialChi (const MultiPlasma& multi_plasma) +{ + HIPACE_PROFILE("MultiLaser::SetInitialChi()"); + + for ( amrex::MFIter mfi(m_slices, DfltMfi); mfi.isValid(); ++mfi ){ + Array3 laser_arr_chi = m_slices.array(mfi, WhichLaserSlice::chi_initial); + + for (auto& plasma : multi_plasma.multi_plasma) { + + 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, auto interp_order) 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()"); - static constexpr int interp_order = 1; for ( amrex::MFIter mfi(m_slices, DfltMfi); mfi.isValid(); ++mfi ){ - - Array2 laser_arr_chi = m_slices.array(mfi, WhichLaserSlice::chi); - Array2 field_arr_chi = fields.getSlices(0).array(mfi, Comps[WhichSlice::This]["chi"]); + 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()); @@ -671,16 +710,31 @@ MultiLaser::InterpolateChi (const Fields& fields, amrex::Geometry const& geom_fi 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); - const int x_lo = fields.getSlices(0)[mfi].box().smallEnd(0); - const int x_hi = fields.getSlices(0)[mfi].box().bigEnd(0); - const int y_lo = fields.getSlices(0)[mfi].box().smallEnd(1); - const int y_hi = fields.getSlices(0)[mfi].box().bigEnd(1); + amrex::Box field_box = fields.getSlices(0)[mfi].box(); + field_box.grow(-2*Fields::m_slices_nguards); - amrex::ParallelFor(mfi.growntilebox(), - [=] AMREX_GPU_DEVICE(int i, int j, int) noexcept { + 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; + + 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; @@ -689,20 +743,22 @@ MultiLaser::InterpolateChi (const Fields& fields, amrex::Geometry const& geom_fi amrex::Real chi = 0; - 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 <= i && i <= x_hi && y_lo <= j && j <= y_hi) { + 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) { chi += shape_x*shape_y*field_arr_chi(cell_x, cell_y); } } + } else { + chi = laser_arr(i, j, WhichLaserSlice::chi_initial); } - laser_arr_chi(i, j) = chi; + laser_arr(i, j, WhichLaserSlice::chi) = chi; }); } } From 72da836707888ad35b08e185ff4f35cde59adca4 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Mon, 29 Apr 2024 04:55:28 +0200 Subject: [PATCH 16/24] try fix includes --- src/laser/MultiLaser.H | 4 ++-- src/laser/MultiLaser.cpp | 2 ++ 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/laser/MultiLaser.H b/src/laser/MultiLaser.H index d1ea9d44f8..b08bc529dc 100644 --- a/src/laser/MultiLaser.H +++ b/src/laser/MultiLaser.H @@ -11,9 +11,7 @@ #define MULTILASER_H_ #include "Laser.H" -#include "fields/Fields.H" #include "mg_solver/HpMultiGrid.H" -#include "particles/plasma/MultiPlasma.H" #include #include @@ -98,6 +96,8 @@ namespace WhichLaserSlice { class Fields; +class MultiPlasma; + class MultiLaser { diff --git a/src/laser/MultiLaser.cpp b/src/laser/MultiLaser.cpp index 580ab3f930..914c724c47 100644 --- a/src/laser/MultiLaser.cpp +++ b/src/laser/MultiLaser.cpp @@ -9,7 +9,9 @@ #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" From 62b07a74e34343848027ecbee9c9b81085225843 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Mon, 29 Apr 2024 04:58:54 +0200 Subject: [PATCH 17/24] fix Array2 --- src/laser/MultiLaser.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/laser/MultiLaser.cpp b/src/laser/MultiLaser.cpp index 914c724c47..6c602b4bcd 100644 --- a/src/laser/MultiLaser.cpp +++ b/src/laser/MultiLaser.cpp @@ -665,7 +665,7 @@ MultiLaser::SetInitialChi (const MultiPlasma& multi_plasma) HIPACE_PROFILE("MultiLaser::SetInitialChi()"); for ( amrex::MFIter mfi(m_slices, DfltMfi); mfi.isValid(); ++mfi ){ - Array3 laser_arr_chi = m_slices.array(mfi, WhichLaserSlice::chi_initial); + Array2 laser_arr_chi = m_slices.array(mfi, WhichLaserSlice::chi_initial); for (auto& plasma : multi_plasma.multi_plasma) { From 813c1b8c89ff056241f0d7440ce821b300a7760d Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Mon, 29 Apr 2024 05:01:03 +0200 Subject: [PATCH 18/24] fix m_all_plasmas --- src/laser/MultiLaser.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/laser/MultiLaser.cpp b/src/laser/MultiLaser.cpp index 6c602b4bcd..00e70c33c7 100644 --- a/src/laser/MultiLaser.cpp +++ b/src/laser/MultiLaser.cpp @@ -667,7 +667,7 @@ MultiLaser::SetInitialChi (const MultiPlasma& multi_plasma) for ( amrex::MFIter mfi(m_slices, DfltMfi); mfi.isValid(); ++mfi ){ Array2 laser_arr_chi = m_slices.array(mfi, WhichLaserSlice::chi_initial); - for (auto& plasma : multi_plasma.multi_plasma) { + 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; From cc8af73d9206f9f410e7b09cf57ad9619a8cb36f Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Thu, 2 May 2024 17:12:17 +0200 Subject: [PATCH 19/24] remove unused variable --- src/laser/MultiLaser.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/laser/MultiLaser.cpp b/src/laser/MultiLaser.cpp index 00e70c33c7..0ef5e36831 100644 --- a/src/laser/MultiLaser.cpp +++ b/src/laser/MultiLaser.cpp @@ -685,7 +685,7 @@ MultiLaser::SetInitialChi (const MultiPlasma& multi_plasma) const amrex::Real dy_laser = m_laser_geom_3D.CellSize(1); amrex::ParallelFor(mfi.growntilebox(), - [=] AMREX_GPU_DEVICE(int i, int j, int, auto interp_order) noexcept { + [=] 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; From fc0d3ab3ab9200dcc28ca938add9a3c15a5613c6 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Thu, 2 May 2024 19:49:35 +0200 Subject: [PATCH 20/24] add code comments --- src/Hipace.cpp | 2 +- src/diagnostics/Diagnostic.cpp | 76 ++++++++++++++++++++++------------ src/laser/MultiLaser.H | 15 ++++++- src/laser/MultiLaser.cpp | 23 ++++++++-- 4 files changed, 84 insertions(+), 32 deletions(-) diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 29e8c28472..cc1ac41576 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -520,7 +520,7 @@ Hipace::SolveOneSlice (int islice, int step) // Advance laser slice by 1 step using chi // no MR for laser - m_multi_laser.AdvanceSlice(islice, m_fields, m_dt, step); + m_multi_laser.AdvanceSlice(islice, m_fields, m_dt, step, m_3D_geom[0]); if (islice-1 >= m_3D_geom[0].Domain().smallEnd(2)) { m_multi_buffer.get_data(islice-1, m_multi_beam, m_multi_laser, WhichBeamSlice::Next); diff --git a/src/diagnostics/Diagnostic.cpp b/src/diagnostics/Diagnostic.cpp index c0025b1632..f9f0783858 100644 --- a/src/diagnostics/Diagnostic.cpp +++ b/src/diagnostics/Diagnostic.cpp @@ -21,6 +21,7 @@ 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 lev = 0; lev 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) { @@ -181,41 +190,52 @@ Diagnostic::Initialize (int nlev, bool use_laser) { queryWithParser(ppd, "field_data", use_comps); } + // set to store all used components to avoid duplicates std::set comps_set{}; if (use_comps.empty()) { - comps_set.insert(geometry_name_to_output_comps[base_geom_name].begin(), - geometry_name_to_output_comps[base_geom_name].end()); - } else { - for (const std::string& comp_name : use_comps) { - if (comp_name == "all" || comp_name == "All") { - is_global_comp_used[comp_name] = true; - 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; - comps_set.clear(); - } else if (geometry_name_to_output_comps[base_geom_name].count(comp_name) > 0) { - is_global_comp_used[comp_name] = true; - 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; - comps_set.erase(comp_name.substr(std::string("remove_").size(), comp_name.size())); - } else if (use_local_comps) { - amrex::Abort("Unknown diagnostics field_data '" + comp_name + - "' in base_geometry '" + base_geom_name + "'!\n" + - all_comps_error_str.str()); - } else { - is_global_comp_used.try_emplace(comp_name, false); - } + // 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 { + // if field_data was specified through diagnostic, + // check later that all components are at least used by on 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(); + // 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) { @@ -225,6 +245,7 @@ Diagnostic::Initialize (int nlev, bool use_laser) { } } + // check that all components are at least used by on of the diagnostics for (auto& [key, val] : is_global_comp_used) { if (!val) { amrex::Abort("Unknown or unused component in diagnostic.field_data.\n'" + @@ -233,6 +254,8 @@ Diagnostic::Initialize (int nlev, bool use_laser) { } } + // if there are multiple diagnositc objects with the same m_base_geom_type (colliding component + // names), append the name of the diagnositc 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; @@ -283,6 +306,7 @@ Diagnostic::ResizeFDiagFAB (amrex::Vector& field_geom, amrex::Geometry geom; + // choose the geometry of the diagnostic switch (fd.m_base_geom_type) { case FieldDiagnosticData::geom_type::field: geom = field_geom[fd.m_level]; diff --git a/src/laser/MultiLaser.H b/src/laser/MultiLaser.H index b08bc529dc..bef6149abd 100644 --- a/src/laser/MultiLaser.H +++ b/src/laser/MultiLaser.H @@ -164,8 +164,15 @@ public: 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 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. @@ -173,8 +180,10 @@ public: * \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 int islice, 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. @@ -232,8 +241,12 @@ public: 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: diff --git a/src/laser/MultiLaser.cpp b/src/laser/MultiLaser.cpp index 0ef5e36831..ce69f72f16 100644 --- a/src/laser/MultiLaser.cpp +++ b/src/laser/MultiLaser.cpp @@ -65,7 +65,7 @@ MultiLaser::ReadParameters () 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); + 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); @@ -93,6 +93,7 @@ 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 { @@ -104,10 +105,12 @@ MultiLaser::MakeLaserGeometry (const amrex::Geometry& field_geom_3D) 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), @@ -121,6 +124,7 @@ MultiLaser::MakeLaserGeometry (const amrex::Geometry& field_geom_3D) 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)}; @@ -128,6 +132,7 @@ MultiLaser::MakeLaserGeometry (const amrex::Geometry& field_geom_3D) 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; @@ -593,6 +598,8 @@ MultiLaser::UpdateLaserAabs (const int islice, const int current_N_level, Fields 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_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(); @@ -689,7 +700,7 @@ MultiLaser::SetInitialChi (const MultiPlasma& multi_plasma) 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; + laser_arr_chi(i, j) += density_func(x, y, c_t) * chi_factor; }); } } @@ -727,6 +738,7 @@ MultiLaser::InterpolateChi (const Fields& fields, amrex::Geometry const& geom_fi 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); @@ -746,6 +758,7 @@ MultiLaser::InterpolateChi (const Fields& fields, amrex::Geometry const& geom_fi 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] = @@ -757,6 +770,7 @@ MultiLaser::InterpolateChi (const Fields& fields, amrex::Geometry const& geom_fi } } } else { + // get initial chi outside the fields box chi = laser_arr(i, j, WhichLaserSlice::chi_initial); } @@ -766,12 +780,13 @@ MultiLaser::InterpolateChi (const Fields& fields, amrex::Geometry const& geom_fi } void -MultiLaser::AdvanceSlice (const int islice, const Fields& fields, amrex::Real dt, int step) +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, m_laser_geom_3D); + InterpolateChi(fields, geom_field_lev0); if (m_solver_type == "multigrid") { AdvanceSliceMG(dt, step); From 50478a928af6120960d3f3accc03c2597ebd02a0 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Thu, 2 May 2024 20:02:18 +0200 Subject: [PATCH 21/24] fix format --- src/diagnostics/Diagnostic.cpp | 8 ++++---- src/laser/MultiLaser.cpp | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/diagnostics/Diagnostic.cpp b/src/diagnostics/Diagnostic.cpp index f9f0783858..8a9e3a0958 100644 --- a/src/diagnostics/Diagnostic.cpp +++ b/src/diagnostics/Diagnostic.cpp @@ -123,7 +123,7 @@ void Diagnostic::Initialize (int nlev, bool use_laser) { amrex::ParmParse ppd("diagnostic"); - // for each diagnostic object, choose a geometry and assign field_data + // for each diagnostic object, choose a geometry and assign field_data // for the default diagnostics, what is the default geometry std::map diag_name_to_default_geometry{}; @@ -220,13 +220,13 @@ Diagnostic::Initialize (int nlev, bool use_laser) { // 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 , + // 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 { - // if field_data was specified through diagnostic, + // if field_data was specified through diagnostic, // check later that all components are at least used by on of the diagnostics is_global_comp_used.try_emplace(comp_name, false); } @@ -254,7 +254,7 @@ Diagnostic::Initialize (int nlev, bool use_laser) { } } - // if there are multiple diagnositc objects with the same m_base_geom_type (colliding component + // if there are multiple diagnositc objects with the same m_base_geom_type (colliding component // names), append the name of the diagnositc 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) { diff --git a/src/laser/MultiLaser.cpp b/src/laser/MultiLaser.cpp index 849407741d..979c2a91ab 100644 --- a/src/laser/MultiLaser.cpp +++ b/src/laser/MultiLaser.cpp @@ -93,7 +93,7 @@ MultiLaser::MakeLaserGeometry (const amrex::Geometry& field_geom_3D) if (!m_use_laser) return; amrex::ParmParse pp("lasers"); - // use field_geom_3D as the default + // 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 { @@ -132,7 +132,7 @@ MultiLaser::MakeLaserGeometry (const amrex::Geometry& field_geom_3D) AMREX_ALWAYS_ASSERT_WITH_MESSAGE(real_box.volume() > 0., "Laser box must have positive volume"); - // make the geometry, slice box and ba and dm + // 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; From f6c4d5cc65cd36bd910b4559e2c2886d18d383df Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Thu, 2 May 2024 20:32:07 +0200 Subject: [PATCH 22/24] add doc --- docs/source/run/parameters.rst | 36 +++++++++++++++++++++++++++------- 1 file changed, 29 insertions(+), 7 deletions(-) diff --git a/docs/source/run/parameters.rst b/docs/source/run/parameters.rst index 79dd078913..01f9049ac1 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 --------- @@ -789,6 +799,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. @@ -882,11 +896,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`) + On 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``. @@ -919,8 +937,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. From 82f12f04045d993bffb8c6f0621d075b75187b35 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Mon, 27 May 2024 15:18:19 +0200 Subject: [PATCH 23/24] add suggestions --- docs/source/run/parameters.rst | 2 +- src/diagnostics/Diagnostic.cpp | 6 +++--- src/laser/MultiLaser.cpp | 2 ++ 3 files changed, 6 insertions(+), 4 deletions(-) diff --git a/docs/source/run/parameters.rst b/docs/source/run/parameters.rst index f8db98acd2..9bd1d99ce2 100644 --- a/docs/source/run/parameters.rst +++ b/docs/source/run/parameters.rst @@ -916,7 +916,7 @@ Field diagnostics a subset of ``lev0 lev1 lev2 laser_diag``. * `` or diagnostic.base_geometry`` (`string`) optional (default `level_0`) - On which geometry the diagnostics should be based on. + 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 diff --git a/src/diagnostics/Diagnostic.cpp b/src/diagnostics/Diagnostic.cpp index 8a9e3a0958..24b1021618 100644 --- a/src/diagnostics/Diagnostic.cpp +++ b/src/diagnostics/Diagnostic.cpp @@ -227,7 +227,7 @@ Diagnostic::Initialize (int nlev, bool use_laser) { all_comps_error_str.str()); } else { // if field_data was specified through diagnostic, - // check later that all components are at least used by on of the diagnostics + // check later that all components are at least used by one of the diagnostics is_global_comp_used.try_emplace(comp_name, false); } } @@ -254,8 +254,8 @@ Diagnostic::Initialize (int nlev, bool use_laser) { } } - // if there are multiple diagnositc objects with the same m_base_geom_type (colliding component - // names), append the name of the diagnositc object to the component name in the output + // 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; diff --git a/src/laser/MultiLaser.cpp b/src/laser/MultiLaser.cpp index 79e271137f..ac84172959 100644 --- a/src/laser/MultiLaser.cpp +++ b/src/laser/MultiLaser.cpp @@ -667,6 +667,8 @@ MultiLaser::InterpolateChi (const Fields& fields, amrex::Geometry const& geom_fi 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; From 798c8deecdb1cb294d0f6070a00709a23ce9f216 Mon Sep 17 00:00:00 2001 From: AlexanderSinn Date: Mon, 27 May 2024 15:25:32 +0200 Subject: [PATCH 24/24] add suggestions 2 --- src/diagnostics/Diagnostic.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/diagnostics/Diagnostic.cpp b/src/diagnostics/Diagnostic.cpp index 24b1021618..963bb1b381 100644 --- a/src/diagnostics/Diagnostic.cpp +++ b/src/diagnostics/Diagnostic.cpp @@ -245,7 +245,7 @@ Diagnostic::Initialize (int nlev, bool use_laser) { } } - // check that all components are at least used by on of the diagnostics + // 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'" +