Update BackTransformed diagnostics to take into account arbitrary moving window velocity (BLAST-WarpX#5341)

In the `development` branch, the `BackTransformed` diagnostics assume
that the moving window moves exactly at the speed of light. This PR
generalizes the code for arbitrary moving window velocity.

This PR does not add an automated test, but the upcoming PR BLAST-WarpX#5337 will
add a test which features a moving window with a speed different than
`c`.

This is a follow-up of BLAST-WarpX#5226, which modified the transformation of the
simulation box coordinates for arbitrary moving window velocity, but did
not yet update the `BackTransformed` diagnostic code.

Author: RemiLehe

Docs/source/usage/faq.rst

@@ -74,10 +74,10 @@ Several BTD quantities differ slightly from the lab frame domain described in th
 In the following discussion, we will use a subscript input (e.g. :math:`\Delta z_{\rm input}`) to denote properties of the lab frame domain.
 
 
-- The first back-transformed diagnostic (BTD) snapshot may not occur at :math:`t=0`. Rather, it occurs at :math:`t_0=\frac{z_{max}}c \beta(1+\beta)\gamma^2`. This is the first time when the boosted frame can complete the snapshot.
+- The first back-transformed diagnostic (BTD) snapshot may not occur at :math:`t=0`. Rather, it occurs at :math:`t_0=\frac{z_{max}}c \beta/(1 - \beta \beta_{mw})`, where :math:`\beta_{mw}` represents the speed of the moving window. This is the first time when the boosted frame can complete the snapshot.
 - The grid spacing of the BTD snapshot is different from the grid spacing indicated in the input script. It is given by :math:`\Delta z_{\rm grid,snapshot}=\frac{c\Delta t_{\rm boost}}{\gamma\beta}`.  For a CFL-limited time step, :math:`\Delta z_{\rm grid,snapshot}\approx \frac{1+\beta}{\beta} \Delta z_{\rm input}\approx 2 \Delta z_{\rm input}`. Hence in many common use cases at large boost, it is expected that the BTD snapshot has a grid spacing twice what is expressed in the input script.
 - The effective length of the BTD snapshot may be longer than anticipated from the input script because the grid spacing is different. Additionally, the number of grid points in the BTD snapshot is a multiple of ``<BTD>.buffer_size`` whereas the number of grid cells specified in the input deck may not be.
-- The code may require longer than anticipated to complete a BTD snapshot. The code starts filling the :math:`i^{th}` snapshot around step :math:`j_{\rm BTD start}={\rm ceil}\left( i\gamma(1-\beta)\frac{\Delta t_{\rm snapshot}}{\Delta t_{\rm boost}}\right)`. The code then saves information for one BTD cell every time step in the boosted frame simulation. The :math:`i^{th}` snapshot is completed and saved :math:`n_{z,{\rm snapshot}}=n_{\rm buffers}\cdot ({\rm buffer\ size})` time steps after it begins, which is when the effective snapshot length is covered by the simulation.
+- The code may require longer than anticipated to complete a BTD snapshot. The code starts filling the :math:`i^{th}` snapshot around step :math:`j_{\rm BTD start}={\rm ceil}\left( i\gamma(1-\beta\beta_{mw})\frac{\Delta t_{\rm snapshot}}{\Delta t_{\rm boost}}\right)`. The code then saves information for one BTD cell every time step in the boosted frame simulation. The :math:`i^{th}` snapshot is completed and saved :math:`n_{z,{\rm snapshot}}=n_{\rm buffers}\cdot ({\rm buffer\ size})` time steps after it begins, which is when the effective snapshot length is covered by the simulation.
 
 What kinds of RZ output do you support?
 ---------------------------------------

Source/Diagnostics/BTDiagnostics.H

@@ -161,6 +161,7 @@ private:
      *  in z-direction for both 2D and 3D simulations in the Cartesian frame of reference.
      */
     int m_moving_window_dir;
+    amrex::Real m_moving_window_beta;
 
     /** Number of back-transformed snapshots in the lab-frame requested by the user */
     int m_num_snapshots_lab = std::numeric_limits<int>::lowest();

Source/Diagnostics/BTDiagnostics.cpp

@@ -69,6 +69,7 @@ void BTDiagnostics::DerivedInitData ()
     m_gamma_boost = WarpX::gamma_boost;
     m_beta_boost = std::sqrt( 1._rt - 1._rt/( m_gamma_boost * m_gamma_boost) );
     m_moving_window_dir = WarpX::moving_window_dir;
+    m_moving_window_beta = WarpX::moving_window_v/PhysConst::c;
     // Currently, for BTD, all the data is averaged+coarsened to coarsest level
     // and then sliced+back-transformed+filled_to_buffer.
     // The number of levels to be output is nlev_output.
@@ -138,7 +139,7 @@ void BTDiagnostics::DerivedInitData ()
     const int lev = 0;
     const amrex::Real dt_boosted_frame = warpx.getdt(lev);
     const int moving_dir = WarpX::moving_window_dir;
-    const amrex::Real Lz_lab = warpx.Geom(lev).ProbLength(moving_dir) / WarpX::gamma_boost / (1._rt+WarpX::beta_boost);
+    const amrex::Real Lz_lab = warpx.Geom(lev).ProbLength(moving_dir) * WarpX::gamma_boost * (1._rt - WarpX::beta_boost*m_moving_window_beta);
     const int ref_ratio = 1;
     const amrex::Real dz_snapshot_grid = dz_lab(dt_boosted_frame, ref_ratio);
     // Need enough buffers so the snapshot length is longer than the lab frame length
@@ -149,22 +150,21 @@ void BTDiagnostics::DerivedInitData ()
     // the final snapshot starts filling when the
     // right edge of the moving window intersects the final snapshot
     // time of final snapshot : t_sn = t0 + i*dt_snapshot
-    // where t0 is the time of first BTD snapshot, t0 = zmax / c  * beta / (1-beta)
+    // where t0 is the time of first BTD snapshot, t0 = zmax / c  * beta / (1-beta*beta_mw)
     //
     // the right edge of the moving window at the time of the final snapshot
     // has space time coordinates
-    // time t_intersect = t_sn, position  z_intersect=zmax + c*t_sn
+    // time t_intersect = t_sn, position  z_intersect=zmax + v_mw*t_sn
     // the boosted time of this space time pair is
     // t_intersect_boost = gamma * (t_intersect - beta * z_intersect_boost/c)
-    //                   = gamma * (t_sn * (1 - beta) - beta * zmax / c)
-    //                   = gamma * (zmax*beta/c + i*dt_snapshot*(1-beta) - beta*zmax/c)
-    //                   = gamma * i * dt_snapshot * (1-beta)
-    //                   = i * dt_snapshot / gamma / (1+beta)
+    //                   = gamma * (t_sn * (1 - beta*beta_mw) - beta * zmax / c)
+    //                   = gamma * (zmax*beta/c + i*dt_snapshot*(1-beta*beta_mw) - beta*zmax/c)
+    //                   = gamma * (1-beta*beta_mw) * i * dt_snapshot
     //
     // if j = final snapshot starting step, then we want to solve
-    // j dt_boosted_frame >= t_intersect_boost = i * dt_snapshot / gamma / (1+beta)
-    // j >= i / gamma / (1+beta) * dt_snapshot / dt_boosted_frame
-    const int final_snapshot_starting_step = static_cast<int>(std::ceil(final_snapshot_iteration / WarpX::gamma_boost / (1._rt+WarpX::beta_boost) * m_dt_snapshots_lab / dt_boosted_frame));
+    // j dt_boosted_frame >= t_intersect_boost = i * gamma * (1-beta*beta_mw) * dt_snapshot
+    // j >= i * gamma * (1-beta*beta_mw) * dt_snapshot / dt_boosted_frame
+    const int final_snapshot_starting_step = static_cast<int>(std::ceil(final_snapshot_iteration * WarpX::gamma_boost * (1._rt - WarpX::beta_boost*m_moving_window_beta) * m_dt_snapshots_lab / dt_boosted_frame));
     const int final_snapshot_fill_iteration = final_snapshot_starting_step + num_buffers * m_buffer_size - 1;
     const amrex::Real final_snapshot_fill_time = final_snapshot_fill_iteration * dt_boosted_frame;
     if (WarpX::compute_max_step_from_btd) {
@@ -256,7 +256,7 @@ BTDiagnostics::ReadParameters ()
     bool snapshot_interval_is_specified = utils::parser::queryWithParser(
         pp_diag_name, "dt_snapshots_lab", m_dt_snapshots_lab);
     if ( utils::parser::queryWithParser(pp_diag_name, "dz_snapshots_lab", m_dz_snapshots_lab) ) {
-        m_dt_snapshots_lab = m_dz_snapshots_lab/PhysConst::c;
+        m_dt_snapshots_lab = m_dz_snapshots_lab/WarpX::moving_window_v;
         snapshot_interval_is_specified = true;
     }
     WARPX_ALWAYS_ASSERT_WITH_MESSAGE(snapshot_interval_is_specified,
@@ -338,13 +338,15 @@ BTDiagnostics::InitializeBufferData ( int i_buffer , int lev, bool restart)
     // When restarting boosted simulations, the code below needs to take
     // into account the fact that the position of the box at the beginning
     // of the simulation, is not the one that we had at t=0 (because of the moving window)
-    const amrex::Real boosted_moving_window_v = (WarpX::moving_window_v - m_beta_boost*PhysConst::c)
-                                        / (1._rt - m_beta_boost * WarpX::moving_window_v/PhysConst::c);
+    const amrex::Real boosted_moving_window_v = (m_moving_window_beta - m_beta_boost)
+                                        / (1._rt - m_beta_boost*m_moving_window_beta);
     // Lab-frame time for the i^th snapshot
     if (!restart) {
-        const amrex::Real zmax_0 = warpx.Geom(lev).ProbHi(m_moving_window_dir);
+        const amrex::Real zmax_boost = warpx.Geom(lev).ProbHi(m_moving_window_dir);
         m_t_lab.at(i_buffer) = m_intervals.GetBTDIteration(i_buffer) * m_dt_snapshots_lab
-            + m_gamma_boost*m_beta_boost*zmax_0/PhysConst::c;
+            + m_gamma_boost*m_beta_boost*zmax_boost/PhysConst::c;
+            // Note: gamma_boost*beta_boost*zmax_boost is equal to
+            // beta_boost*zmax_lab/(1-beta_boost*beta_moving_window)
     }
 
     // Define buffer domain in boosted frame at level, lev, with user-defined lo and hi
@@ -403,9 +405,9 @@ BTDiagnostics::InitializeBufferData ( int i_buffer , int lev, bool restart)
     // Define buffer_domain in lab-frame for the i^th snapshot.
     // Replace z-dimension with lab-frame co-ordinates.
     const amrex::Real zmin_buffer_lab = ( diag_dom.lo(m_moving_window_dir) - boosted_moving_window_v * warpx.gett_new(0) )
-                                / ( (1.0_rt + m_beta_boost) * m_gamma_boost);
+                                * (1.0_rt - m_beta_boost*m_moving_window_beta) * m_gamma_boost;
     const amrex::Real zmax_buffer_lab = ( diag_dom.hi(m_moving_window_dir) - boosted_moving_window_v * warpx.gett_new(0) )
-                                / ( (1.0_rt + m_beta_boost) * m_gamma_boost);
+                                * (1.0_rt - m_beta_boost*m_moving_window_beta) * m_gamma_boost;
 
     // Initialize buffer counter and z-positions of the  i^th snapshot in
     // boosted-frame and lab-frame