[WIP] Set time step with parser

src/Hipace.H

@@ -257,7 +257,13 @@ public:
      * or radiation reaction in normalized units */
     inline static amrex::Real m_background_density_SI = 0.;
     /** Time step for the beam evolution */
-    amrex::Real m_dt = +0.0;
+    amrex::Parser m_parser_dt;/**< owns data for dt_str */
+    /** stores the output function of makeFunctionWithParser for profile_real_str */
+    amrex::ParserExecutor<1> m_exe_dt;
+    /** stores time step of current time steps */
+    amrex::Real m_dt = 0;
+    /** stores time step of previous time steps */
+    amrex::Real m_dt_pre = 0;
     /** Number of binary collisions */
     inline static int m_ncollisions = 0;
     /** Adaptive time step instance */

src/Hipace.cpp

@@ -97,12 +97,6 @@ Hipace::Hipace () :
 
     amrex::ParmParse pph("hipace");
 
-    std::string str_dt {""};
-    queryWithParser(pph, "dt", str_dt);
-    if (str_dt != "adaptive") {
-        queryWithParser(pph, "dt", m_dt);
-        m_max_time = std::copysign(m_max_time, m_dt);
-    }
     queryWithParser(pph, "max_time", m_max_time);
     queryWithParser(pph, "verbose", m_verbose);
     m_numprocs = amrex::ParallelDescriptor::NProcs();
@@ -407,30 +401,69 @@ Hipace::Evolve ()
     HIPACE_PROFILE("Hipace::Evolve()");
     const double start_time = amrex::second();
     const int rank = amrex::ParallelDescriptor::MyProc();
-
+    m_comm = amrex::ParallelDescriptor::Communicator();
+    amrex::ParmParse pph("hipace");
     // now each rank starts with its own time step and writes to its own file. The first rank starts with step 0
     for (int step = rank; step <= m_max_step; step += m_numprocs)
     {
         ResetAllQuantities();
+        m_physical_time = step == 0 ? m_initial_time : m_multi_buffer.get_time();
 
         const amrex::Box& bx = m_3D_ba[0][0];
 
         if (m_multi_laser.UseLaser()) {
             AMREX_ALWAYS_ASSERT(!m_adaptive_time_step.m_do_adaptive_time_step);
         }
 
-        m_physical_time = step == 0 ? m_initial_time : m_multi_buffer.get_time();
-
         if (m_physical_time == std::numeric_limits<amrex::Real>::infinity()) {
             if (step+1 <= m_max_step && !m_has_last_step) {
                 m_multi_buffer.put_time(m_physical_time);
             }
             break;
         }
+        std::string str_dt {"0."};
+        queryWithParser(pph, "dt", str_dt);
+        constexpr int tag_time_step = 100;  // Use an integer for the tag
+
+        if (str_dt != "adaptive") {
+            m_exe_dt = makeFunctionWithParser<1>(str_dt, m_parser_dt, {"t"});
+
+            if (step == 0) {
+                m_dt = m_exe_dt(m_physical_time);
+            } else {
+                MPI_Request recv_request;
+                MPI_Irecv(
+                    &m_dt_pre,
+                    1,
+                    amrex::Real,
+                    rank - 1,
+                    tag_time_step,
+                    amrex::ParallelDescriptor::m_comm,
+                    &recv_request);
+                // You should later call MPI_Wait or MPI_Test to complete it
+            }
 
-        m_adaptive_time_step.CalculateFromDensity(m_physical_time, m_dt, m_multi_plasma);
+            m_dt = m_exe_dt(m_physical_time);
+
+            MPI_Request send_request;
+            MPI_Isend(
+                &m_dt,
+                1,
+                amrex::Real,
+                rank + 1,
+                tag_time_step,
+                amrex::ParallelDescriptor::m_comm,
+                &send_request);
+            // You should also wait or test this one
+            MPI_Wait(&recv_request, MPI_STATUS_IGNORE);
+            MPI_Wait(&send_request, MPI_STATUS_IGNORE);
+            m_max_time = std::copysign(m_max_time, m_dt);
+        }
+        else{
+            m_adaptive_time_step.CalculateFromDensity(m_physical_time, m_dt, m_multi_plasma);
+        }
 
-        amrex::Real next_time = 0.;
+        amrex::Real next_time = 0. ;
 
         // adjust time step to reach max_time
         if (m_physical_time == m_max_time) {
@@ -645,7 +678,7 @@ Hipace::SolveOneSlice (int islice, int step)
 
     // Advance laser slice by 1 step using chi
     // no MR for laser
-    m_multi_laser.AdvanceSlice(islice, m_fields, m_dt, step, m_3D_geom[0]);
+    m_multi_laser.AdvanceSlice(islice, m_fields, m_dt, m_dt_pre, 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);

src/laser/MultiLaser.H

@@ -116,28 +116,31 @@ public:
     /** Wrapper function to advance a laser slice by 1 time step.
      * \param[in] islice slice index
      * \param[in] fields Field object
-     * \param[in] dt time step of the simulation
+     * \param[in] dt time step of the simulation of the current rank
      * \param[in] step current iteration. Needed because step 0 needs a specific treatment.
      * \param[in] geom_field_lev0 Geometry of the fields on MR level 0
+     * \param[in] dt_pre  time step of the previous rank
      */
-    void AdvanceSlice (const int islice, const Fields& fields, amrex::Real dt, int step,
+    void AdvanceSlice (const int islice, const Fields& fields, const amrex::Real dt, const amrex::Real dt_pre, int step,
                        amrex::Geometry const& geom_field_lev0);
 
     /** Advance a laser slice by 1 time step using a multigrid solver.
      * The complex phase of the envelope is evaluated on-axis only, but can be generalized to everywhere.
      *
-     * \param[in] dt time step of the simulation
+     * \param[in] dt time step of the simulation of the current rank
      * \param[in] step current iteration. Needed because step 0 needs a specific treatment.
+     * \param[in] dt_pre  time step of the previous rank
      */
-    void AdvanceSliceMG (amrex::Real dt, int step);
+    void AdvanceSliceMG (const amrex::Real dt, const amrex::Real dt_pre, int step);
 
     /** Advance a laser slice by 1 time step using a FFT solver.
      * The complex phase of the envelope is evaluated on-axis only.
      *
-     * \param[in] dt time step of the simulation
+     * \param[in] dt time step of the simulation of the current rank
      * \param[in] step current iteration. Needed because step 0 needs a specific treatment.
+     * \param[in] dt_pre  time step of the previous rank
      */
-    void AdvanceSliceFFT (amrex::Real dt, int step);
+    void AdvanceSliceFFT (const amrex::Real dt, const amrex::Real dt_pre, int step);
 
     /** Initialize 1 longitudinal slice of the laser, and store it in n00j00 (current time step)
      * and nm1j00 (previous time step).

src/laser/MultiLaser.cpp

@@ -415,7 +415,7 @@ MultiLaser::InterpolateChi (const Fields& fields, amrex::Geometry const& geom_fi
 }
 
 void
-MultiLaser::AdvanceSlice (const int islice, const Fields& fields, amrex::Real dt, int step,
+MultiLaser::AdvanceSlice (const int islice, const Fields& fields,const amrex::Real dt, const amrex::Real dt_pre, int step,
                           amrex::Geometry const& geom_field_lev0)
 {
 
@@ -426,16 +426,16 @@ MultiLaser::AdvanceSlice (const int islice, const Fields& fields, amrex::Real dt
     InterpolateChi(fields, geom_field_lev0);
 
     if (m_solver_type == "multigrid") {
-        AdvanceSliceMG(dt, step);
+        AdvanceSliceMG(dt, dt_pre, step);
     } else if (m_solver_type == "fft") {
-        AdvanceSliceFFT(dt, step);
+        AdvanceSliceFFT(dt, dt_pre, step);
     } else {
         amrex::Abort("laser.solver_type must be fft or multigrid");
     }
 }
 
 void
-MultiLaser::AdvanceSliceMG (amrex::Real dt, int step)
+MultiLaser::AdvanceSliceMG (const amrex::Real dt, const amrex::Real dt_pre, int step)
 {
 
     HIPACE_PROFILE("MultiLaser::AdvanceSliceMG()");
@@ -543,7 +543,7 @@ MultiLaser::AdvanceSliceMG (amrex::Real dt, int step)
                 using namespace WhichLaserSlice;
                 // Transverse Laplacian of real and imaginary parts of A_j^n-1
                 amrex::Real lapR, lapI;
-                if (step == 0) {
+                if (step == 0 || dt_pre!= dt) {
                     lapR = i>imin && i<imax && j>jmin && j<jmax ?
                         (arr(i+1, j, n00j00_r)+arr(i-1, j, n00j00_r)-2._rt*arr(i, j, n00j00_r))/(dx*dx) +
                         (arr(i, j+1, n00j00_r)+arr(i, j-1, n00j00_r)-2._rt*arr(i, j, n00j00_r))/(dy*dy) : 0._rt;
@@ -566,7 +566,7 @@ MultiLaser::AdvanceSliceMG (amrex::Real dt, int step)
                     acoeff_real_scalar + arr(i, j, chi) : acoeff_real_scalar;
 
                 Complex rhs;
-                if (step == 0) {
+                if (step == 0 || dt_pre!= dt) {
                     // First time step: non-centered push to go
                     // from step 0 to step 1 without knowing -1.
                     const Complex an00jp1 = arr(i, j, n00jp1_r) + I * arr(i, j, n00jp1_i);
@@ -615,7 +615,7 @@ MultiLaser::AdvanceSliceMG (amrex::Real dt, int step)
 }
 
 void
-MultiLaser::AdvanceSliceFFT (const amrex::Real dt, int step)
+MultiLaser::AdvanceSliceFFT (const amrex::Real dt, const amrex::Real dt_pre, int step)
 {
 
     HIPACE_PROFILE("MultiLaser::AdvanceSliceFFT()");