Separate grid for Laser

docs/source/run/parameters.rst

@@ -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`)
@@ -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
 ---------
 
@@ -789,6 +799,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.
@@ -882,11 +896,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`)
+    On 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``.
@@ -919,8 +937,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.

src/Hipace.H

@@ -279,7 +279,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 ();

src/Hipace.cpp

@@ -68,7 +68,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);
@@ -155,7 +155,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 */
@@ -189,15 +192,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()) {
@@ -210,10 +212,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;
@@ -332,7 +331,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);
         }
 
@@ -376,6 +375,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) {
@@ -418,7 +419,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
@@ -481,7 +482,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) {
@@ -518,7 +519,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);
@@ -569,12 +570,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) {
@@ -610,7 +609,7 @@ Hipace::SolveOneSlice (int islice, int step)
 
     m_multi_beam.shiftBeamSlices();
 
-    m_multi_laser.ShiftLaserSlices();
+    m_multi_laser.ShiftLaserSlices(islice);
 }
 
 void
@@ -961,18 +960,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);
         }
     }
 }

src/diagnostics/Diagnostic.H

@@ -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)*/
@@ -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";
 };
 
 
@@ -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; }
@@ -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);
     }
@@ -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: