Skip to content

Commit

Permalink
Better defaults for production (in particular adaptive time step) (#893)
Browse files Browse the repository at this point in the history 
* better defaults

* slight value change

* lack of newline messes up with post-processing routine

* export threshold uz to user for adaptive time stepg

* reduce nt_per_betatron default

* fix amd compilation issue

* revert to 1 subcycle

* update benchmark

* fix doc

---------

Co-authored-by: Tools <tools@desy.de>
  • Loading branch information
@MaxThevenet
MaxThevenet and Tools authored Mar 14, 2023
1 parent f788b3a commit 9c73bff
Show file tree
Hide file tree
Showing 4 changed files with 66 additions and 53 deletions.
16 changes: 10 additions & 6 deletions docs/source/run/parameters.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,35 +69,39 @@ General parameters
Only used if ``hipace.dt = adaptive``. Upper bound of the adaptive time step: if the computed adaptive time step is is larger than ``dt_max``, then ``dt_max`` is used instead.
Useful when the plasma profile starts with a very low density (e.g. in the presence of a realistic density ramp), to avoid unreasonably large time steps.

* ``hipace.nt_per_betatron`` (`Real`) optional (default `40.`)
* ``hipace.nt_per_betatron`` (`Real`) optional (default `20.`)
Only used when using adaptive time step (see ``hipace.dt`` above).
Number of time steps per betatron period (of the full blowout regime).
The time step is given by :math:`\omega_{\beta}\Delta t = 2 \pi/N`
(:math:`N` is ``nt_per_betatron``) where :math:`\omega_{\beta}=\omega_p/\sqrt{2\gamma}` with
:math:`\omega_p` the plasma angular frequency and :math:`\gamma` is an average of Lorentz
factors of the slowest particles in all beams.

* ``hipace.adaptive_predict_step`` (`bool`) optional (default `0`)
* ``hipace.adaptive_predict_step`` (`bool`) optional (default `1`)
Only used when using adaptive time step (see ``hipace.dt`` above).
If true, the current Lorentz factor and accelerating field on the beams are used to predict the (adaptive) ``dt`` of the next time steps.
This prediction is used to better estimate the betatron frequency at the beginning of the next step performed by the current rank.
It improves accuracy for parallel simulations (with significant deceleration and/or z-dependent plasma profile).
Note: should be on by default once good defaults are determined.

* ``hipace.adaptive_control_phase_advance`` (`bool`) optional (default `0`)
* ``hipace.adaptive_control_phase_advance`` (`bool`) optional (default `1`)
Only used when using adaptive time step (see ``hipace.dt`` above).
If true, a test on the phase advance sets the time step so it matches the phase advance expected for a uniform plasma (to a certain tolerance).
This should improve the accuracy in the presence of density gradients.
Note: should be on by default once good defaults are determined.

* ``hipace.adaptive_phase_tolerance`` (`Real`) optional (default `5.e-4`)
* ``hipace.adaptive_phase_tolerance`` (`Real`) optional (default `4.e-4`)
Only used when using adaptive time step (see ``hipace.dt`` above) and ``adaptive_control_phase_advance``.
Tolerance for the controlled phase advance described above (lower is more accurate, but should result in more time steps).

* ``hipace.adaptive_phase_substeps`` (`int`) optional (default `100`)
* ``hipace.adaptive_phase_substeps`` (`int`) optional (default `2000`)
Only used when using adaptive time step (see ``hipace.dt`` above) and ``adaptive_control_phase_advance``.
Number of sub-steps in the controlled phase advance described above (higher is more accurate, but should be slower).

* ``hipace.adaptive_threshold_uz`` (`Real`) optional (default `2.`)
Only used when using adaptive time step (see ``hipace.dt`` above).
Threshold beam momentum, below which the time step is not decreased (to avoid arbitrarily small time steps).

* ``hipace.normalized_units`` (`bool`) optional (default `0`)
Using normalized units in the simulation.

Expand Down Expand Up @@ -656,4 +660,4 @@ Diagnostic parameters
Lower limit for the diagnostic grid.

* ``diagnostic.patch_hi`` (3 `float`) optional (default `infinity infinity infinity`)
Upper limit for the diagnostic grid.
Upper limit for the diagnostic grid.
12 changes: 7 additions & 5 deletions src/utils/AdaptiveTimeStep.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,22 +22,24 @@ private:
amrex::Vector<amrex::Vector<amrex::Real>> m_timestep_data;

/** Number of time steps per betatron period for the adaptive time step */
amrex::Real m_nt_per_betatron = 40.;
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();
/** 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 */
bool m_adaptive_predict_step = false;
bool m_adaptive_predict_step = true;
/** If true, a test on the phase advance sets the time step so it matches the phase advance
* expected for a uniform plasma. Relevant in the presence of density gradients.
* The tolerance on the phase advance difference is controlled by m_adaptive_phase_tolerance */
bool m_adaptive_control_phase_advance = false;
bool m_adaptive_control_phase_advance = true;
/** Phase shift tolerance. Relevant when density gradients are present.
* Lower is more accurate. */
amrex::Real m_adaptive_phase_tolerance = 5.e-4;
amrex::Real m_adaptive_phase_tolerance = 4.e-4;
/** Number of substeps on which the phase advance is monitored. */
int m_adaptive_phase_substeps = 100;
int m_adaptive_phase_substeps = 2000;

public:
/** Whether to use an adaptive time step */
Expand Down
23 changes: 15 additions & 8 deletions src/utils/AdaptiveTimeStep.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ AdaptiveTimeStep::AdaptiveTimeStep (const int nbeams)
m_do_adaptive_time_step = true;
queryWithParser(ppa, "nt_per_betatron", m_nt_per_betatron);
queryWithParser(ppa, "dt_max", m_dt_max);
queryWithParser(ppa, "adaptive_threshold_uz", m_threshold_uz);
queryWithParser(ppa, "adaptive_phase_tolerance", m_adaptive_phase_tolerance);
queryWithParser(ppa, "adaptive_predict_step", m_adaptive_predict_step);
queryWithParser(ppa, "adaptive_control_phase_advance", m_adaptive_control_phase_advance);
Expand All @@ -49,8 +50,11 @@ AdaptiveTimeStep::BroadcastTimeStep (amrex::Real& dt, MPI_Comm a_comm_z, int a_n
{
if (m_do_adaptive_time_step == 0) return;

// Broadcast time step
MPI_Bcast(&dt, 1, amrex::ParallelDescriptor::Mpi_typemap<amrex::Real>::type(),
a_numprocs_z - 1, a_comm_z);

// Broadcast minimum uz
MPI_Bcast(&m_min_uz, 1, amrex::ParallelDescriptor::Mpi_typemap<amrex::Real>::type(),
a_numprocs_z - 1, a_comm_z);
}
Expand Down Expand Up @@ -168,7 +172,7 @@ AdaptiveTimeStep::Calculate (
const auto& beam = beams.getBeam(ibeam);

AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_timestep_data[ibeam][WhichDouble::SumWeights] != 0,
"The sum of all weights is 0! Probably no beam particles are initialized");
"The sum of all weights is 0! Probably no beam particles are initialized\n");
const amrex::Real mean_uz = m_timestep_data[ibeam][WhichDouble::SumWeightsTimesUz]
/m_timestep_data[ibeam][WhichDouble::SumWeights];
const amrex::Real sigma_uz = std::sqrt(std::abs(m_timestep_data[ibeam][WhichDouble::SumWeightsTimesUzSquared]
Expand All @@ -181,17 +185,17 @@ AdaptiveTimeStep::Calculate (
max_supported_uz);
m_min_uz = std::min(m_min_uz,
chosen_min_uz*beam.m_mass*beam.m_mass/m_e/m_e);
m_min_uz = std::max(m_min_uz, 1._rt);

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

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

Expand All @@ -209,11 +213,13 @@ AdaptiveTimeStep::Calculate (
{
plasma_density = plasmas.maxDensity(c * new_time);
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);
new_dt = 2. * MathConst::pi / omega_b / m_nt_per_betatron;
new_time += new_dt;
if (min_uz > 1) new_dts[ibeam] = new_dt;
if (min_uz > m_threshold_uz) new_dts[ibeam] = new_dt;
}
}
/* set the new time step */
Expand All @@ -226,7 +232,7 @@ AdaptiveTimeStep::Calculate (
void
AdaptiveTimeStep::CalculateFromDensity (amrex::Real t, amrex::Real& dt, MultiPlasma& plasmas)
{
HIPACE_PROFILE("AdaptiveTimeStep::Calculate()");
HIPACE_PROFILE("AdaptiveTimeStep::CalculateFromDensity()");

using namespace amrex::literals;

Expand Down Expand Up @@ -254,8 +260,9 @@ AdaptiveTimeStep::CalculateFromDensity (amrex::Real t, amrex::Real& dt, MultiPla
phase_advance += omgb * dt_sub;
phase_advance0 += omgb0 * dt_sub;
if(std::abs(phase_advance - phase_advance0) >
2.*MathConst::pi*m_adaptive_phase_tolerance)
2.*MathConst::pi*m_adaptive_phase_tolerance/m_nt_per_betatron)
{
if (i==0) amrex::AllPrint()<<"WARNING: adaptive time step exits at first substep. Consider increasing hipace.adaptive_phase_substeps!\n";
dt = i*dt_sub;
return;
}
Expand Down
68 changes: 34 additions & 34 deletions tests/checksum/benchmarks_json/production.SI.2Rank_pwfa.json
View file
Original file line number Diff line number Diff line change
@@ -1,46 +1,46 @@
{
"lev=0": {
"Bx": 0.02197551604223,
"By": 8966.1351421771,
"Bz": 0.031374152137171,
"ExmBy": 3100334618693.1,
"EypBx": 8141775.7529099,
"Ez": 2397186218239.8,
"Psi": 119994784.5044,
"Sx": 77320436549663.0,
"Sy": 128903797.14308,
"chi": 4382910301271.3,
"jx": 159741118147660.0,
"jx_beam": 159899383408.72,
"jy": 1022627377.0168,
"jy_beam": 111931.8410149,
"jz": 135745902812450.0,
"jz_beam": 484032338457010.0,
"rho": 2254468.0648512,
"rho_beam": 0.0
},
"driver": {
"charge": 1.602176634e-13,
"id": 499935519084,
"mass": 9.1093837015e-25,
"ux": 1037910.4608001,
"uy": 1036315.4362785,
"uz": 983566459.0206799,
"w": 3744396137.536,
"x": 3.12584205475840,
"y": 3.11138512226390018,
"z": 23.93368200407
},
"lev=0": {
"Bx": 0.021977466274843,
"By": 8966.2786868988,
"Bz": 0.031374687381087,
"ExmBy": 3100310233053.2,
"EypBx": 8141258.83133830,
"Ez": 2397160352859.7998,
"Psi": 119993348.4861,
"Sx": 77328010812514.0,
"Sy": 128946289.75447,
"chi": 4382917081936.4,
"jx": 159740165479810.0,
"jx_beam": 160257422032.839996,
"jy": 1022704886.79560005664,
"jy_beam": 106573.87189609,
"jz": 135746558371620.0,
"jz_beam": 484051765931080.0,
"rho": 2254515.36896609980,
"rho_beam": 0.0
"x": 3.1261792894625,
"y": 3.1117169618229,
"z": 23.933678838624,
"ux": 1037252.131241,
"uy": 1035672.4179983,
"uz": 983582204.84356,
"w": 3744396137.536
},
"witness": {
"charge": 1.602176634e-13,
"id": 750000500000,
"mass": 9.1093837015e-25,
"ux": 1131666.1776384999975,
"uy": 1136188.4040961,
"uz": 1032049455.7666,
"w": 1248301814.8922,
"x": 2.4396528434158,
"y": 2.4453187801483,
"z": 160.0073173040800043
"x": 2.4414326107363,
"y": 2.4470915062159,
"z": 160.00730897915,
"ux": 1130244.8645185,
"uy": 1134758.5097234,
"uz": 1032012856.823,
"w": 1248301814.8922
}
}

0 comments on commit 9c73bff

Please sign in to comment.