WarpXEvolve.cpp

/* Copyright 2019-2020 Andrew Myers, Ann Almgren, Aurore Blelly
 *                     Axel Huebl, Burlen Loring, David Grote
 *                     Glenn Richardson, Jean-Luc Vay, Luca Fedeli
 *                     Maxence Thevenet, Remi Lehe, Revathi Jambunathan
 *                     Weiqun Zhang, Yinjian Zhao
 *
 * This file is part of WarpX.
 *
 * License: BSD-3-Clause-LBNL
 */
#include "WarpX.H"

#include "BoundaryConditions/PML.H"
#include "Diagnostics/MultiDiagnostics.H"
#include "Diagnostics/ReducedDiags/MultiReducedDiags.H"
#include "EmbeddedBoundary/Enabled.H"
#include "Fields.H"
#include "FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H"
#ifdef WARPX_USE_FFT
#   ifdef WARPX_DIM_RZ
#       include "FieldSolver/SpectralSolver/SpectralSolverRZ.H"
#   else
#       include "FieldSolver/SpectralSolver/SpectralSolver.H"
#   endif
#endif
#include "FieldSolver/ImplicitSolvers/ImplicitSolver.H"
#include "Parallelization/GuardCellManager.H"
#include "Particles/MultiParticleContainer.H"
#include "Fluids/MultiFluidContainer.H"
#include "Fluids/WarpXFluidContainer.H"
#include "Particles/ParticleBoundaryBuffer.H"
#include "Python/callbacks.H"
#include "Utils/TextMsg.H"
#include "Utils/WarpXAlgorithmSelection.H"
#include "Utils/WarpXUtil.H"
#include "Utils/WarpXConst.H"

#include <ablastr/profiler/ProfilerWrapper.H>
#include <ablastr/utils/SignalHandling.H>
#include <ablastr/warn_manager/WarnManager.H>

#include <AMReX.H>
#include <AMReX_Array.H>
#include <AMReX_BLassert.H>
#include <AMReX_Geometry.H>
#include <AMReX_IntVect.H>
#include <AMReX_LayoutData.H>
#include <AMReX_MultiFab.H>
#include <AMReX_ParmParse.H>
#include <AMReX_Print.H>
#include <AMReX_REAL.H>
#include <AMReX_RealVect.H>
#include <AMReX_Utility.H>
#include <AMReX_Vector.H>

#include <algorithm>
#include <cmath>
#include <array>
#include <memory>
#include <ostream>
#include <string>
#include <vector>

using namespace amrex;
using ablastr::utils::SignalHandling;

namespace
{
    /** Print Unused Parameter Warnings after Step 1
     *
     * Instead of waiting for a simulation to end, we already do an early "unused parameter check"
     * after step 1 to inform users early of potential issues with their simulation setup.
     */
    void checkEarlyUnusedParams ()
    {
        amrex::Print() << "\n"; // better: conditional \n based on return value
        amrex::ParmParse::QueryUnusedInputs();

        // Print the warning list right after the first step.
        amrex::Print() << ablastr::warn_manager::GetWMInstance().PrintGlobalWarnings("FIRST STEP");
    }

    void StoreCurrent (int lev, ablastr::fields::MultiFabRegister& fields)
    {
        using ablastr::fields::Direction;
        using warpx::fields::FieldType;

        for (int idim = 0; idim < 3; ++idim) {
            const auto dir = Direction{idim};
            if (fields.has(FieldType::current_store, dir,lev)) {
                MultiFab::Copy(*fields.get(FieldType::current_store, dir, lev),
                               *fields.get(FieldType::current_fp, dir, lev),
                               0, 0, 1, fields.get(FieldType::current_store, dir, lev)->nGrowVect());
            }
        }
    }

    void RestoreCurrent (int lev, ablastr::fields::MultiFabRegister& fields)
    {
        using ablastr::fields::Direction;
        using warpx::fields::FieldType;

        for (int idim = 0; idim < 3; ++idim) {
            const auto dir = Direction{idim};
            if (fields.has(FieldType::current_store, dir, lev)) {
                std::swap(
                    *fields.get(FieldType::current_fp, dir, lev),
                    *fields.get(FieldType::current_store, dir, lev)
                );
            }
        }
    }
}

void
WarpX::SynchronizeVelocityWithPosition () {
    using ablastr::fields::Direction;
    using warpx::fields::FieldType;

    if (!m_is_synchronized) {
        // This assumes that the particle boundary conditions have been checked
        // so that the field gather in PushP will be correct.
        FillBoundaryE(guard_cells.ng_FieldGather);
        FillBoundaryB(guard_cells.ng_FieldGather);
        if (fft_do_time_averaging)
        {
            FillBoundaryE_avg(guard_cells.ng_FieldGather);
            FillBoundaryB_avg(guard_cells.ng_FieldGather);
        }
        UpdateAuxiliaryData();
        FillBoundaryAux(guard_cells.ng_UpdateAux);
        for (int lev = 0; lev <= finest_level; ++lev) {
            mypc->PushP(
                lev,
                0.5_rt*dt[lev],
                *m_fields.get(FieldType::Efield_aux, Direction{0}, lev),
                *m_fields.get(FieldType::Efield_aux, Direction{1}, lev),
                *m_fields.get(FieldType::Efield_aux, Direction{2}, lev),
                *m_fields.get(FieldType::Bfield_aux, Direction{0}, lev),
                *m_fields.get(FieldType::Bfield_aux, Direction{1}, lev),
                *m_fields.get(FieldType::Bfield_aux, Direction{2}, lev),
                MomentumPushType::Full
            );
        }
        m_is_synchronized = true;
    }
}

void
WarpX::Evolve (int numsteps)
{
    ABLASTR_PROFILE_REGION("WarpX::Evolve()");
    ABLASTR_PROFILE("WarpX::Evolve()");

    using ablastr::fields::Direction;

    Real cur_time = t_new[0];

    // Note that the default argument is numsteps = -1
    const int numsteps_max = (numsteps < 0)?(max_step):(istep[0] + numsteps);

    // check typos in inputs after step 1
    bool early_params_checked = false;

    static Real evolve_time = 0;

    const int step_begin = istep[0];
    for (int step = istep[0]; step < numsteps_max && cur_time < stop_time; ++step)
    {
        ABLASTR_PROFILE("WarpX::Evolve::step");
        const auto evolve_time_beg_step = static_cast<Real>(amrex::second());

        // Check and clear signal flags and asynchronously broadcast them from process 0
        SignalHandling::CheckSignals();

        multi_diags->NewIteration();

        bool verbose_step = (bool)verbose;
        if (verbose && m_limit_verbose_step) {

            int verbose_step_interval = 1;
            if (step<10) { verbose_step_interval = 1; }
            else if (step<100) { verbose_step_interval = 10; }
            else { verbose_step_interval = 100; }

            verbose_step = !((step+1)%verbose_step_interval);

        }

        // Start loop on time steps
        if (verbose_step) {
            amrex::Print() << "STEP " << step+1 << " starts ...\n";
        }
        ExecutePythonCallback("beforestep");

        CheckLoadBalance(step);

        // Update the timestep for solvers that support adaptive timestepping
        // (electrostatic and theta-implicit EM), provided const_dt is not specified.
        if (m_dt_update_interval.contains(step+1) || (step == 0 && m_max_dt.has_value())) {
            SynchronizeVelocityWithPosition();
            ApplyDtLimiters();
            if (verbose_step) {
                std::ostringstream oss;
                oss << "updating timestep to DT = " << std::scientific << std::setprecision(6) << dt[0];
                amrex::Print() << Utils::TextMsg::Info(oss.str());
            }
        }

        // If position and velocity are synchronized, push velocity backward one half step
        if (evolve_scheme == EvolveScheme::Explicit)
        {
            ExplicitFillBoundaryEBUpdateAux();
        }

        // If needed, deposit the initial ion charge and current densities that
        // will be used to update the E-field in Ohm's law.
        if (step == step_begin &&
            electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC
        ) {
            HybridPICInitializeRhoJandB();
        }

        // multi-physics: field ionization
        doFieldIonization();

#ifdef WARPX_QED
        // multi-physics: QED effects
        doQEDEvents();
        mypc->doQEDSchwinger();
#endif

        // perform particle injection
        ExecutePythonCallback("particleinjection");

        // perform collisions and advance fields and particles by one time step
        OneStep(cur_time, dt[0], step);

        // Resample particles
        // +1 is necessary here because value of step seen by user (first step is 1) is different than
        // value of step in code (first step is 0)
        mypc->doResampling(Geom(), istep[0]+1, verbose_step);

        if (evolve_scheme == EvolveScheme::Explicit) {
            applyMirrors(cur_time);
            // E : guard cells are NOT up-to-date
            // B : guard cells are NOT up-to-date
        }

        for (int lev = 0; lev <= max_level; ++lev) {
            ++istep[lev];
        }

        cur_time += dt[0];

        ShiftGalileanBoundary();

        // sync up time
        for (int i = 0; i <= max_level; ++i) {
            t_old[i] = t_new[i];
            t_new[i] = cur_time;
        }
        multi_diags->FilterComputePackFlush( step, false, true );

        const bool move_j = m_is_synchronized;
        // If m_is_synchronized we need to shift j too so that next step we can evolve E by dt/2.
        // We might need to move j because we are going to make a plotfile.
        const int num_moved = MoveWindow(step+1, move_j);

        // Update the accelerator lattice element finder if the window has moved,
        // from either a moving window or a boosted frame
        if (num_moved != 0 || gamma_boost > 1) {
            for (int lev = 0; lev <= finest_level; ++lev) {
                m_accelerator_lattice[lev]->UpdateElementFinder(lev, gett_new());
            }
        }

        HandleParticlesAtBoundaries(step, cur_time, num_moved);

        // Apply particle thermalizer (no-op until implemented)
        if (m_particle_thermalizer.defined()) {
            m_particle_thermalizer.applyThermalizer(*mypc);
        }

        if (m_implicit_solver) {
            ExecutePythonCallback("beforecollisions");
            mypc->doCollisions(step, cur_time, dt[0]);
            ExecutePythonCallback("aftercollisions");
        }

        // Electrostatic field solve step for electrostatic or Darwin solvers
        if( electrostatic_solver_id != ElectrostaticSolverAlgo::None )
        {
            ExecutePythonCallback("beforeEsolve");

            // Electrostatic solver:
            // The E-field is always reset to hold just the electrostatic component
            bool const reset_E_field = true;
            // The B-field is also reset unless the Darwin solver is used
            bool const reset_B_field = (evolve_scheme != EvolveScheme::Semi_Implicit_Darwin);

            // For each species: deposit charge and add the associated space-charge
            // E and B field to the grid ; this is done at the end of the PIC
            // loop (i.e. immediately after a `Redistribute` and before particle
            // positions are next pushed) so that the particles do not deposit out of bounds
            // and so that the fields are at the correct time in the output.
            ComputeSpaceChargeField( reset_E_field, reset_B_field );
            if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic) {
                // Call Magnetostatic Solver to solve for the vector potential A and compute the
                // B field.  Time varying A contribution to E field is neglected.
                // This is currently a lab frame calculation.
                ComputeMagnetostaticField();
            }

            // The external fields are added back on to the fine patch fields
            // (which were overwritten by electrostatic / magnetostatic solvers)
            // so that the net fields are the sum of the field solutions and any
            // external fields.
            // This is skipped for Darwin since in that case the "external" fields
            // are just treated as initial conditions (as for other EM solvers).
            if (evolve_scheme != EvolveScheme::Semi_Implicit_Darwin) {
                for (int lev = 0; lev <= max_level; ++lev) {
                    AddExternalFields(lev);
                }
            }
            ExecutePythonCallback("afterEsolve");
        }

        // Hybrid-PIC case. With the theta-implicit hybrid evolve scheme the
        // fields are already advanced self-consistently with the particles in
        // the implicit solver's OneStep, so this explicit field update is skipped.
        if (electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC &&
            !m_implicit_solver) {
            ExecutePythonCallback("beforeEsolve");
            // The particles are now at p^{n+1/2} and x^{n+1}. The fields
            // are updated according to the hybrid-PIC scheme (Ohm's law
            // and Ampere's law).
            HybridPICEvolveFields();
            ExecutePythonCallback("afterEsolve");
        }

        bool const do_diagnostic = (multi_diags->DoComputeAndPack(step) || reduced_diags->DoDiags(step));
        bool const end_of_step_loop = (step == numsteps_max - 1) || (cur_time + dt[0] >= stop_time - 1.e-3*dt[0]);
        if (synchronize_velocity_for_diagnostics &&
            (do_diagnostic || end_of_step_loop)) {
            // When the diagnostics require synchronization, push p by 0.5*dt to synchronize.
            // Note that this will be undone at the start of the next step by the half v-push
            // backwards.
            SynchronizeVelocityWithPosition();
        }

        // afterstep callback runs with the updated global time. It is included
        // in the evolve timing.
        ExecutePythonCallback("afterstep");

        /// reduced diags
        if (reduced_diags->m_plot_rd != 0)
        {
            reduced_diags->LoadBalance();
            reduced_diags->ComputeDiags(step);
            reduced_diags->WriteToFile(step);
        }
        multi_diags->FilterComputePackFlush( step );

        // execute afterdiagnostic callbacks
        ExecutePythonCallback("afterdiagnostics");

        // inputs: unused parameters (e.g. typos) check after step 1 has finished
        if (!early_params_checked) {
            ::checkEarlyUnusedParams();
            early_params_checked = true;
        }

        // create ending time stamp for calculating elapsed time each iteration
        const auto evolve_time_end_step = static_cast<Real>(amrex::second());
        evolve_time += evolve_time_end_step - evolve_time_beg_step;

        HandleSignals();

        if (verbose_step) {
            amrex::Print()<< "STEP " << step+1 << " ends." << " TIME = " << cur_time
                        << " DT = " << dt[0] << "\n";
            amrex::Print()<< "Evolve time = " << evolve_time
                      << " s; This step = " << evolve_time_end_step-evolve_time_beg_step
                      << " s; Avg. per step = " << evolve_time/(step-step_begin+1) << " s\n\n";
        }

        if (checkStopSimulation(cur_time)) {
            break;
        }
    } // End loop on time steps

    // This if statement is needed for PICMI, which allows the Evolve routine to be
    // called multiple times, otherwise diagnostics will be done at every call,
    // regardless of the diagnostic period parameter provided in the inputs.
    bool const final_time_step = (istep[0] == max_step)
                                || (cur_time >= stop_time - 1.e-3*dt[0]
                                 && cur_time < stop_time + dt[0]);
    if (final_time_step || m_exit_loop_due_to_interrupt_signal) {
        multi_diags->FilterComputePackFlushLastTimestep( istep[0] );
        if (m_exit_loop_due_to_interrupt_signal) { ExecutePythonCallback("onbreaksignal"); }
    }

    amrex::Print() <<
        ablastr::warn_manager::GetWMInstance().PrintGlobalWarnings("THE END");
}

void WarpX::OneStep (
    amrex::Real a_cur_time,
    amrex::Real a_dt,
    int a_step
)
{
    ABLASTR_PROFILE("WarpX::OneStep()");

    // implicit solver
    if (m_implicit_solver) {
        // advance fields and particles by one time step
        const int exit_status = m_implicit_solver->OneStep(a_cur_time, a_dt, a_step);
        if (exit_status < 0) {
            std::stringstream solverMsg;
            solverMsg << "ImplicitSolver::OneStep() failed at step = " << a_step
                      << " using dt = " << a_dt << ".\n"
                      << "Nonlinear solver failed to converge: exit status = " << exit_status;
            WARPX_ABORT_WITH_MESSAGE(solverMsg.str());
        }
    }
    // explicit solver
    else {
        // electrostatic solver or hybrid solver
        if (electromagnetic_solver_id == ElectromagneticSolverAlgo::None ||
            electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC) {
            // with collisions placed in the middle of the momentum push
            if (m_collisions_split_momentum_push) {
                // push particles (half momentum)
                PushParticlesandDeposit(
                    a_cur_time,
                    /*skip_deposition=*/true,
                    PositionPushType::None,
                    MomentumPushType::FirstHalf
                );

                // perform particle collisions
                ExecutePythonCallback("beforecollisions");
                mypc->doCollisions(a_step, a_cur_time, a_dt);
                ExecutePythonCallback("aftercollisions");

                // push particles (full position and half momentum)
                PushParticlesandDeposit(
                    a_cur_time,
                    /*skip_deposition=*/true,
                    PositionPushType::Full,
                    MomentumPushType::SecondHalf
                );
            }
            // with collisions placed before the position and momentum push, or without collisions
            else {
                // perform particle collisions
                ExecutePythonCallback("beforecollisions");
                mypc->doCollisions(a_step, a_cur_time, a_dt);
                ExecutePythonCallback("aftercollisions");

                // push particles (full position and full momentum)
                PushParticlesandDeposit(
                    a_cur_time,
                    /*skip_deposition=*/true,
                    PositionPushType::Full,
                    MomentumPushType::Full
                );
            }
        }
        // electromagnetic solver
        else {
            // without mesh refinement
            if (finest_level == 0) {
                // standard PIC loop
                if (!m_JRhom) {
                    OneStep_nosub(a_cur_time, a_dt, a_step);
                }
                // JRhom PIC loop
                else {
                    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
                        m_collisions_split_momentum_push == 0,
                        "Collisions with split momentum push not yet implemented for JRhom PIC loop."
                        "Set `collisions.split_momentum_push=0` to use JRhom with standard (pre-v-push collisions placement) collisions model."
                    );
                    // perform particle collisions
                    ExecutePythonCallback("beforecollisions");
                    mypc->doCollisions(a_step, a_cur_time, a_dt);
                    ExecutePythonCallback("aftercollisions");

                    OneStep_JRhom(a_cur_time);
                }
            }
            // with mesh refinement
            else {
                // without subcycling
                if (!m_do_subcycling) {
                    OneStep_nosub(a_cur_time, a_dt, a_step);
                }
                // with subcycling
                else {
                    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
                        finest_level == 1,
                        "Subcycling not implemented with more than 1 mesh refinement level"
                    );
                    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
                        m_collisions_split_momentum_push == 0,
                        "Collisions with split momentum push not yet implemented with subcycling."
                        "Set `collisions.split_momentum_push=0` to use subcycling with standard (pre-v-push collisions placement) collisions model."
                    );
                    // perform particle collisions
                    ExecutePythonCallback("beforecollisions");
                    mypc->doCollisions(a_step, a_cur_time, a_dt);
                    ExecutePythonCallback("aftercollisions");

                    OneStep_sub1(a_cur_time);
                }
            }
        }
    }
}

/**
 * \brief Perform one PIC iteration, without subcycling
 * i.e. all levels/patches use the same timestep (that of the finest level)
 * for the field advance and particle pusher.
 */
void
WarpX::OneStep_nosub (
    amrex::Real a_cur_time,
    amrex::Real a_dt,
    int a_step
)
{
    ABLASTR_PROFILE("WarpX::OneStep_nosub()");

    // Push particle from x^{n} to x^{n+1}
    //               from p^{n-1/2} to p^{n+1/2}
    // Deposit current j^{n+1/2}
    // Deposit charge density rho^{n}

    ExecutePythonCallback("beforedeposition");

    // with collisions placed in the middle of the momentum push
    if (m_collisions_split_momentum_push) {
        // push particles (half momentum)
        PushParticlesandDeposit(
            a_cur_time,
            /*skip_deposition=*/true,
            PositionPushType::None,
            MomentumPushType::FirstHalf
        );
        // perform particle collisions
        ExecutePythonCallback("beforecollisions");
        mypc->doCollisions(a_step, a_cur_time, a_dt);
        ExecutePythonCallback("aftercollisions");

        // push particles (full position and half momentum)
        PushParticlesandDeposit(
            a_cur_time,
            /*skip_deposition=*/false,
            PositionPushType::Full,
            MomentumPushType::SecondHalf
        );
    }
    else {
        // perform particle collisions
        ExecutePythonCallback("beforecollisions");
        mypc->doCollisions(a_step, a_cur_time, a_dt);
        ExecutePythonCallback("aftercollisions");

        // push particles (full position and full momentum)
        PushParticlesandDeposit(
            a_cur_time,
            /*skip_deposition=*/false,
            PositionPushType::Full,
            MomentumPushType::Full
        );
    }

    ExecutePythonCallback("afterdeposition");

    // Synchronize J and rho:
    // filter (if used), exchange guard cells, interpolate across MR levels
    // and apply boundary conditions
    SyncCurrentAndRho();

    // At this point, J is up-to-date inside the domain, and E and B are
    // up-to-date including enough guard cells for first step of the field
    // solve.

    // For extended PML: copy J from regular grid to PML, and damp J in PML
    if (do_pml && pml_has_particles) { CopyJPML(); }
    if (do_pml && do_pml_j_damping) { DampJPML(); }

    ExecutePythonCallback("beforeEsolve");

    // Push E and B from {n} to {n+1}
    // (And update guard cells immediately afterwards)
    if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) {
        if (use_hybrid_QED)
        {
            WarpX::Hybrid_QED_Push(dt);
            FillBoundaryE(guard_cells.ng_alloc_EB);
        }
        PushPSATD(a_cur_time);

        if (do_pml) {
            DampPML();
        }

        if (use_hybrid_QED) {
            FillBoundaryE(guard_cells.ng_alloc_EB);
            FillBoundaryB(guard_cells.ng_alloc_EB, WarpX::sync_nodal_points);
            WarpX::Hybrid_QED_Push(dt);
            FillBoundaryE(guard_cells.ng_afterPushPSATD, WarpX::sync_nodal_points);
        }
        else {
            FillBoundaryE(guard_cells.ng_afterPushPSATD, WarpX::sync_nodal_points);
            FillBoundaryB(guard_cells.ng_afterPushPSATD, WarpX::sync_nodal_points);
            if (WarpX::do_dive_cleaning || WarpX::do_pml_dive_cleaning) {
                FillBoundaryF(guard_cells.ng_alloc_F, WarpX::sync_nodal_points);
            }
            if (WarpX::do_divb_cleaning || WarpX::do_pml_divb_cleaning) {
                FillBoundaryG(guard_cells.ng_alloc_G, WarpX::sync_nodal_points);
            }
        }
    } else {
        EvolveF(0.5_rt * dt[0], /*rho_comp=*/0);
        EvolveG(0.5_rt * dt[0]);
        FillBoundaryF(guard_cells.ng_FieldSolverF);
        FillBoundaryG(guard_cells.ng_FieldSolverG);

        EvolveB(0.5_rt * dt[0], SubcyclingHalf::FirstHalf, a_cur_time); // We now have B^{n+1/2}
        FillBoundaryB(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);

        if (m_em_solver_medium == MediumForEM::Vacuum) {
            // vacuum medium
            EvolveE(dt[0], a_cur_time); // We now have E^{n+1}
        } else if (m_em_solver_medium == MediumForEM::Macroscopic) {
            // macroscopic medium
            MacroscopicEvolveE(dt[0], a_cur_time); // We now have E^{n+1}
        } else {
            WARPX_ABORT_WITH_MESSAGE("Medium for EM is unknown");
        }
        FillBoundaryE(guard_cells.ng_FieldSolver, WarpX::sync_nodal_points);

        EvolveF(0.5_rt * dt[0], /*rho_comp=*/1);
        EvolveG(0.5_rt * dt[0]);
        EvolveB(0.5_rt * dt[0], SubcyclingHalf::SecondHalf, a_cur_time + 0.5_rt * dt[0]); // We now have B^{n+1}

        if (do_pml) {
            DampPML();
            FillBoundaryE(guard_cells.ng_MovingWindow, WarpX::sync_nodal_points);
            FillBoundaryB(guard_cells.ng_MovingWindow, WarpX::sync_nodal_points);
            FillBoundaryF(guard_cells.ng_MovingWindow, WarpX::sync_nodal_points);
            FillBoundaryG(guard_cells.ng_MovingWindow, WarpX::sync_nodal_points);
        }

        // E and B are up-to-date in the domain, but all guard cells are
        // outdated.
        if (m_safe_guard_cells) {
            FillBoundaryB(guard_cells.ng_alloc_EB);
        }
    } // !PSATD

    ExecutePythonCallback("afterEsolve");
}

bool WarpX::checkStopSimulation (amrex::Real cur_time)
{
    m_exit_loop_due_to_interrupt_signal = SignalHandling::TestAndResetActionRequestFlag(SignalHandling::SIGNAL_REQUESTS_BREAK);
    return (cur_time >= stop_time - 1.e-3*dt[0])  ||
        m_exit_loop_due_to_interrupt_signal;
}

void WarpX::ExplicitFillBoundaryEBUpdateAux ()
{
    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(evolve_scheme == EvolveScheme::Explicit,
        "Cannot call WarpX::ExplicitFillBoundaryEBUpdateAux without Explicit evolve scheme set!");

    using ablastr::fields::Direction;
    using warpx::fields::FieldType;

    // At the beginning, we have B^{n} and E^{n}.
    // Particles have p^{n} and x^{n}.
    // m_is_synchronized is true.

    if (m_is_synchronized) {
        // Not called at each iteration, so exchange all guard cells
        FillBoundaryE(guard_cells.ng_alloc_EB);
        FillBoundaryB(guard_cells.ng_alloc_EB);

        UpdateAuxiliaryData();
        FillBoundaryAux(guard_cells.ng_UpdateAux);
        // on first step, push p by -0.5*dt
        for (int lev = 0; lev <= finest_level; ++lev)
        {
            mypc->PushP(
                lev,
                -0.5_rt*dt[lev],
                *m_fields.get(FieldType::Efield_aux, Direction{0}, lev),
                *m_fields.get(FieldType::Efield_aux, Direction{1}, lev),
                *m_fields.get(FieldType::Efield_aux, Direction{2}, lev),
                *m_fields.get(FieldType::Bfield_aux, Direction{0}, lev),
                *m_fields.get(FieldType::Bfield_aux, Direction{1}, lev),
                *m_fields.get(FieldType::Bfield_aux, Direction{2}, lev),
                MomentumPushType::Full
            );
        }
        m_is_synchronized = false;

    } else {
        // Beyond one step, we have E^{n} and B^{n}.
        // Particles have p^{n-1/2} and x^{n}.
        // E and B: enough guard cells to update Aux or call Field Gather in fp and cp
        // Need to update Aux on lower levels, to interpolate to higher levels.

        // E and B are up-to-date inside the domain only
        FillBoundaryE(guard_cells.ng_FieldGather);
        FillBoundaryB(guard_cells.ng_FieldGather);
        if (electrostatic_solver_id == ElectrostaticSolverAlgo::None) {
            if (fft_do_time_averaging)
            {
                FillBoundaryE_avg(guard_cells.ng_FieldGather);
                FillBoundaryB_avg(guard_cells.ng_FieldGather);
            }
            // TODO Remove call to FillBoundaryAux before UpdateAuxiliaryData?
            if (WarpX::electromagnetic_solver_id != ElectromagneticSolverAlgo::PSATD) {
                FillBoundaryAux(guard_cells.ng_UpdateAux);
            }
        }
        UpdateAuxiliaryData();
        FillBoundaryAux(guard_cells.ng_UpdateAux);
    }
}

void WarpX::HandleParticlesAtBoundaries (int step, amrex::Real cur_time, int num_moved)
{
    mypc->ContinuousFluxInjection(cur_time, dt[0]);

    ExecutePythonCallback("particlescraper");

    mypc->ApplyBoundaryConditions();
    m_particle_boundary_buffer->gatherParticlesFromDomainBoundaries(*mypc, cur_time);

    // Without mesh refinement, use a local redistribute when particles can only
    // have moved by a small number of cells; otherwise fall back to a global one.
    if (finest_level == 0) {
        // Estimate, per direction, the maximum distance a particle may have
        // travelled during this step, expressed in number of cells.
        // (Geom().CellSizeArray() is indexed by active dimension 0..SPACEDIM-1.)
        const amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> dx = Geom(0).CellSizeArray();

        // Particles cannot travel faster than the speed of light, so c * dt / dx
        // is a physical upper bound on the number of cells crossed per direction.
        amrex::RealVect max_distance_relative_to_grid;
        for (int d = 0; d < AMREX_SPACEDIM; ++d) {
            max_distance_relative_to_grid[d] = PhysConst::c * dt[0] / dx[d];
        }

        // Moving window: particles can additionally move by the number of cells
        // that the window was shifted, along the moving-window direction.
        if (moving_window_dir >= 0) {
            max_distance_relative_to_grid[moving_window_dir] += static_cast<amrex::Real>(num_moved);
        }


        // Galilean algorithm: account for the extra grid shift due to the moving
        // Galilean frame. m_v_galilean is indexed by x/y/z, so map its components
        // onto the active simulation dimensions.
#if defined(WARPX_DIM_3D)
        const amrex::RealVect v_galilean = {m_v_galilean[0], m_v_galilean[1], m_v_galilean[2]};
#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
        const amrex::RealVect v_galilean = {m_v_galilean[0], m_v_galilean[2]};
#elif defined(WARPX_DIM_1D_Z)
        const amrex::RealVect v_galilean(m_v_galilean[2]);
#else // WARPX_DIM_RCYLINDER, WARPX_DIM_RSPHERE: no Galilean shift
        const amrex::RealVect v_galilean = amrex::RealVect::TheZeroVector();
#endif
        for (int d = 0; d < AMREX_SPACEDIM; ++d) {
            max_distance_relative_to_grid[d] += std::abs(v_galilean[d]) * dt[0] / dx[d];
        }

        // Convert to an integer number of cells (rounding up), per direction.
        amrex::IntVect max_cells_travelled;
        for (int d = 0; d < AMREX_SPACEDIM; ++d) {
            max_cells_travelled[d] =
                static_cast<int>(std::ceil(max_distance_relative_to_grid[d]));
        }

        // If, in any direction, max_cells_travelled reaches the domain size, the
        // local search is no longer more efficient than (and may crash in lieu
        // of) a full redistribute, so fall back in that case.
        const amrex::IntVect domain_length = Geom(0).Domain().length();
        bool use_local_redistribute = true;
        for (int d = 0; d < AMREX_SPACEDIM; ++d) {
            if (max_cells_travelled[d] >= domain_length[d]) { use_local_redistribute = false; }
        }
        if (use_local_redistribute) {
            mypc->RedistributeLocal(max_cells_travelled);
        } else {
            mypc->Redistribute();
        }
    }
    else {
        mypc->Redistribute();
    }

    // interact the particles with EB walls (if present)
    if (EB::enabled()) {
        using warpx::fields::FieldType;
        mypc->ScrapeParticlesAtEB(m_fields.get_mr_levels(FieldType::distance_to_eb, finest_level));
        m_particle_boundary_buffer->gatherParticlesFromEmbeddedBoundaries(
            *mypc, m_fields.get_mr_levels(FieldType::distance_to_eb, finest_level), cur_time);
        if (eb_particle_boundary == ParticleBoundaryType::Absorbing) {
            // If particles are simply absorbed, no need for a full Redistribute.
            // Instead: simply delete the absorbed particles
            mypc->deleteInvalidParticles();
        } else {
            // For other particle boundary conditions (e.g. reflecting),
            // particles can move to a different sub-domain, so we need a full Redistribute
            mypc->Redistribute();
        }
    }

    if (sort_intervals.contains(step+1)) {
        if (verbose && !m_limit_verbose_step) {
            amrex::Print() << Utils::TextMsg::Info("re-sorting particles");
        }
        mypc->SortParticlesByBin(
            sort_bin_size, m_sort_particles_for_deposition, m_sort_idx_type);
    }
}

void WarpX::SyncCurrentAndRho ()
{
    using ablastr::fields::Direction;
    using warpx::fields::FieldType;

    if (electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD)
    {
        if (fft_periodic_single_box)
        {
            // With periodic single box, synchronize J and rho here,
            // even with current correction or Vay deposition
            std::string const current_fp_string = (current_deposition_algo == CurrentDepositionAlgo::Vay)
                ? "current_fp_vay" : "current_fp";
            // TODO Replace current_cp with current_cp_vay once Vay deposition is implemented with MR

            SyncCurrent(current_fp_string);
            SyncRho();

        }
        else // no periodic single box
        {
            // Without periodic single box, synchronize J and rho here,
            // except with current correction or Vay deposition:
            // in these cases, synchronize later (in WarpX::PushPSATD)
            if (!current_correction &&
                current_deposition_algo != CurrentDepositionAlgo::Vay)
            {
                SyncCurrent("current_fp");
                SyncRho();
            }

            if (current_deposition_algo == CurrentDepositionAlgo::Vay)
            {
                // TODO This works only without mesh refinement
                const int lev = 0;
                if (use_filter) {
                    ApplyFilterJ(m_fields.get_mr_levels_alldirs(FieldType::current_fp_vay, finest_level), lev);
                }
            }
        }
    }
    else // FDTD
    {
        SyncCurrent("current_fp");
        SyncRho();
    }

    // Reflect charge and current density over PEC boundaries, if needed.
    for (int lev = 0; lev <= finest_level; ++lev)
    {
        if (m_fields.has(FieldType::rho_fp, lev)) {
            ApplyRhofieldBoundary(lev, m_fields.get(FieldType::rho_fp,lev), PatchType::fine);
        }
        ApplyJfieldBoundary(lev,
            m_fields.get(FieldType::current_fp, Direction{0}, lev),
            m_fields.get(FieldType::current_fp, Direction{1}, lev),
            m_fields.get(FieldType::current_fp, Direction{2}, lev),
            PatchType::fine);
        if (lev > 0) {
            if (m_fields.has(FieldType::rho_cp, lev)) {
                ApplyRhofieldBoundary(lev, m_fields.get(FieldType::rho_cp,lev), PatchType::coarse);
            }
            ApplyJfieldBoundary(lev,
                m_fields.get(FieldType::current_cp, Direction{0}, lev),
                m_fields.get(FieldType::current_cp, Direction{1}, lev),
                m_fields.get(FieldType::current_cp, Direction{2}, lev),
                PatchType::coarse);
        }
    }
}

void
WarpX::OneStep_JRhom (const amrex::Real cur_time)
{
#ifdef WARPX_USE_FFT

    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
        WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD,
        "JRhom algorithm not implemented with the FDTD solver"
    );

    using warpx::fields::FieldType;

    bool const skip_lev0_coarse_patch = true;

    const int rho_mid = spectral_solver_fp[0]->m_spectral_index.rho_mid;
    const int rho_new = spectral_solver_fp[0]->m_spectral_index.rho_new;

    // Push particle from x^{n} to x^{n+1}
    //               from p^{n-1/2} to p^{n+1/2}
    const bool skip_deposition = true;
    PushParticlesandDeposit(cur_time, skip_deposition);

    // Initialize PSATD-JRhom loop:

    // 1) Prepare E,B,F,G fields in spectral space
    PSATDForwardTransformEB();
    if (WarpX::do_dive_cleaning) { PSATDForwardTransformF(); }
    if (WarpX::do_divb_cleaning) { PSATDForwardTransformG(); }

    // 2) Set the averaged fields to zero
    if (WarpX::fft_do_time_averaging) { PSATDEraseAverageFields(); }

    // 3) Deposit rho (in rho_new, since it will be moved during the loop)
    //    (after checking that pointer to rho_fp on MR level 0 is not null)
    if (m_fields.has(FieldType::rho_fp, 0) && time_dependency_rho != TimeDependencyRho::Constant)
    {
        ablastr::fields::MultiLevelScalarField const rho_fp = m_fields.get_mr_levels(FieldType::rho_fp, finest_level);

        std::string const rho_fp_string = "rho_fp";
        std::string const rho_cp_string = "rho_cp";

        // Deposit rho at relative time -dt
        // (dt[0] denotes the time step on mesh refinement level 0)
        mypc->DepositCharge(rho_fp, -dt[0]);
        // Filter, exchange boundary, and interpolate across levels
        SyncRho();
        // Forward FFT of rho
        PSATDForwardTransformRho(rho_fp_string, rho_cp_string, 0, rho_new);
    }

    // 4) Deposit J at relative time -dt with time step dt
    //    (dt[0] denotes the time step on mesh refinement level 0)
    if (time_dependency_J != TimeDependencyJ::Constant)
    {
        std::string const current_string = (do_current_centering) ? "current_fp_nodal" : "current_fp";
        mypc->DepositCurrent( m_fields.get_mr_levels_alldirs(current_string, finest_level), dt[0], -dt[0]);
        // Synchronize J: filter, exchange boundary, and interpolate across levels.
        // With current centering, the nodal current is deposited in 'current',
        // namely 'current_fp_nodal': SyncCurrent stores the result of its centering
        // into 'current_fp' and then performs both filtering, if used, and exchange
        // of guard cells.
        SyncCurrent("current_fp");
        // Forward FFT of J
        PSATDForwardTransformJ("current_fp", "current_cp");
    }

    // Number of depositions for multi-J scheme
    const int n_deposit = WarpX::m_JRhom_subintervals;
    // Time sub-step for each multi-J deposition
    const amrex::Real sub_dt = dt[0] / static_cast<amrex::Real>(n_deposit);
    // Whether to perform PSATD-JRhom depositions on a time interval that spans
    // one or two full time steps (from n*dt to (n+1)*dt, or from n*dt to (n+2)*dt)
    const int n_loop = (WarpX::fft_do_time_averaging) ? 2*n_deposit : n_deposit;

    // Loop over PSATD-JRhom depositions
    for (int i_deposit = 0; i_deposit < n_loop; i_deposit++)
    {
        // Move J from new to old if J is linear or quadratic in time
        if (time_dependency_J != TimeDependencyJ::Constant) { PSATDMoveJNewToJOld(); }

        const amrex::Real t_deposit_current = (time_dependency_J == TimeDependencyJ::Linear) ?
            (i_deposit-n_deposit+1)*sub_dt : (i_deposit-n_deposit+0.5_rt)*sub_dt;

        const amrex::Real t_deposit_charge = (time_dependency_rho == TimeDependencyRho::Linear) ?
            (i_deposit-n_deposit+1)*sub_dt : (i_deposit-n_deposit+0.5_rt)*sub_dt;

        // Deposit new J at relative time t_deposit_current with time step dt
        // (dt[0] denotes the time step on mesh refinement level 0)
        std::string const current_string = (do_current_centering) ? "current_fp_nodal" : "current_fp";
        mypc->DepositCurrent( m_fields.get_mr_levels_alldirs(current_string, finest_level), dt[0], t_deposit_current);
        // Synchronize J: filter, exchange boundary, and interpolate across levels.
        // With current centering, the nodal current is deposited in 'current',
        // namely 'current_fp_nodal': SyncCurrent stores the result of its centering
        // into 'current_fp' and then performs both filtering, if used, and exchange
        // of guard cells.
        SyncCurrent("current_fp");
        // Forward FFT of J
        PSATDForwardTransformJ("current_fp", "current_cp");

        if (time_dependency_J == TimeDependencyJ::Quadratic)
        {
            PSATDMoveJNewToJMid();
            mypc->DepositCurrent( m_fields.get_mr_levels_alldirs(current_string, finest_level),  dt[0], t_deposit_current + 0.5_rt*sub_dt);
            SyncCurrent("current_fp");
            PSATDForwardTransformJ("current_fp", "current_cp");
        }

        // Deposit new rho
        // (after checking that pointer to rho_fp on MR level 0 is not null)
        if (m_fields.has(FieldType::rho_fp, 0))
        {
            ablastr::fields::MultiLevelScalarField const rho_fp = m_fields.get_mr_levels(FieldType::rho_fp, finest_level);

            std::string const rho_fp_string = "rho_fp";
            std::string const rho_cp_string = "rho_cp";

            // Move rho from new to old if rho is linear in time
            if (time_dependency_rho != TimeDependencyRho::Constant) { PSATDMoveRhoNewToRhoOld(); }

            // Deposit rho at relative time t_deposit_charge
            mypc->DepositCharge(rho_fp, t_deposit_charge);
            // Filter, exchange boundary, and interpolate across levels
            SyncRho();
            // Forward FFT of rho
            const int rho_idx = (time_dependency_rho != TimeDependencyRho::Constant) ? rho_new : rho_mid;
            PSATDForwardTransformRho(rho_fp_string, rho_cp_string, 0, rho_idx);

            if (time_dependency_rho == TimeDependencyRho::Quadratic)
            {
                PSATDMoveRhoNewToRhoMid();
                mypc->DepositCharge(rho_fp, t_deposit_charge + 0.5_rt*sub_dt);
                SyncRho();
                PSATDForwardTransformRho(rho_fp_string, rho_cp_string, 0, rho_new);
            }
        }

        if (WarpX::current_correction)
        {
            WARPX_ABORT_WITH_MESSAGE(
                "Current correction not implemented for PSATD-JRhom algorithm.");
        }

        // Advance E,B,F,G fields in time and update the average fields
        PSATDPushSpectralFields();

        // Transform non-average fields E,B,F,G after n_deposit pushes
        // (the relative time reached here coincides with an integer full time step)
        if (i_deposit == n_deposit-1)
        {
            PSATDBackwardTransformEB();
            if (WarpX::do_dive_cleaning) { PSATDBackwardTransformF(); }
            if (WarpX::do_divb_cleaning) { PSATDBackwardTransformG(); }
        }
    }

    // Transform fields back to real space
    if (WarpX::fft_do_time_averaging)
    {
        // We summed the integral of the field over 2*dt
        PSATDScaleAverageFields(1._rt / (2._rt*dt[0]));
        PSATDBackwardTransformEBavg(
            m_fields.get_mr_levels_alldirs(FieldType::Efield_avg_fp, finest_level),
            m_fields.get_mr_levels_alldirs(FieldType::Bfield_avg_fp, finest_level),
            m_fields.get_mr_levels_alldirs(FieldType::Efield_avg_cp, finest_level, skip_lev0_coarse_patch),
            m_fields.get_mr_levels_alldirs(FieldType::Bfield_avg_cp, finest_level, skip_lev0_coarse_patch)
        );
    }

    // Evolve fields in PML
    for (int lev = 0; lev <= finest_level; ++lev)
    {
        if (do_pml && pml[lev]->ok())
        {
            pml[lev]->PushPSATD(m_fields, lev);
        }
        ApplyEfieldBoundary(lev, PatchType::fine, cur_time + dt[0]);
        if (lev > 0) { ApplyEfieldBoundary(lev, PatchType::coarse, cur_time + dt[0]); }
        ApplyBfieldBoundary(lev, PatchType::fine, SubcyclingHalf::FirstHalf, cur_time + dt[0]);
        if (lev > 0) { ApplyBfieldBoundary(lev, PatchType::coarse, SubcyclingHalf::FirstHalf, cur_time + dt[0]); }
    }

    // Damp fields in PML before exchanging guard cells
    if (do_pml)
    {
        DampPML();
    }

    // Exchange guard cells and synchronize nodal points
    FillBoundaryE(guard_cells.ng_alloc_EB, WarpX::sync_nodal_points);
    FillBoundaryB(guard_cells.ng_alloc_EB, WarpX::sync_nodal_points);
    if (WarpX::do_dive_cleaning || WarpX::do_pml_dive_cleaning) {
        FillBoundaryF(guard_cells.ng_alloc_F, WarpX::sync_nodal_points);
    }
    if (WarpX::do_divb_cleaning || WarpX::do_pml_divb_cleaning) {
        FillBoundaryG(guard_cells.ng_alloc_G, WarpX::sync_nodal_points);
    }

#else
    amrex::ignore_unused(cur_time);
    WARPX_ABORT_WITH_MESSAGE(
        "JRhom algorithm not implemented with the FDTD solver");
#endif // WARPX_USE_FFT
}

/**
 *  \brief Perform one PIC iteration, with subcycling
 *  i.e. The fine patch uses a smaller timestep (and steps more often)
 *  than the coarse patch, for the field advance and particle pusher.
 *
 * This version of subcycling only works for 2 levels and with a refinement
 * ratio of 2.
 * The particles and fields of the fine patch are pushed twice
 * (with dt[coarse]/2) in this routine.
 * The particles of the coarse patch and mother grid are pushed only once
 * (with dt[coarse]). The fields on the coarse patch and mother grid
 * are pushed in a way which is equivalent to pushing once only, with
 * a current which is the average of the coarse + fine current at the 2
 * steps of the fine grid.
 *
 */
void
WarpX::OneStep_sub1 (Real cur_time)
{
    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
        electrostatic_solver_id == ElectrostaticSolverAlgo::None,
        "Electrostatic solver cannot be used with sub-cycling."
    );

    // TODO: we could save some charge depositions

    // check that MR has exactly two levels and a refinement ratio of 2 in all directions when subcycling is used
    const int fine_lev = 1;
    const int coarse_lev = 0;
    const amrex::IntVect& crse_to_fine_ref_ratio = this->refRatio(coarse_lev);
    const bool ref_ratio_is_uniform_two = (crse_to_fine_ref_ratio == amrex::IntVect(2));
    if (finest_level != 1 || !ref_ratio_is_uniform_two) {
        std::ostringstream msg;
        msg << "MR with subcycling algorithm requires exactly two levels of MR "
            << "and a refinement ratio of 2 in all directions. "
            << "Found: finest_level = " << finest_level
            << ", refRatio(" << coarse_lev << ") = " << crse_to_fine_ref_ratio;
        WARPX_ABORT_WITH_MESSAGE(msg.str());
    }

    using warpx::fields::FieldType;

    bool const skip_lev0_coarse_patch = true;

    // i) Push particles and fields on the fine patch (first fine step)
    PushParticlesandDeposit(fine_lev, cur_time, SubcyclingHalf::FirstHalf);
    RestrictCurrentFromFineToCoarsePatch(
        m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level),
        m_fields.get_mr_levels_alldirs(FieldType::current_cp, finest_level, skip_lev0_coarse_patch), fine_lev);
    RestrictRhoFromFineToCoarsePatch(fine_lev);
    if (use_filter) {
        ApplyFilterJ( m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level), fine_lev);
    }
    SumBoundaryJ(
        m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level),
        fine_lev, Geom(fine_lev).periodicity());

    if (m_fields.has(FieldType::rho_fp, finest_level) &&
        m_fields.has(FieldType::rho_cp, finest_level)) {
        ApplyFilterandSumBoundaryRho(
            m_fields.get_mr_levels(FieldType::rho_fp, finest_level),
            m_fields.get_mr_levels(FieldType::rho_cp, finest_level, skip_lev0_coarse_patch),
            fine_lev, PatchType::fine, 0, 2*ncomps);
    }

    EvolveB(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], SubcyclingHalf::FirstHalf, cur_time);
    EvolveF(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], /*rho_comp=*/0);
    FillBoundaryB(fine_lev, PatchType::fine, guard_cells.ng_FieldSolver,
                  WarpX::sync_nodal_points);
    FillBoundaryF(fine_lev, PatchType::fine, guard_cells.ng_alloc_F,
                  WarpX::sync_nodal_points);

    EvolveE(fine_lev, PatchType::fine, dt[fine_lev], cur_time);
    FillBoundaryE(fine_lev, PatchType::fine, guard_cells.ng_FieldGather);

    EvolveB(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], SubcyclingHalf::SecondHalf, cur_time + 0.5_rt * dt[fine_lev]);
    EvolveF(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], /*rho_comp=*/1);

    if (do_pml) {
        FillBoundaryF(fine_lev, PatchType::fine, guard_cells.ng_alloc_F);
        DampPML(fine_lev, PatchType::fine);
        FillBoundaryE(fine_lev, PatchType::fine, guard_cells.ng_FieldGather);
    }

    FillBoundaryB(fine_lev, PatchType::fine, guard_cells.ng_FieldGather);

    // ii) Push particles on the coarse patch and mother grid.
    // Push the fields on the coarse patch and mother grid
    // by only half a coarse step (first half)
    PushParticlesandDeposit(coarse_lev, cur_time, SubcyclingHalf::None);
    ::StoreCurrent(coarse_lev, m_fields);
    AddCurrentFromFineLevelandSumBoundary(
        m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level),
        m_fields.get_mr_levels_alldirs(FieldType::current_cp, finest_level, skip_lev0_coarse_patch),
        m_fields.get_mr_levels_alldirs(FieldType::current_buf, finest_level, skip_lev0_coarse_patch), coarse_lev);

    if (m_fields.has(FieldType::rho_fp, finest_level) &&
        m_fields.has(FieldType::rho_cp, finest_level) &&
        m_fields.has(FieldType::rho_buf, finest_level)) {
        AddRhoFromFineLevelandSumBoundary(
            m_fields.get_mr_levels(FieldType::rho_fp, finest_level),
            m_fields.get_mr_levels(FieldType::rho_cp, finest_level, skip_lev0_coarse_patch),
            m_fields.get_mr_levels(FieldType::rho_buf, finest_level, skip_lev0_coarse_patch),
            coarse_lev, 0, ncomps);
    }

    EvolveB(fine_lev, PatchType::coarse, dt[fine_lev], SubcyclingHalf::FirstHalf, cur_time);
    EvolveF(fine_lev, PatchType::coarse, dt[fine_lev], /*rho_comp=*/0);
    FillBoundaryB(fine_lev, PatchType::coarse, guard_cells.ng_FieldGather);
    FillBoundaryF(fine_lev, PatchType::coarse, guard_cells.ng_FieldSolverF);

    EvolveE(fine_lev, PatchType::coarse, dt[fine_lev], cur_time);
    FillBoundaryE(fine_lev, PatchType::coarse, guard_cells.ng_FieldGather);

    EvolveB(coarse_lev, PatchType::fine, 0.5_rt*dt[coarse_lev], SubcyclingHalf::FirstHalf, cur_time);
    EvolveF(coarse_lev, PatchType::fine, 0.5_rt*dt[coarse_lev], /*rho_comp=*/0);
    FillBoundaryB(coarse_lev, PatchType::fine, guard_cells.ng_FieldGather,
                    WarpX::sync_nodal_points);
    FillBoundaryF(coarse_lev, PatchType::fine, guard_cells.ng_FieldSolverF,
                    WarpX::sync_nodal_points);

    EvolveE(coarse_lev, PatchType::fine, 0.5_rt*dt[coarse_lev], cur_time);
    FillBoundaryE(coarse_lev, PatchType::fine, guard_cells.ng_FieldGather);

    // TODO Remove call to FillBoundaryAux before UpdateAuxiliaryData?
    FillBoundaryAux(guard_cells.ng_UpdateAux);
    // iii) Get auxiliary fields on the fine grid, at dt[fine_lev]
    UpdateAuxiliaryData();
    FillBoundaryAux(guard_cells.ng_UpdateAux);

    // iv) Push particles and fields on the fine patch (second fine step)
    PushParticlesandDeposit(fine_lev, cur_time + dt[fine_lev], SubcyclingHalf::SecondHalf);
    RestrictCurrentFromFineToCoarsePatch(
        m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level),
        m_fields.get_mr_levels_alldirs(FieldType::current_cp, finest_level, skip_lev0_coarse_patch), fine_lev);
    RestrictRhoFromFineToCoarsePatch(fine_lev);
    if (use_filter) {
        ApplyFilterJ( m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level), fine_lev);
    }
    SumBoundaryJ( m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level), fine_lev, Geom(fine_lev).periodicity());

    if (m_fields.has(FieldType::rho_fp, finest_level) &&
        m_fields.has(FieldType::rho_cp, finest_level)) {
        ApplyFilterandSumBoundaryRho(
            m_fields.get_mr_levels(FieldType::rho_fp, finest_level),
            m_fields.get_mr_levels(FieldType::rho_cp, finest_level, skip_lev0_coarse_patch),
            fine_lev, PatchType::fine, 0, ncomps);
    }

    EvolveB(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], SubcyclingHalf::FirstHalf, cur_time + dt[fine_lev]);
    EvolveF(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], /*rho_comp=*/0);
    FillBoundaryB(fine_lev, PatchType::fine, guard_cells.ng_FieldSolver);
    FillBoundaryF(fine_lev, PatchType::fine, guard_cells.ng_FieldSolverF);

    EvolveE(fine_lev, PatchType::fine, dt[fine_lev], cur_time + dt[fine_lev]);
    FillBoundaryE(fine_lev, PatchType::fine, guard_cells.ng_FieldSolver,
                    WarpX::sync_nodal_points);

    EvolveB(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], SubcyclingHalf::SecondHalf, cur_time + 1.5_rt*dt[fine_lev]);
    EvolveF(fine_lev, PatchType::fine, 0.5_rt*dt[fine_lev], /*rho_comp=*/1);

    if (do_pml) {
        DampPML(fine_lev, PatchType::fine);
        FillBoundaryE(fine_lev, PatchType::fine, guard_cells.ng_FieldSolver);
    }

    if ( m_safe_guard_cells ) {
        FillBoundaryF(fine_lev, PatchType::fine, guard_cells.ng_FieldSolver);
    }
    FillBoundaryB(fine_lev, PatchType::fine, guard_cells.ng_FieldSolver);

    // v) Push the fields on the coarse patch and mother grid
    // by only half a coarse step (second half)
    ::RestoreCurrent(coarse_lev, m_fields);
    AddCurrentFromFineLevelandSumBoundary(
        m_fields.get_mr_levels_alldirs(FieldType::current_fp, finest_level),
        m_fields.get_mr_levels_alldirs(FieldType::current_cp, finest_level, skip_lev0_coarse_patch),
        m_fields.get_mr_levels_alldirs(FieldType::current_buf, finest_level, skip_lev0_coarse_patch),
        coarse_lev);

    if (m_fields.has(FieldType::rho_fp, finest_level) &&
        m_fields.has(FieldType::rho_cp, finest_level) &&
        m_fields.has(FieldType::rho_buf, finest_level)) {
        AddRhoFromFineLevelandSumBoundary(
            m_fields.get_mr_levels(FieldType::rho_fp, finest_level),
            m_fields.get_mr_levels(FieldType::rho_cp, finest_level, skip_lev0_coarse_patch),
            m_fields.get_mr_levels(FieldType::rho_buf, finest_level, skip_lev0_coarse_patch),
            coarse_lev, ncomps, ncomps);
    }

    EvolveE(fine_lev, PatchType::coarse, dt[fine_lev], cur_time + 0.5_rt * dt[fine_lev]);
    FillBoundaryE(fine_lev, PatchType::coarse, guard_cells.ng_FieldSolver,
                  WarpX::sync_nodal_points);

    EvolveB(fine_lev, PatchType::coarse, dt[fine_lev], SubcyclingHalf::SecondHalf, cur_time + 0.5_rt * dt[fine_lev]);
    EvolveF(fine_lev, PatchType::coarse, dt[fine_lev], /*rho_comp=*/1);

    if (do_pml) {
        FillBoundaryF(fine_lev, PatchType::fine, guard_cells.ng_FieldSolverF);
        DampPML(fine_lev, PatchType::coarse); // do it twice
        DampPML(fine_lev, PatchType::coarse);
        FillBoundaryE(fine_lev, PatchType::coarse, guard_cells.ng_alloc_EB);
    }

    FillBoundaryB(fine_lev, PatchType::coarse, guard_cells.ng_FieldSolver,
                  WarpX::sync_nodal_points);
    FillBoundaryF(fine_lev, PatchType::coarse, guard_cells.ng_FieldSolverF,
                  WarpX::sync_nodal_points);

    EvolveE(coarse_lev, PatchType::fine, 0.5_rt*dt[coarse_lev], cur_time + 0.5_rt*dt[coarse_lev]);
    FillBoundaryE(coarse_lev, PatchType::fine, guard_cells.ng_FieldSolver,
                  WarpX::sync_nodal_points);

    EvolveB(coarse_lev, PatchType::fine, 0.5_rt*dt[coarse_lev], SubcyclingHalf::SecondHalf, cur_time + 0.5_rt*dt[coarse_lev]);
    EvolveF(coarse_lev, PatchType::fine, 0.5_rt*dt[coarse_lev], /*rho_comp=*/1);

    if (do_pml) {
        if (moving_window_active(istep[0]+1)){
            // Exchange guard cells of PMLs only (0 cells are exchanged for the
            // regular B field MultiFab). This is required as B and F have just been
            // evolved.
            FillBoundaryB(coarse_lev, PatchType::fine, IntVect::TheZeroVector(),
                          WarpX::sync_nodal_points);
            FillBoundaryF(coarse_lev, PatchType::fine, IntVect::TheZeroVector(),
                          WarpX::sync_nodal_points);
        }
        DampPML(coarse_lev, PatchType::fine);
        if ( m_safe_guard_cells ) {
            FillBoundaryE(coarse_lev, PatchType::fine, guard_cells.ng_FieldSolver,
                          WarpX::sync_nodal_points);
        }
    }
    if ( m_safe_guard_cells ) {
        FillBoundaryB(coarse_lev, PatchType::fine, guard_cells.ng_FieldSolver,
                      WarpX::sync_nodal_points);
    }
}

void
WarpX::doFieldIonization ()
{
    using ablastr::fields::Direction;
    using warpx::fields::FieldType;

    for (int lev = 0; lev <= finest_level; ++lev) {
        mypc->doFieldIonization(
            lev,
            *m_fields.get(FieldType::Efield_aux, Direction{0}, lev),
            *m_fields.get(FieldType::Efield_aux, Direction{1}, lev),
            *m_fields.get(FieldType::Efield_aux, Direction{2}, lev),
            *m_fields.get(FieldType::Bfield_aux, Direction{0}, lev),
            *m_fields.get(FieldType::Bfield_aux, Direction{1}, lev),
            *m_fields.get(FieldType::Bfield_aux, Direction{2}, lev)
        );
    }
}

#ifdef WARPX_QED
void
WarpX::doQEDEvents ()
{
    using ablastr::fields::Direction;
    using warpx::fields::FieldType;

    for (int lev = 0; lev <= finest_level; ++lev) {
        mypc->doQedEvents(
            lev,
            *m_fields.get(FieldType::Efield_aux, Direction{0}, lev),
            *m_fields.get(FieldType::Efield_aux, Direction{1}, lev),
            *m_fields.get(FieldType::Efield_aux, Direction{2}, lev),
            *m_fields.get(FieldType::Bfield_aux, Direction{0}, lev),
            *m_fields.get(FieldType::Bfield_aux, Direction{1}, lev),
            *m_fields.get(FieldType::Bfield_aux, Direction{2}, lev)
        );
    }
}
#endif

void
WarpX::PushParticlesandDeposit (
    amrex::Real cur_time,
    bool skip_deposition,
    PositionPushType position_push_type,
    MomentumPushType momentum_push_type,
    ImplicitOptions const * implicit_options
)
{
    // Evolve particles to p^{n+1/2} and x^{n+1}
    // Deposit current, j^{n+1/2}
    for (int lev = 0; lev <= finest_level; ++lev) {
        PushParticlesandDeposit(
            lev,
            cur_time,
            SubcyclingHalf::None,
            skip_deposition,
            position_push_type,
            momentum_push_type,
            implicit_options
        );
    }
}

void
WarpX::PushParticlesandDeposit (
    int lev,
    amrex::Real cur_time,
    SubcyclingHalf subcycling_half,
    bool skip_deposition,
    PositionPushType position_push_type,
    MomentumPushType momentum_push_type,
    ImplicitOptions const * implicit_options
)
{
    using ablastr::fields::Direction;
    using warpx::fields::FieldType;

    std::string current_fp_string;

    if (WarpX::do_current_centering)
    {
        current_fp_string = "current_fp_nodal";
    }
    else if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay)
    {
        current_fp_string = "current_fp_vay";
    }
    else
    {
        current_fp_string = "current_fp";
    }

    mypc->Evolve(
        m_fields,
        lev,
        current_fp_string,
        cur_time,
        dt[lev],
        subcycling_half,
        skip_deposition,
        position_push_type,
        momentum_push_type,
        implicit_options
    );

    if (!skip_deposition) {
#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
        // This is called after all particles have deposited their current and charge.
        if (!implicit_options) {
            // Skip scaling J here for the implicit solvers: the total current is
            // accumulated from multiple containers after this call (see CumulateJ()
            // and ComputeJfromMassMatrices()), and is scaled in PreRHSOp().
            ApplyInverseVolumeScalingToCurrentDensity(
                m_fields.get(FieldType::current_fp, Direction{0}, lev),
                m_fields.get(FieldType::current_fp, Direction{1}, lev),
                m_fields.get(FieldType::current_fp, Direction{2}, lev),
                lev);
            if (m_fields.has_vector(FieldType::current_buf, lev)) {
                ApplyInverseVolumeScalingToCurrentDensity(
                    m_fields.get(FieldType::current_buf, Direction{0}, lev),
                    m_fields.get(FieldType::current_buf, Direction{1}, lev),
                    m_fields.get(FieldType::current_buf, Direction{2}, lev),
                    lev-1);
            }
        }
        // Unlike J, the charge density has no post-deposition accumulation step:
        // rho is reset and fully deposited within this call on both the explicit
        // and implicit paths, so it is scaled here in all cases.
        if (m_fields.has(FieldType::rho_fp, lev)) {
            ApplyInverseVolumeScalingToChargeDensity(m_fields.get(FieldType::rho_fp, lev), lev);
            if (m_fields.has(FieldType::rho_buf, lev)) {
                ApplyInverseVolumeScalingToChargeDensity(m_fields.get(FieldType::rho_buf, lev), lev-1);
            }
        }
// #else
        // I left this comment here as a reminder that currently the
        // boundary handling for cartesian grids are not matching the RZ handling
        // (done in the ApplyInverseScalingToChargeDensity function). The
        // Cartesian grid code had to be moved from here to after the application
        // of the filter to avoid incorrect results (moved to `SyncCurrentAndRho()`).
        // Might this be related to issue #1943?
#endif
        if (do_fluid_species && !implicit_options) {
            myfl->Evolve(m_fields,
                         lev,
                         current_fp_string,
                         cur_time,
                         skip_deposition
            );
        }
    }
}

/* \brief Apply perfect mirror condition inside the box (not at a boundary).
 * In practice, set all fields to 0 on a section of the simulation domain
 * (as for a perfect conductor with a given thickness).
 * The mirror normal direction has to be parallel to the z axis.
 */
void
WarpX::applyMirrors (Real time)
{
    using ablastr::fields::Direction;

    // something to do?
    if (m_num_mirrors == 0) {
        return;
    }

    // Loop over the mirrors
    for(int i_mirror=0; i_mirror<m_num_mirrors; ++i_mirror)
    {
        // Get mirror properties (lower and upper z bounds)
        amrex::Real z_min = m_mirror_z[i_mirror];
        amrex::Real z_max_tmp = z_min + m_mirror_z_width[i_mirror];

        // Boost quantities for boosted frame simulations
        if (gamma_boost>1)
        {
            z_min = z_min/gamma_boost - PhysConst::c*beta_boost*time;
            z_max_tmp = z_max_tmp/gamma_boost - PhysConst::c*beta_boost*time;
        }

        // Loop over levels
        for(int lev=0; lev<=finest_level; lev++)
        {
            // Mirror must contain at least m_mirror_z_npoints[i_mirror] cells
            const amrex::Real dz = WarpX::CellSize(lev)[2];
            const amrex::Real z_max = std::max(z_max_tmp, z_min+m_mirror_z_npoints[i_mirror]*dz);

            // Set each field on the fine patch to zero between z_min and z_max
            NullifyMF(m_fields, "Efield_fp", Direction{0}, lev, z_min, z_max);
            NullifyMF(m_fields, "Efield_fp", Direction{1}, lev, z_min, z_max);
            NullifyMF(m_fields, "Efield_fp", Direction{2}, lev, z_min, z_max);
            NullifyMF(m_fields, "Bfield_fp", Direction{0}, lev, z_min, z_max);
            NullifyMF(m_fields, "Bfield_fp", Direction{1}, lev, z_min, z_max);
            NullifyMF(m_fields, "Bfield_fp", Direction{2}, lev, z_min, z_max);

            // If div(E)/div(B) cleaning are used, set F/G field to zero
            NullifyMF(m_fields, "F_fp", lev, z_min, z_max);
            NullifyMF(m_fields, "G_fp", lev, z_min, z_max);

            if (lev>0)
            {
                // Set each field on the coarse patch to zero between z_min and z_max
                NullifyMF(m_fields, "Efield_cp", Direction{0}, lev, z_min, z_max);
                NullifyMF(m_fields, "Efield_cp", Direction{1}, lev, z_min, z_max);
                NullifyMF(m_fields, "Efield_cp", Direction{2}, lev, z_min, z_max);
                NullifyMF(m_fields, "Bfield_cp", Direction{0}, lev, z_min, z_max);
                NullifyMF(m_fields, "Bfield_cp", Direction{1}, lev, z_min, z_max);
                NullifyMF(m_fields, "Bfield_cp", Direction{2}, lev, z_min, z_max);

                // If div(E)/div(B) cleaning are used, set F/G field to zero
                NullifyMF(m_fields, "F_cp", lev, z_min, z_max);
                NullifyMF(m_fields, "G_cp", lev, z_min, z_max);
            }
        }
    }
}

void
WarpX::HandleSignals()
{
    SignalHandling::WaitSignals();

    // SIGNAL_REQUESTS_BREAK is handled directly in WarpX::Evolve

    if (SignalHandling::TestAndResetActionRequestFlag(SignalHandling::SIGNAL_REQUESTS_CHECKPOINT)) {
        multi_diags->FilterComputePackFlushLastTimestep( istep[0] );
        ExecutePythonCallback("oncheckpointsignal");
    }
}