Skip to content

Commit

Permalink
Fix adaptive time step for non-electron beams (#1092)
Browse files Browse the repository at this point in the history
  • Loading branch information
@AlexanderSinn
AlexanderSinn authored Apr 28, 2024
1 parent 482e3b6 commit cb94a31
Show file tree
Hide file tree
Showing 4 changed files with 37 additions and 33 deletions.
4 changes: 2 additions & 2 deletions src/particles/plasma/MultiPlasma.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,12 +63,12 @@ public:
void ExplicitDeposition (Fields& fields, amrex::Vector<amrex::Geometry> const& gm,
int const lev);

/** \brief Return max density, to compute the adaptive time step.
/** \brief Return max charge density, to compute the adaptive time step.
*
* the max is taken across species AND include m_adaptive_density, giving a way to
* specify a density to the adaptive time step calculator even with no plasma species.
*/
amrex::Real maxDensity (amrex::Real z);
amrex::Real maxChargeDensity (amrex::Real z);

/** \brief Gather field values and push particles
*
Expand Down
9 changes: 5 additions & 4 deletions src/particles/plasma/MultiPlasma.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,14 +63,15 @@ MultiPlasma::InitData (amrex::Vector<amrex::BoxArray> slice_ba,
}

amrex::Real
MultiPlasma::maxDensity (amrex::Real z)
MultiPlasma::maxChargeDensity (amrex::Real z)
{
amrex::Real max_density = 0;
amrex::Real max_density = std::abs(m_adaptive_density * get_phys_const().q_e);
for (auto& plasma : m_all_plasmas) {
plasma.UpdateDensityFunction(z);
max_density = amrex::max<amrex::Real>(max_density, plasma.m_density_func(0., 0., z));
max_density = amrex::max<amrex::Real>(
max_density, std::abs(plasma.GetCharge() * plasma.m_density_func(0., 0., z)));
}
return amrex::max(max_density, m_adaptive_density);
return max_density;
}

void
Expand Down
4 changes: 2 additions & 2 deletions src/utils/AdaptiveTimeStep.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,8 +25,8 @@ private:
amrex::Real m_nt_per_betatron = 20.;
/** Upper bound of the time step. Avoid gigantic time step(s) when beam starts near vacuum */
amrex::Real m_dt_max = std::numeric_limits<amrex::Real>::infinity();
/** uz of the slowest particles */
amrex::Real m_min_uz = std::numeric_limits<amrex::Real>::max();
/** uz*mass/charge of the slowest particles */
amrex::Real m_min_uz_mq = std::numeric_limits<amrex::Real>::max();
/** Threshold beam momentum, below which the time step is not decreased */
amrex::Real m_threshold_uz = 2.;
/** Whether to predict the next time steps. More accurate for parallel simulations */
Expand Down
53 changes: 28 additions & 25 deletions src/utils/AdaptiveTimeStep.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,7 @@ AdaptiveTimeStep::BroadcastTimeStep (amrex::Real& dt)
amrex::ParallelDescriptor::Communicator());

// Broadcast minimum uz
MPI_Bcast(&m_min_uz,
MPI_Bcast(&m_min_uz_mq,
1,
amrex::ParallelDescriptor::Mpi_typemap<amrex::Real>::type(),
Hipace::HeadRankID(),
Expand Down Expand Up @@ -171,22 +171,25 @@ AdaptiveTimeStep::CalculateFromMinUz (

const PhysConst phys_const = get_phys_const();
const amrex::Real c = phys_const.c;
const amrex::Real q_e = phys_const.q_e;
const amrex::Real m_e = phys_const.m_e;
const amrex::Real ep0 = phys_const.ep0;

// Extract properties associated with physical size of the box
const int nbeams = beams.get_nbeams();
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(nbeams >= 1,
"Must have at least one beam to use adaptive time step");
const int numprocs = Hipace::m_numprocs;

amrex::Vector<amrex::Real> new_dts;
new_dts.resize(nbeams);
amrex::Vector<amrex::Real> beams_min_uz;
beams_min_uz.resize(nbeams, std::numeric_limits<amrex::Real>::max());
amrex::Vector<amrex::Real> beams_min_uz_mq;
beams_min_uz_mq.resize(nbeams, std::numeric_limits<amrex::Real>::max());

for (int ibeam = 0; ibeam < nbeams; ibeam++) {
new_dts[ibeam] = dt;

const auto& beam = beams.getBeam(ibeam);
if (beam.m_charge == 0.) { continue; }
const amrex::Real mass_charge_ratio = beam.m_mass / beam.m_charge;

AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
m_timestep_data[ibeam][WhichDouble::SumWeights] != 0,
Expand All @@ -203,19 +206,17 @@ AdaptiveTimeStep::CalculateFromMinUz (
std::min(std::max(sigma_uz_dev, m_timestep_data[ibeam][WhichDouble::MinUz]),
max_supported_uz);

beams_min_uz[ibeam] = chosen_min_uz*beam.m_mass*beam.m_mass/m_e/m_e;

if (Hipace::m_verbose >=2 ){
amrex::Print()<<"Minimum gamma of beam " << ibeam <<
" to calculate new time step: " << beams_min_uz[ibeam] << "\n";
" to calculate new time step: " << chosen_min_uz << "\n";
}

if (beams_min_uz[ibeam] < m_threshold_uz) {
if (chosen_min_uz < m_threshold_uz) {
amrex::Print()<<"WARNING: beam particles of beam "<< ibeam <<
" have non-relativistic velocities!\n";
}
beams_min_uz[ibeam] = std::max(beams_min_uz[ibeam], m_threshold_uz);
new_dts[ibeam] = dt;
chosen_min_uz = std::max(chosen_min_uz, m_threshold_uz);
beams_min_uz_mq[ibeam] = std::abs(chosen_min_uz * mass_charge_ratio);

/*
Calculate the time step for this beam used in the next time iteration
Expand All @@ -228,25 +229,29 @@ AdaptiveTimeStep::CalculateFromMinUz (
*/
amrex::Real new_dt = dt;
amrex::Real new_time = t;
amrex::Real min_uz = beams_min_uz[ibeam];
amrex::Real min_uz = chosen_min_uz;
const int niter = m_adaptive_predict_step ? numprocs : 1;
for (int i = 0; i < niter; i++)
{
amrex::Real plasma_density = plasmas.maxDensity(c * new_time);
AMREX_ALWAYS_ASSERT_WITH_MESSAGE( plasma_density > 0.,
amrex::Real plasma_charge_density = plasmas.maxChargeDensity(c * new_time);
AMREX_ALWAYS_ASSERT_WITH_MESSAGE( plasma_charge_density > 0.,
"A >0 plasma density must be specified to use an adaptive time step.");
min_uz += m_timestep_data[ibeam][WhichDouble::MinAcc] * new_dt;
if (m_adaptive_gather_ez) {
min_uz += m_timestep_data[ibeam][WhichDouble::MinAcc] * new_dt;
}
// Just make sure min_uz is >0, to avoid nans below.
min_uz = std::max(min_uz, 0.001_rt*m_threshold_uz);
const amrex::Real omega_p = std::sqrt(plasma_density * q_e * q_e / (ep0 * m_e));
amrex::Real omega_b = omega_p / std::sqrt(2. * min_uz);
amrex::Real omega_b = std::sqrt(plasma_charge_density /
(2. * std::abs(min_uz * mass_charge_ratio) * ep0));
new_dt = 2. * MathConst::pi / omega_b / m_nt_per_betatron;
new_time += new_dt;
if (min_uz > m_threshold_uz) new_dts[ibeam] = new_dt;
if (min_uz > m_threshold_uz) {
new_dts[ibeam] = new_dt;
}
}
}
// Store min uz across beams, used in the phase advance method
m_min_uz = *std::min_element(beams_min_uz.begin(), beams_min_uz.end());
m_min_uz_mq = *std::min_element(beams_min_uz_mq.begin(), beams_min_uz_mq.end());
/* set the new time step */
dt = *std::min_element(new_dts.begin(), new_dts.end());
// Make sure the new time step is smaller than the upper bound
Expand Down Expand Up @@ -341,17 +346,15 @@ AdaptiveTimeStep::CalculateFromDensity (amrex::Real t, amrex::Real& dt, MultiPla
amrex::Real phase_advance0 = 0.;

// Get plasma density at beginning of step
const amrex::Real plasma_density0 = plasmas.maxDensity(pc.c * t);
const amrex::Real omgp0 = std::sqrt(plasma_density0 * pc.q_e * pc.q_e / (pc.ep0 * pc.m_e));
amrex::Real omgb0 = omgp0 / std::sqrt(2. *m_min_uz);
const amrex::Real plasma_charge_density0 = plasmas.maxChargeDensity(pc.c * t);
const amrex::Real omgb0 = std::sqrt(plasma_charge_density0 / (2. * m_min_uz_mq * pc.ep0));

// Numerically integrate the phase advance from t to t+dt. The time step is reduced such that
// the expected phase advance equals that of a uniform plasma up to a tolerance level.
for (int i = 0; i < m_adaptive_phase_substeps; i++)
{
const amrex::Real plasma_density = plasmas.maxDensity(pc.c * (t+i*dt_sub));
const amrex::Real omgp = std::sqrt(plasma_density * pc.q_e * pc.q_e / (pc.ep0 * pc.m_e));
amrex::Real omgb = omgp / std::sqrt(2. *m_min_uz);
const amrex::Real plasma_charge_density = plasmas.maxChargeDensity(pc.c * (t+i*dt_sub));
const amrex::Real omgb = std::sqrt(plasma_charge_density / (2. * m_min_uz_mq * pc.ep0));
phase_advance += omgb * dt_sub;
phase_advance0 += omgb0 * dt_sub;
if(std::abs(phase_advance - phase_advance0) >
Expand Down

0 comments on commit cb94a31

Please sign in to comment.