Skip to content

Commit

Permalink
Remove zeta refinement (#911)
Browse files Browse the repository at this point in the history 
* remove zeta refinement

* add assert
  • Loading branch information
@AlexanderSinn
AlexanderSinn authored Apr 6, 2023
1 parent 7ceaad0 commit 4177aac
Show file tree
Hide file tree
Showing 4 changed files with 51 additions and 134 deletions.
4 changes: 2 additions & 2 deletions src/Hipace.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<amrex::Vector<BeamBins>>& bins);

/** \brief Full evolve on 1 subslice with explicit solver
Expand Down
96 changes: 40 additions & 56 deletions src/Hipace.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -186,6 +186,8 @@ Hipace::Hipace () :
for (int idim=0; idim<AMREX_SPACEDIM; ++idim) patch_lo[idim] = loc_array[idim];
getWithParser(pph, "patch_hi", loc_array);
for (int idim=0; idim<AMREX_SPACEDIM; ++idim) patch_hi[idim] = loc_array[idim];
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(GetRefRatio(1)[2] == 1,
"Mesh refinement in the zeta direction is not supported");
}
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!m_do_MR || !m_multi_beam.AnySpeciesSalame(),
"Cannot use SALAME algorithm with mesh refinement");
Expand Down Expand Up @@ -566,24 +568,24 @@ Hipace::Evolve ()
}

void
Hipace::SolveOneSlice (int islice_coarse, const int ibox, int step,
Hipace::SolveOneSlice (int islice, const int ibox, int step,
const amrex::Vector<amrex::Vector<BeamBins>>& 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) {

if (lev == 1) { // skip all slices which are not existing on level 1
// 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;
}

Expand All @@ -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]);
}
}


Expand Down Expand Up @@ -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);

Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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
Expand All @@ -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 */
Expand Down Expand Up @@ -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 */
Expand All @@ -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
Expand Down Expand Up @@ -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),
Expand Down
3 changes: 1 addition & 2 deletions src/fields/Fields.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<amrex::Geometry> 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
Expand Down
Loading

0 comments on commit 4177aac

Please sign in to comment.