From 4177aac284510122ff90d835a47c7511c64a3ecf Mon Sep 17 00:00:00 2001 From: AlexanderSinn <64009254+AlexanderSinn@users.noreply.github.com> Date: Thu, 6 Apr 2023 17:42:28 +0200 Subject: [PATCH] Remove zeta refinement (#911) * remove zeta refinement * add assert --- src/Hipace.H | 4 +- src/Hipace.cpp | 96 ++++++++++++++++++------------------------- src/fields/Fields.H | 3 +- src/fields/Fields.cpp | 82 ++++-------------------------------- 4 files changed, 51 insertions(+), 134 deletions(-) diff --git a/src/Hipace.H b/src/Hipace.H index 0a01a138d3..7f03bdd96d 100644 --- a/src/Hipace.H +++ b/src/Hipace.H @@ -170,12 +170,12 @@ public: /** \brief Full evolve on 1 slice * - * \param[in] islice_coarse slice number + * \param[in] islice slice number * \param[in] ibox index of the current box to be calculated * \param[in] step current time step * \param[in] bins an amrex::DenseBins object that orders particles by slice */ - void SolveOneSlice (int islice_coarse, const int ibox, int step, + void SolveOneSlice (int islice, const int ibox, int step, const amrex::Vector>& bins); /** \brief Full evolve on 1 subslice with explicit solver diff --git a/src/Hipace.cpp b/src/Hipace.cpp index a1575f5b14..588d3b7ff7 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -186,6 +186,8 @@ Hipace::Hipace () : for (int idim=0; idim>& bins) { HIPACE_PROFILE("Hipace::SolveOneSlice()"); - m_multi_beam.InSituComputeDiags(step, islice_coarse, bins[0], + m_multi_beam.InSituComputeDiags(step, islice, bins[0], boxArray(0)[ibox].smallEnd(Direction::z), m_box_sorters, ibox); // Get this laser slice from the 3D array - m_multi_laser.Copy(islice_coarse, false); + m_multi_laser.Copy(islice, false); for (int lev = 0; lev <= finestLevel(); ++lev) { @@ -583,7 +585,7 @@ Hipace::SolveOneSlice (int islice_coarse, const int ibox, int step, // use geometry of coarse grid to determine whether slice is to be solved const amrex::Real* problo = Geom(0).ProbLo(); const amrex::Real* dx = Geom(0).CellSize(); - amrex::Real pos = (islice_coarse+0.5)*dx[2]+problo[2]; + amrex::Real pos = (islice+0.5)*dx[2]+problo[2]; if (pos < patch_lo[2] || pos > patch_hi[2]) continue; } @@ -593,58 +595,49 @@ Hipace::SolveOneSlice (int islice_coarse, const int ibox, int step, const amrex::Box& bx = boxArray(lev)[ibox]; - const int nsubslice = GetRefRatio(lev)[Direction::z]; + // calculate correct slice for refined level + const int islice_local = islice - bx.smallEnd(Direction::z); - for (int isubslice = nsubslice-1; isubslice >= 0; --isubslice) { + // reorder plasma before TileSort + m_multi_plasma.ReorderParticles(islice); - // calculate correct slice for refined level - const int islice = nsubslice*islice_coarse + isubslice; - const int islice_local = islice - bx.smallEnd(Direction::z); - - // reorder plasma before TileSort - m_multi_plasma.ReorderParticles(islice); - - if (m_explicit) { - ExplicitSolveOneSubSlice(lev, step, ibox, bx, islice, islice_local, bins[lev]); - } else { - PredictorCorrectorSolveOneSubSlice(lev, step, ibox, bx, islice, - islice_local, bins[lev]); - } + if (m_explicit) { + ExplicitSolveOneSubSlice(lev, step, ibox, bx, islice, islice_local, bins[lev]); + } else { + PredictorCorrectorSolveOneSubSlice(lev, step, ibox, bx, islice, + islice_local, bins[lev]); + } - FillDiagnostics(lev, islice); + FillDiagnostics(lev, islice); - m_multi_plasma.DoFieldIonization(lev, geom[lev], m_fields); + m_multi_plasma.DoFieldIonization(lev, geom[lev], m_fields); - if (m_multi_plasma.IonizationOn() && m_do_tiling) m_multi_plasma.TileSort(bx, geom[lev]); + if (m_multi_plasma.IonizationOn() && m_do_tiling) m_multi_plasma.TileSort(bx, geom[lev]); - // Push plasma particles - m_multi_plasma.AdvanceParticles(m_fields, m_multi_laser, geom[lev], false, lev); + // Push plasma particles + m_multi_plasma.AdvanceParticles(m_fields, m_multi_laser, geom[lev], false, lev); - m_multi_plasma.doCoulombCollision(lev, bx, geom[lev]); + m_multi_plasma.doCoulombCollision(lev, bx, geom[lev]); - // Push beam particles - m_multi_beam.AdvanceBeamParticlesSlice(m_fields, geom[lev], lev, islice_local, bx, - bins[lev], m_box_sorters, ibox); + // Push beam particles + m_multi_beam.AdvanceBeamParticlesSlice(m_fields, geom[lev], lev, islice_local, bx, + bins[lev], m_box_sorters, ibox); - if (lev == 0) { - // Advance laser slice by 1 step and store result to 3D array - m_multi_laser.AdvanceSlice(m_fields, Geom(0), m_dt, step); - m_multi_laser.Copy(islice_coarse, true); - } - - if (lev != 0) { - // shift slices of level 1 - m_fields.ShiftSlices(lev, islice_coarse, Geom(0), patch_lo[2], patch_hi[2]); - } - } // end for (int isubslice = nsubslice-1; isubslice >= 0; --isubslice) + if (lev == 0) { + // Advance laser slice by 1 step and store result to 3D array + m_multi_laser.AdvanceSlice(m_fields, Geom(0), m_dt, step); + m_multi_laser.Copy(islice, true); + } // After this, the parallel context is the full 3D communicator again amrex::ParallelContext::pop(); } // end for (int lev = 0; lev <= finestLevel(); ++lev) - // shift level 0 - m_fields.ShiftSlices(0, islice_coarse, Geom(0), patch_lo[2], patch_hi[2]); + // shift all levels + for (int lev = 0; lev <= finestLevel(); ++lev) { + m_fields.ShiftSlices(lev, islice, Geom(0), patch_lo[2], patch_hi[2]); + } } @@ -681,7 +674,7 @@ Hipace::ExplicitSolveOneSubSlice (const int lev, const int step, const int ibox, FillBoundaryChargeCurrents(lev); - m_fields.SolvePoissonExmByAndEypBx(Geom(), m_comm_xy, lev, islice); + m_fields.SolvePoissonExmByAndEypBx(Geom(), lev, islice); m_fields.SolvePoissonEz(Geom(), lev, islice); m_fields.SolvePoissonBz(Geom(), lev, islice); @@ -731,7 +724,7 @@ Hipace::PredictorCorrectorSolveOneSubSlice (const int lev, const int step, const FillBoundaryChargeCurrents(lev); if (!m_do_beam_jz_minus_rho) { - m_fields.SolvePoissonExmByAndEypBx(Geom(), m_comm_xy, lev, islice); + m_fields.SolvePoissonExmByAndEypBx(Geom(), lev, islice); } // deposit jx jy jz and maybe rho on This slice @@ -740,7 +733,7 @@ Hipace::PredictorCorrectorSolveOneSubSlice (const int lev, const int step, const m_do_beam_jz_minus_rho, WhichSlice::This); if (m_do_beam_jz_minus_rho) { - m_fields.SolvePoissonExmByAndEypBx(Geom(), m_comm_xy, lev, islice); + m_fields.SolvePoissonExmByAndEypBx(Geom(), lev, islice); } // deposit grid current into jz_beam @@ -842,7 +835,6 @@ void Hipace::ExplicitMGSolveBxBy (const int lev, const int which_slice, const int islice) { HIPACE_PROFILE("Hipace::ExplicitMGSolveBxBy()"); - amrex::ParallelContext::push(m_comm_xy); // always get chi from WhichSlice::This const int which_slice_chi = WhichSlice::This; @@ -939,7 +931,6 @@ Hipace::ExplicitMGSolveBxBy (const int lev, const int which_slice, const int isl m_hpmg[lev]->solve1(BxBy[0], SySx[0], Mult[0], m_MG_tolerance_rel, m_MG_tolerance_abs, max_iters, m_MG_verbose); } - amrex::ParallelContext::pop(); } void @@ -960,11 +951,9 @@ Hipace::PredictorCorrectorLoopToSolveBxBy (const int islice_local, const int lev m_fields.InitialBfieldGuess(relative_Bfield_error, m_predcorr_B_error_tolerance, lev); if (!m_fields.m_extended_solve) { - amrex::ParallelContext::push(m_comm_xy); // exchange ExmBy EypBx Ez Bx By Bz m_fields.FillBoundary(Geom(lev).periodicity(), lev, WhichSlice::This, "ExmBy", "EypBx", "Ez", "Bx", "By", "Bz", "jx", "jy", "jz", "rho", "Psi"); - amrex::ParallelContext::pop(); } /* creating temporary Bx and By arrays for the current and previous iteration */ @@ -1012,11 +1001,9 @@ Hipace::PredictorCorrectorLoopToSolveBxBy (const int islice_local, const int lev m_box_sorters, ibox, m_do_beam_jx_jy_deposition, false, false, WhichSlice::Next); if (!m_fields.m_extended_solve) { - amrex::ParallelContext::push(m_comm_xy); // need to exchange jx jy jx_beam jy_beam m_fields.FillBoundary(Geom(lev).periodicity(), lev, WhichSlice::Next, "jx", "jy"); - amrex::ParallelContext::pop(); } /* Calculate Bx and By */ @@ -1043,11 +1030,9 @@ Hipace::PredictorCorrectorLoopToSolveBxBy (const int islice_local, const int lev m_fields.setVal(0., lev, WhichSlice::Next, "jx", "jy"); if (!m_fields.m_extended_solve) { - amrex::ParallelContext::push(m_comm_xy); // exchange Bx By m_fields.FillBoundary(Geom(lev).periodicity(), lev, WhichSlice::This, "ExmBy", "EypBx", "Ez", "Bx", "By", "Bz", "jx", "jy", "jz", "rho", "Psi"); - amrex::ParallelContext::pop(); } // Shift relative_Bfield_error values @@ -1496,16 +1481,15 @@ Hipace::ResizeFDiagFAB (const int it, const int step) // boxArray(1) is not correct in z direction. We need to manually enforce a // parent/child relationship between lev_0 and lev_1 boxes in z const amrex::Box& bx_lev0 = boxArray(0)[it]; - const int ref_ratio_z = GetRefRatio(lev)[Direction::z]; // This seems to be required for some reason - domain.setBig(Direction::z, domain.bigEnd(Direction::z) - ref_ratio_z); + domain.setBig(Direction::z, domain.bigEnd(Direction::z) - 1); // Ensuring the IO boxes on level 1 are aligned with the boxes on level 0 local_box.setSmall(Direction::z, amrex::max(domain.smallEnd(Direction::z), - ref_ratio_z*bx_lev0.smallEnd(Direction::z))); + bx_lev0.smallEnd(Direction::z))); local_box.setBig (Direction::z, amrex::min(domain.bigEnd(Direction::z), - ref_ratio_z*bx_lev0.bigEnd(Direction::z)+(ref_ratio_z-1))); + bx_lev0.bigEnd(Direction::z))); } m_diags.ResizeFDiagFAB(local_box, domain, lev, Geom(lev), diff --git a/src/fields/Fields.H b/src/fields/Fields.H index 5ce7a35eff..a5af8666f5 100644 --- a/src/fields/Fields.H +++ b/src/fields/Fields.H @@ -204,12 +204,11 @@ public: * ExmBy and EypBx are solved in the same function because both rely on Psi. * * \param[in] geom Geometry - * \param[in] m_comm_xy transverse communicator on the slice * \param[in] lev current level * \param[in] islice longitudinal slice */ void SolvePoissonExmByAndEypBx (amrex::Vector const& geom, - const MPI_Comm& m_comm_xy, const int lev, const int islice); + const int lev, const int islice); /** \brief Compute Ez on the slice container from J by solving a Poisson equation * * \param[in] geom Geometry diff --git a/src/fields/Fields.cpp b/src/fields/Fields.cpp index d5737494fd..caeee69b40 100644 --- a/src/fields/Fields.cpp +++ b/src/fields/Fields.cpp @@ -63,7 +63,6 @@ Fields::AllocData ( m_explicit = Hipace::GetInstance().m_explicit; m_any_neutral_background = Hipace::GetInstance().m_multi_plasma.AnySpeciesNeutralizeBackground(); - const bool mesh_refinement = Hipace::m_do_MR; const bool any_salame = Hipace::GetInstance().m_multi_beam.AnySpeciesSalame(); if (m_explicit) { @@ -85,10 +84,6 @@ Fields::AllocData ( "jx_beam", "jy_beam", "jz_beam", "rho_beam", "jx", "jy", "jz", "rho"); isl = WhichSlice::Previous1; - if (mesh_refinement) { - // for interpolating boundary conditions to level 1 - Comps[isl].multi_emplace(N_Comps, "Ez", "Bx", "By", "Bz", "Psi", "rho"); - } Comps[isl].multi_emplace(N_Comps, "jx_beam", "jy_beam"); isl = WhichSlice::Previous2; @@ -122,10 +117,6 @@ Fields::AllocData ( } isl = WhichSlice::Previous1; - if (mesh_refinement) { - // for interpolating boundary conditions to level 1 - Comps[isl].multi_emplace(N_Comps, "Ez", "Bz", "Psi", "rho"); - } Comps[isl].multi_emplace(N_Comps, "Bx", "By", "jx", "jy"); isl = WhichSlice::Previous2; @@ -287,45 +278,13 @@ struct interpolated_field_xy { // use .array(mfi) like with amrex::MultiFab auto array (amrex::MFIter& mfi) const { - auto mfab_array = mfab.array(mfi); + auto mfab_array = to_array2(mfab.array(mfi)); return interpolated_field_xy_inner{ mfab_array, 1._rt/geom.CellSize(0), 1._rt/geom.CellSize(1), GetPosOffset(0, geom, geom.Domain()), GetPosOffset(1, geom, geom.Domain())}; } }; -/** \brief inner version of interpolated_field_z */ -struct interpolated_field_z_inner { - // captured variables for GPU - Array2 arr_this; - Array2 arr_prev; - amrex::Real rel_z; - - // linear longitudinal field interpolation - AMREX_GPU_DEVICE amrex::Real operator() (int i, int j) const noexcept { - return (1._rt-rel_z)*arr_this(i, j) + rel_z*arr_prev(i, j); - } -}; - -/** \brief linear longitudinal field interpolation */ -struct interpolated_field_z { - // use brace initialization as constructor - amrex::MultiFab mfab_this; // field to interpolate on this slice - amrex::MultiFab mfab_prev; // field to interpolate on previous slice - amrex::Real rel_z; // mixing factor between f_view_this and f_view_prev for z interpolation - - // use .array(mfi) like with amrex::MultiFab - interpolated_field_z_inner array (amrex::MFIter& mfi) const { - return interpolated_field_z_inner{ - mfab_this.array(mfi), mfab_prev.array(mfi), rel_z}; - } -}; - -/** \brief interpolate field in interp_order_xy order transversely and - * first order (linear) longitudinally */ -template -using interpolated_field_xyz = interpolated_field_xy; - /** \brief inner version of guarded_field_xy */ struct guarded_field_xy_inner { // captured variables for GPU @@ -522,26 +481,15 @@ Fields::ShiftSlices (int lev, int islice, amrex::Geometry geom, amrex::Real patc } const bool explicit_solve = Hipace::GetInstance().m_explicit; - const bool mesh_refinement = Hipace::m_do_MR; // only shift the slices that are allocated if (explicit_solve) { - if (mesh_refinement) { - shift(lev, WhichSlice::Previous1, WhichSlice::This, - "Ez", "Bx", "By", "Bz", "rho", "Psi", "jx_beam", "jy_beam"); - } else { - shift(lev, WhichSlice::Previous1, WhichSlice::This, "jx_beam", "jy_beam"); - } + shift(lev, WhichSlice::Previous1, WhichSlice::This, "jx_beam", "jy_beam"); duplicate(lev, WhichSlice::This, {"jx_beam", "jy_beam", "jx" , "jy" }, WhichSlice::Next, {"jx_beam", "jy_beam", "jx_beam", "jy_beam"}); } else { shift(lev, WhichSlice::Previous2, WhichSlice::Previous1, "Bx", "By"); - if (mesh_refinement) { - shift(lev, WhichSlice::Previous1, WhichSlice::This, - "Ez", "Bx", "By", "Bz", "jx", "jy", "rho", "Psi"); - } else { - shift(lev, WhichSlice::Previous1, WhichSlice::This, "Bx", "By", "jx", "jy"); - } + shift(lev, WhichSlice::Previous1, WhichSlice::This, "Bx", "By", "jx", "jy"); } } @@ -686,14 +634,8 @@ Fields::SetBoundaryCondition (amrex::Vector const& geom, const // Fine level: interpolate solution from coarser level to get Dirichlet boundary conditions constexpr int interp_order = 2; - const amrex::Real ref_ratio_z = Hipace::GetRefRatio(lev)[2]; - const amrex::Real islice_coarse = (islice + 0.5_rt) / ref_ratio_z; - const amrex::Real rel_z = islice_coarse-static_cast(amrex::Math::floor(islice_coarse)); - - auto solution_interp = interpolated_field_xyz{ - {getField(lev-1, WhichSlice::This, component), - getField(lev-1, WhichSlice::Previous1, component), - rel_z}, geom[lev-1]}; + auto solution_interp = interpolated_field_xy{ + getField(lev-1, WhichSlice::This, component), geom[lev-1]}; for (amrex::MFIter mfi(staging_area, DfltMfi); mfi.isValid(); ++mfi) { @@ -728,14 +670,8 @@ Fields::InterpolateFromLev0toLev1 (amrex::Vector const& geom, c HIPACE_PROFILE("Fields::InterpolateFromLev0toLev1()"); constexpr int interp_order = 2; - const amrex::Real ref_ratio_z = Hipace::GetRefRatio(lev)[2]; - const amrex::Real islice_coarse = (islice + 0.5_rt) / ref_ratio_z; - const amrex::Real rel_z = islice_coarse - static_cast(amrex::Math::floor(islice_coarse)); - - auto field_coarse_interp = interpolated_field_xyz{ - {getField(lev-1, WhichSlice::This, component), - getField(lev-1, WhichSlice::Previous1, component), - rel_z}, geom[lev-1]}; + auto field_coarse_interp = interpolated_field_xy{ + getField(lev-1, WhichSlice::This, component), geom[lev-1]}; amrex::MultiFab field_fine = getField(lev, WhichSlice::This, component); for (amrex::MFIter mfi( field_fine, DfltMfi); mfi.isValid(); ++mfi) @@ -772,7 +708,7 @@ Fields::InterpolateFromLev0toLev1 (amrex::Vector const& geom, c void Fields::SolvePoissonExmByAndEypBx (amrex::Vector const& geom, - const MPI_Comm& m_comm_xy, const int lev, const int islice) + const int lev, const int islice) { /* Solves Laplacian(Psi) = 1/epsilon0 * -(rho-Jz/c) and * calculates Ex-c By, Ey + c Bx from grad(-Psi) @@ -802,9 +738,7 @@ Fields::SolvePoissonExmByAndEypBx (amrex::Vector const& geom, if (!m_extended_solve) { /* ---------- Transverse FillBoundary Psi ---------- */ - amrex::ParallelContext::push(m_comm_xy); lhs.FillBoundary(geom[lev].periodicity()); - amrex::ParallelContext::pop(); } InterpolateFromLev0toLev1(geom, lev, "Psi", islice, m_slices_nguards, m_poisson_nguards);