Skip to content

Commit

Permalink
Separate grid for Laser (#1106)
Browse files Browse the repository at this point in the history
  • Loading branch information
@AlexanderSinn
AlexanderSinn authored May 29, 2024
1 parent d8efbed commit cb4858b
Show file tree
Hide file tree
Showing 12 changed files with 649 additions and 302 deletions.
36 changes: 29 additions & 7 deletions docs/source/run/parameters.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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`)
Expand All @@ -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
---------

Expand Down Expand Up @@ -805,6 +815,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.
Expand Down Expand Up @@ -898,11 +912,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``.

* ``<diag name> or diagnostic.level`` (`integer`) optional (default `0`)
From which mesh refinement level the diagnostics should be collected.
If ``<diag name>`` is equal to ``lev1``, the default for this parameter becomes 1 etc.
* ``<diag name> or diagnostic.base_geometry`` (`string`) optional (default `level_0`)
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 ``<diag name>`` is equal to ``lev0 lev1 lev2 laser_diag``, the default for this parameter
becomes ``level_0 level_1 level_2 laser``respectively.
* ``<diag name>.output_period`` (`integer`) optional (default `0`)
Output period for fields. No output is given for ``<diag name>.output_period = 0``.
Expand Down Expand Up @@ -935,8 +953,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_<plasma name>`` 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_<field name>`` after it has been added. If a field is added and removed
multiple times, the last occurrence takes precedence.

* ``<diag name> or diagnostic.patch_lo`` (3 `float`) optional (default `-infinity -infinity -infinity`)
Lower limit for the diagnostic grid.
Expand Down
2 changes: 1 addition & 1 deletion src/Hipace.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 ();
Expand Down
51 changes: 24 additions & 27 deletions src/Hipace.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,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);
Expand Down Expand Up @@ -162,7 +162,10 @@ Hipace::Hipace () :

MakeGeometry();

m_use_laser = m_multi_laser.m_use_laser;
// use level 0 as default for laser geometry
m_multi_laser.MakeLaserGeometry(m_3D_geom[0]);

m_use_laser = m_multi_laser.UseLaser();

queryWithParser(pph, "collisions", m_collision_names);
/** Initialize the collision objects */
Expand Down Expand Up @@ -196,15 +199,14 @@ Hipace::InitData ()
amrex::Print() << "using the leapfrog plasma particle pusher\n";
#endif

m_multi_laser.InitData();

for (int lev=0; lev<m_N_level; ++lev) {
m_fields.AllocData(lev, m_3D_geom[lev], m_slice_ba[lev], m_slice_dm[lev]);
if (lev==0) {
// laser inits only on level 0
m_multi_laser.InitData(m_slice_ba[0], m_slice_dm[0], m_3D_geom[0]);
}
m_diags.Initialize(lev, m_multi_laser.m_use_laser);
}

m_diags.Initialize(m_N_level, m_multi_laser.UseLaser());

m_initial_time = m_multi_beam.InitData(m_3D_geom[0]);

if (Hipace::HeadRank()) {
Expand All @@ -217,10 +219,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;
Expand Down Expand Up @@ -339,7 +338,7 @@ Hipace::Evolve ()

const amrex::Box& bx = m_3D_ba[0][0];

if (m_multi_laser.m_use_laser) {
if (m_multi_laser.UseLaser()) {
AMREX_ALWAYS_ASSERT(!m_adaptive_time_step.m_do_adaptive_time_step);
}

Expand Down Expand Up @@ -383,6 +382,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) {
Expand Down Expand Up @@ -425,7 +426,7 @@ Hipace::Evolve ()
m_fields.InSituWriteToFile(step, m_physical_time, m_3D_geom[0], m_max_step, m_max_time);
m_multi_beam.InSituWriteToFile(step, m_physical_time, m_3D_geom[0], m_max_step, m_max_time);
m_multi_plasma.InSituWriteToFile(step, m_physical_time, m_3D_geom[0], m_max_step, m_max_time);
m_multi_laser.InSituWriteToFile(step, m_physical_time, m_3D_geom[0], m_max_step, m_max_time);
m_multi_laser.InSituWriteToFile(step, m_physical_time, m_max_step, m_max_time);

if (!m_explicit) {
// averaging predictor corrector loop diagnostics
Expand Down Expand Up @@ -488,7 +489,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<current_N_level; ++lev) {
Expand Down Expand Up @@ -525,7 +526,7 @@ Hipace::SolveOneSlice (int islice, int step)

// Advance laser slice by 1 step using chi
// no MR for laser
m_multi_laser.AdvanceSlice(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);
Expand Down Expand Up @@ -576,12 +577,10 @@ Hipace::SolveOneSlice (int islice, int step)
m_fields.InSituComputeDiags(step, m_physical_time, islice, m_3D_geom[0], m_max_step, m_max_time);

// get laser insitu diagnostics
m_multi_laser.InSituComputeDiags(step, m_physical_time, islice, m_3D_geom[0], m_max_step, m_max_time);
m_multi_laser.InSituComputeDiags(step, m_physical_time, islice, m_max_step, m_max_time);

// copy fields (and laser) to diagnostic array
for (int lev=0; lev<current_N_level; ++lev) {
FillFieldDiagnostics(lev, islice);
}
FillFieldDiagnostics(current_N_level, islice);

// plasma ionization
for (int lev=0; lev<current_N_level; ++lev) {
Expand Down Expand Up @@ -617,7 +616,7 @@ Hipace::SolveOneSlice (int islice, int step)

m_multi_beam.shiftBeamSlices();

m_multi_laser.ShiftLaserSlices();
m_multi_laser.ShiftLaserSlices(islice);
}

void
Expand Down Expand Up @@ -968,18 +967,16 @@ Hipace::InitDiagnostics (const int step)
m_openpmd_writer.InitBeamData(m_multi_beam, getDiagBeamNames());
}
#endif
for (int lev=0; lev<m_N_level; ++lev) {
m_diags.ResizeFDiagFAB(m_3D_geom[lev].Domain(), lev, m_3D_geom[lev],
step, m_max_step, m_physical_time, m_max_time);
}
m_diags.ResizeFDiagFAB(m_3D_geom, m_multi_laser.GetLaserGeom(),
step, m_max_step, m_physical_time, m_max_time);
}

void
Hipace::FillFieldDiagnostics (const int lev, int islice)
Hipace::FillFieldDiagnostics (const int current_N_level, int islice)
{
for (auto& fd : m_diags.getFieldData()) {
if (fd.m_level == lev && fd.m_has_field) {
m_fields.Copy(lev, islice, fd, m_3D_geom[lev], m_multi_laser);
if (fd.m_has_field) {
m_fields.Copy(current_N_level, islice, fd, m_3D_geom, m_multi_laser);
}
}
}
Expand Down
27 changes: 13 additions & 14 deletions src/diagnostics/Diagnostic.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,11 @@
struct FieldDiagnosticData
{
std::string m_diag_name; /**< name used for input parameters and in the output */
int m_level = 0; /**< MR level */
enum struct geom_type {
field,
laser
} m_base_geom_type = geom_type::field; /**< which geometry the diagnostics is based on */
int m_level = 0; /**< MR level when field_geom is used */
int m_slice_dir; /**< Slicing direction */
bool m_include_ghost_cells = false; /**< if ghost cells are included in output */
bool m_use_custom_size_lo = false; /**< if a user defined diagnostics size should be used (lo)*/
Expand All @@ -30,23 +34,19 @@ struct FieldDiagnosticData
/** 3D array with upper ends of the diagnostics grid */
amrex::RealVect m_diag_hi {0., 0., 0.};
amrex::IntVect m_diag_coarsen; /**< xyz coarsening ratio (positive) */
bool m_do_laser {false}; /**< Whether to output the laser */
int m_nfields; /**< Number of physical fields to write */
amrex::Vector<std::string> m_comps_output; /**< Component names to Write to output file */
/** Component indexes to Write to output file */
amrex::Gpu::DeviceVector<int> m_comps_output_idx;
/** Vector over levels, all fields */
amrex::FArrayBox m_F;
using complex_type = amrex::GpuComplex<amrex::Real>;
/** FAB for laser */
amrex::BaseFab<complex_type> m_F_laser;
amrex::BaseFab<amrex::GpuComplex<amrex::Real>> 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";
};


Expand All @@ -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<std::string>& getBeamNames () { return m_output_beam_names; }
Expand All @@ -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);
}
Expand Down Expand Up @@ -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<amrex::Geometry>& field_geom,
amrex::Geometry const& laser_geom, int output_step, int max_step,
amrex::Real output_time, amrex::Real max_time);

private:
Expand Down
Loading

0 comments on commit cb4858b

Please sign in to comment.