PSATD Runtime Control (#1300)

* Docs: PSATD Runtime Option

* Tests: PSATD Runtime Option

Add new runtime option to PSATD regression test matrix.

* PICMI: PSATD runtime option

* Source: PSATD Runtime Option

39 files changed, 505 insertions, 396 deletions

diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst
index d0d490875..d2b18c628 100644
--- a/Docs/source/running_cpp/parameters.rst
+++ b/Docs/source/running_cpp/parameters.rst
@@ -1078,9 +1078,8 @@ Numerics and algorithms
 * ``algo.current_deposition`` (`string`, optional)
     This parameter selects the algorithm for the deposition of the current density.
     Available options are: ``direct``, ``esirkepov``, and ``vay``. The default choice
-    is ``esirkepov`` if WarpX is compiled with the FDTD solver (that is, with
-    ``USE_PSATD=FALSE``) and ``direct`` if WarpX is compiled with the standard or
-    Galilean PSATD solver (that is, with ``USE_PSATD=TRUE``).
+    is ``esirkepov`` for FDTD maxwell solvers and ``direct`` for standard or
+    Galilean PSATD solver (that is, with ``algo.maxwell_solver = psatd``).
 
     1. ``direct``
 
@@ -1136,9 +1135,9 @@ Numerics and algorithms
      - ``yee``: Yee FDTD solver.
      - ``ckc``: (not available in ``RZ`` geometry) Cole-Karkkainen solver with Cowan
        coefficients (see `Cowan, PRSTAB 16 (2013) <https://journals.aps.org/prab/abstract/10.1103/PhysRevSTAB.16.041303>`__)
+     - ``psatd``: Pseudo-spectral solver (see :ref:`theory <theory-pic-mwsolve-psatd>`)
 
      If ``algo.maxwell_solver`` is not specified, ``yee`` is the default.
-     Note: this option is currently ignored with PSATD.
 
 * ``algo.em_solver_medium`` (`string`, optional)
     The medium for evaluating the Maxwell solver. Available options are :
diff --git a/Docs/source/theory/picsar_theory.rst b/Docs/source/theory/picsar_theory.rst
index 0519778f3..b4b18e1a4 100644
--- a/Docs/source/theory/picsar_theory.rst
+++ b/Docs/source/theory/picsar_theory.rst
@@ -10,6 +10,8 @@
 
    \linenumbers
 
+.. _theory-pic:
+
 The electromagnetic Particle-In-Cell method
 ===========================================
 
@@ -49,6 +51,8 @@ on the grid from the particles’ positions and velocities, while the
 electric and magnetic field components are interpolated from the grid
 to the particles’ positions for the velocity update.
 
+.. _theory-pic-push:
+
 Particle push
 -------------
 
@@ -64,6 +68,8 @@ equations of motion is given by
 In order to close the system, :math:`\bar{\mathbf{v}}^{i}` must be
 expressed as a function of the other quantities. The two implementations that have become the most popular are presented below.
 
+.. _theory-pic-push-boris:
+
 Boris relativistic velocity rotation
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -94,6 +100,8 @@ where :math:`\mathbf{t}=\left(q\Delta t/2m\right)\mathbf{B}^{i}/\bar{\gamma}^{i}
 
 The Boris implementation is second-order accurate, time-reversible and fast. Its implementation is very widespread and used in the vast majority of PIC codes.
 
+.. _theory-pic-push-vay:
+
 Vay Lorentz-invariant formulation
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -125,6 +133,8 @@ charged particle beams, where the accurate account of the cancellation
 of the self-generated electric and magnetic fields is essential, as
 shown in (Vay 2008).
 
+.. _theory-pic-mwsolve:
+
 Field solve
 -----------
 
@@ -139,6 +149,8 @@ analytical time-domain (PSATD) and pseudo-spectral time-domain (PSTD)
 algorithms. Extension to multiresolution (or mesh refinement) PIC
 is described in, e.g. (Vay et al. 2012; Vay, Adam, and Heron 2004).
 
+.. _theory-pic-mwsolve-fdtd:
+
 Finite-Difference Time-Domain (FDTD)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -179,6 +191,8 @@ center of the cell faces. Knowing the current densities at half-integer steps,
 the electric field components are updated alternately with the magnetic
 field components at integer and half-integer steps respectively.
 
+.. _theory-pic-mwsolve-nsfdtd:
+
 Non-Standard Finite-Difference Time-Domain (NSFDTD)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -236,6 +250,8 @@ the same paper that removing the Nyquist component in all the source
 terms using a bilinear filter (see description of the filter below)
 suppresses this instability.
 
+.. _theory-pic-mwsolve-psatd:
+
 Pseudo Spectral Analytical Time Domain (PSATD)
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
diff --git a/Examples/Tests/Maxwell_Hybrid_QED/inputs_2d b/Examples/Tests/Maxwell_Hybrid_QED/inputs_2d
index ba7d9c33e..9905e0fa5 100644
--- a/Examples/Tests/Maxwell_Hybrid_QED/inputs_2d
+++ b/Examples/Tests/Maxwell_Hybrid_QED/inputs_2d
@@ -16,6 +16,7 @@ warpx.quantum_xi = 1.e-23
 #################################
 ############ NUMERICS ###########
 #################################
+algo.maxwell_solver = psatd
 warpx.verbose = 0
 warpx.use_filter = 1
 warpx.cfl = 1.
diff --git a/Examples/Tests/averaged_galilean/inputs_avg_2d b/Examples/Tests/averaged_galilean/inputs_avg_2d
index f553009d3..a9f7cedbb 100644
--- a/Examples/Tests/averaged_galilean/inputs_avg_2d
+++ b/Examples/Tests/averaged_galilean/inputs_avg_2d
@@ -22,7 +22,7 @@ geometry.prob_hi     = 12.3776 49.5104
 #################################
 warpx.verbose = 1
 
-algo.maxwell_solver = ckc
+algo.maxwell_solver = psatd
 algo.current_deposition = direct
 algo.particle_pusher = vay
 
diff --git a/Examples/Tests/averaged_galilean/inputs_avg_3d b/Examples/Tests/averaged_galilean/inputs_avg_3d
index 747d0f46a..0579de4ac 100644
--- a/Examples/Tests/averaged_galilean/inputs_avg_3d
+++ b/Examples/Tests/averaged_galilean/inputs_avg_3d
@@ -20,7 +20,7 @@ geometry.prob_hi     = 9.67 9.67 19.34
 #################################
 warpx.verbose = 1
 
-algo.maxwell_solver = ckc
+algo.maxwell_solver = psatd
 algo.current_deposition = direct
 algo.particle_pusher = vay
 
diff --git a/Examples/Tests/galilean/inputs_2d b/Examples/Tests/galilean/inputs_2d
index c52471983..9b7330b18 100644
--- a/Examples/Tests/galilean/inputs_2d
+++ b/Examples/Tests/galilean/inputs_2d
@@ -20,6 +20,7 @@ geometry.prob_hi     = 38.68 38.68
 warpx.verbose = 1
 
 algo.current_deposition = direct
+algo.maxwell_solver = psatd
 algo.particle_pusher = vay
 
 warpx.cfl = 1.
diff --git a/Examples/Tests/galilean/inputs_3d b/Examples/Tests/galilean/inputs_3d
index 7b2905437..050ed0b6f 100644
--- a/Examples/Tests/galilean/inputs_3d
+++ b/Examples/Tests/galilean/inputs_3d
@@ -21,6 +21,7 @@ geometry.prob_hi     = 9.67 9.67 9.67
 warpx.verbose = 1
 
 algo.current_deposition = direct
+algo.maxwell_solver = psatd
 algo.particle_pusher = vay
 
 warpx.cfl = 1.
diff --git a/Examples/Tests/galilean/inputs_rz b/Examples/Tests/galilean/inputs_rz
index f2805faf0..1a207257f 100644
--- a/Examples/Tests/galilean/inputs_rz
+++ b/Examples/Tests/galilean/inputs_rz
@@ -20,6 +20,7 @@ geometry.prob_hi     = 38.68 38.68
 warpx.verbose = 1
 
 algo.current_deposition = direct
+algo.maxwell_solver = psatd
 algo.particle_pusher = vay
 
 warpx.cfl = 1.
diff --git a/Examples/Tests/spectral_staggered/inputs_hybrid_2d b/Examples/Tests/spectral_staggered/inputs_hybrid_2d
index 32884a88d..ee197ebc9 100644
--- a/Examples/Tests/spectral_staggered/inputs_hybrid_2d
+++ b/Examples/Tests/spectral_staggered/inputs_hybrid_2d
@@ -9,6 +9,7 @@ geometry.is_periodic = 0 0
 geometry.prob_lo = -90.e-6 -70.e-6
 geometry.prob_hi =  90.e-6   0.e-6
 
+algo.maxwell_solver = psatd
 algo.current_deposition = direct
 algo.charge_deposition = standard
 algo.field_gathering = momentum-conserving
diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py
index 784d9af01..05986e785 100644
--- a/Python/pywarpx/picmi.py
+++ b/Python/pywarpx/picmi.py
@@ -503,9 +503,8 @@ class ElectromagneticSolver(picmistandard.PICMI_ElectromagneticSolver):
                     self.galilean_velocity = [self.galilean_velocity[0], 0., self.galilean_velocity[1]]
                 pywarpx.psatd.v_galilean = np.array(self.galilean_velocity)/constants.c
 
-        else:
-            # --- Same method names are used, though mapped to lower case.
-            pywarpx.algo.maxwell_solver = self.method
+        # --- Same method names are used, though mapped to lower case.
+        pywarpx.algo.maxwell_solver = self.method
 
         if self.cfl is not None:
             pywarpx.warpx.cfl = self.cfl
diff --git a/Regression/WarpX-GPU-tests.ini b/Regression/WarpX-GPU-tests.ini
index 48172bbb5..6bef020f3 100644
--- a/Regression/WarpX-GPU-tests.ini
+++ b/Regression/WarpX-GPU-tests.ini
@@ -95,7 +95,7 @@ tolerance = 1e-9
 #[pml_x_psatd]
 #buildDir = .
 #inputFile = Examples/Tests/PML/inputs_2d
-#runtime_params = warpx.do_dynamic_scheduling=0
+#runtime_params = algo.maxwell_solver=psatd warpx.do_dynamic_scheduling=0
 #dim = 2
 #addToCompileString = USE_PSATD=TRUE USE_GPU=TRUE
 #restartTest = 0
@@ -377,6 +377,7 @@ tolerance = 5e-11
 [Langmuir_multi_psatd]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_3d_multi_rt
+runtime_params = algo.maxwell_solver=psatd psatd.fftw_plan_measure=0 warpx.cfl = 0.5773502691896258
 dim = 3
 addToCompileString = USE_PSATD=TRUE USE_GPU=TRUE
 restartTest = 0
@@ -387,7 +388,6 @@ numthreads = 1
 compileTest = 0
 doVis = 0
 compareParticles = 0
-runtime_params = psatd.fftw_plan_measure=0 warpx.cfl = 0.5773502691896258
 particleTypes = electrons positrons
 analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir_multi.py
 analysisOutputImage = langmuir_multi_analysis.png
@@ -396,6 +396,7 @@ tolerance = 2.e-9
 [Langmuir_multi_psatd_nodal]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_3d_multi_rt
+runtime_params = algo.maxwell_solver=psatd psatd.fftw_plan_measure=0 warpx.do_dynamic_scheduling=0 warpx.do_nodal=1 algo.current_deposition=direct warpx.cfl = 0.5773502691896258
 dim = 3
 addToCompileString = USE_PSATD=TRUE USE_GPU=TRUE
 restartTest = 0
@@ -406,7 +407,6 @@ numthreads = 1
 compileTest = 0
 doVis = 0
 compareParticles = 0
-runtime_params = psatd.fftw_plan_measure=0 warpx.do_dynamic_scheduling=0 warpx.do_nodal=1 algo.current_deposition=direct warpx.cfl = 0.5773502691896258
 particleTypes = electrons positrons
 analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir_multi.py
 analysisOutputImage = langmuir_multi_analysis.png
@@ -434,6 +434,7 @@ tolerance = 5e-11
 [Langmuir_multi_2d_psatd]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rt
+runtime_params = algo.maxwell_solver=psatd psatd.fftw_plan_measure=0 diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell warpx.cfl = 0.7071067811865475
 dim = 2
 addToCompileString = USE_PSATD=TRUE USE_GPU=TRUE
 restartTest = 0
@@ -444,7 +445,6 @@ numthreads = 1
 compileTest = 0
 doVis = 0
 compareParticles = 0
-runtime_params = psatd.fftw_plan_measure=0 diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell warpx.cfl = 0.7071067811865475
 particleTypes = electrons positrons
 analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir_multi_2d.py
 analysisOutputImage = langmuir_multi_2d_analysis.png
@@ -453,6 +453,7 @@ tolerance = 5e-11
 # [Langmuir_multi_2d_psatd_nodal]
 # buildDir = .
 # inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rt
+# runtime_params = algo.maxwell_solver=psatd psatd.fftw_plan_measure=0 warpx.do_nodal=1 algo.current_deposition=direct diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell
 # dim = 2
 # addToCompileString = USE_PSATD=TRUE USE_GPU=TRUE
 # restartTest = 0
@@ -463,7 +464,6 @@ tolerance = 5e-11
 # compileTest = 0
 # doVis = 0
 # compareParticles = 0
-# runtime_params =  psatd.fftw_plan_measure=0 warpx.do_nodal=1 algo.current_deposition=direct diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell
 # particleTypes = electrons positrons
 # analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir_multi_2d.py
 # analysisOutputImage = langmuir_multi_2d_analysis.png
diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini
index cdf000a9f..000cab93e 100644
--- a/Regression/WarpX-tests.ini
+++ b/Regression/WarpX-tests.ini
@@ -95,7 +95,7 @@ tolerance = 1.e-14
 [pml_x_psatd]
 buildDir = .
 inputFile = Examples/Tests/PML/inputs_2d
-runtime_params = psatd.update_with_rho=1 warpx.do_dynamic_scheduling=0 diag1.fields_to_plot = Ex Ey Ez Bx By Bz rho divE warpx.cfl = 0.7071067811865475
+runtime_params = algo.maxwell_solver=psatd psatd.update_with_rho=1 warpx.do_dynamic_scheduling=0 diag1.fields_to_plot = Ex Ey Ez Bx By Bz rho divE warpx.cfl = 0.7071067811865475
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -362,7 +362,7 @@ tolerance = 1.e-14
 [Langmuir_multi_psatd]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_3d_multi_rt
-runtime_params = psatd.fftw_plan_measure=0 warpx.cfl = 0.5773502691896258
+runtime_params = algo.maxwell_solver=psatd psatd.fftw_plan_measure=0 warpx.cfl = 0.5773502691896258
 dim = 3
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -381,7 +381,7 @@ tolerance = 5.e-11
 [Langmuir_multi_psatd_current_correction]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_3d_multi_rt
-runtime_params = algo.current_deposition=esirkepov psatd.fftw_plan_measure=0 psatd.periodic_single_box_fft=1 psatd.current_correction=1 diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz part_per_cell rho divE warpx.cfl = 0.5773502691896258
+runtime_params = algo.maxwell_solver=psatd algo.current_deposition=esirkepov psatd.fftw_plan_measure=0 psatd.periodic_single_box_fft=1 psatd.current_correction=1 diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz part_per_cell rho divE warpx.cfl = 0.5773502691896258
 dim = 3
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -400,7 +400,7 @@ analysisOutputImage = langmuir_multi_analysis.png
 [Langmuir_multi_psatd_current_correction_nodal]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_3d_multi_rt
-runtime_params = algo.current_deposition=direct psatd.fftw_plan_measure=0 psatd.periodic_single_box_fft=1 psatd.current_correction=1 warpx.do_nodal=1 diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz part_per_cell rho divE warpx.cfl = 0.5773502691896258
+runtime_params = algo.maxwell_solver=psatd algo.current_deposition=direct psatd.fftw_plan_measure=0 psatd.periodic_single_box_fft=1 psatd.current_correction=1 warpx.do_nodal=1 diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz part_per_cell rho divE warpx.cfl = 0.5773502691896258
 dim = 3
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -419,7 +419,7 @@ analysisOutputImage = langmuir_multi_analysis.png
 [Langmuir_multi_psatd_Vay_deposition]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_3d_multi_rt
-runtime_params = psatd.fftw_plan_measure=0 psatd.periodic_single_box_fft=1 algo.current_deposition=vay diag1.fields_to_plot = Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.5773502691896258
+runtime_params = algo.maxwell_solver=psatd psatd.fftw_plan_measure=0 psatd.periodic_single_box_fft=1 algo.current_deposition=vay diag1.fields_to_plot = Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.5773502691896258
 dim = 3
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -438,7 +438,7 @@ analysisOutputImage = langmuir_multi_analysis.png
 [Langmuir_multi_psatd_Vay_deposition_nodal]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_3d_multi_rt
-runtime_params = warpx.do_nodal=1 psatd.fftw_plan_measure=0 psatd.periodic_single_box_fft=1 algo.current_deposition=vay diag1.fields_to_plot = Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.5773502691896258
+runtime_params = algo.maxwell_solver=psatd warpx.do_nodal=1 psatd.fftw_plan_measure=0 psatd.periodic_single_box_fft=1 algo.current_deposition=vay diag1.fields_to_plot = Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.5773502691896258
 dim = 3
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -457,7 +457,7 @@ analysisOutputImage = langmuir_multi_analysis.png
 [Langmuir_multi_psatd_momentum_conserving]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_3d_multi_rt
-runtime_params = psatd.fftw_plan_measure=0 algo.field_gathering=momentum-conserving warpx.cfl = 0.5773502691896258
+runtime_params = algo.maxwell_solver=psatd psatd.fftw_plan_measure=0 algo.field_gathering=momentum-conserving warpx.cfl = 0.5773502691896258
 dim = 3
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -476,7 +476,7 @@ tolerance = 5.e-11
 [Langmuir_multi_psatd_nodal]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_3d_multi_rt
-runtime_params = psatd.fftw_plan_measure=0 warpx.do_dynamic_scheduling=0 warpx.do_nodal=1 algo.current_deposition=direct warpx.cfl = 0.5773502691896258
+runtime_params = algo.maxwell_solver=psatd psatd.fftw_plan_measure=0 warpx.do_dynamic_scheduling=0 warpx.do_nodal=1 algo.current_deposition=direct warpx.cfl = 0.5773502691896258
 dim = 3
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -495,7 +495,7 @@ tolerance = 5.e-11
 [Langmuir_multi_psatd_single_precision]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_3d_multi_rt
-runtime_params = psatd.fftw_plan_measure=0  warpx.cfl = 0.5773502691896258
+runtime_params = algo.maxwell_solver=psatd psatd.fftw_plan_measure=0  warpx.cfl = 0.5773502691896258
 dim = 3
 addToCompileString = USE_PSATD=TRUE PRECISION=FLOAT USE_SINGLE_PRECISION_PARTICLES=TRUE
 restartTest = 0
@@ -533,7 +533,7 @@ tolerance = 1.e-14
 [Langmuir_multi_2d_psatd]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rt
-runtime_params = psatd.fftw_plan_measure=0 diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell warpx.cfl = 0.7071067811865475
+runtime_params = algo.maxwell_solver=psatd psatd.fftw_plan_measure=0 diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell warpx.cfl = 0.7071067811865475
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -552,7 +552,7 @@ tolerance = 1.e-14
 [Langmuir_multi_2d_psatd_momentum_conserving]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rt
-runtime_params = algo.field_gathering=momentum-conserving psatd.fftw_plan_measure=0 diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell warpx.cfl = 0.7071067811865475
+runtime_params = algo.maxwell_solver=psatd algo.field_gathering=momentum-conserving psatd.fftw_plan_measure=0 diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell warpx.cfl = 0.7071067811865475
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -571,7 +571,7 @@ tolerance = 1.e-14
 [Langmuir_multi_2d_psatd_current_correction]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rt
-runtime_params = amr.max_grid_size=128 algo.current_deposition=esirkepov psatd.fftw_plan_measure=0 psatd.periodic_single_box_fft=1 psatd.current_correction=1 diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot =Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.7071067811865475
+runtime_params = algo.maxwell_solver=psatd amr.max_grid_size=128 algo.current_deposition=esirkepov psatd.fftw_plan_measure=0 psatd.periodic_single_box_fft=1 psatd.current_correction=1 diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot =Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.7071067811865475
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -589,7 +589,7 @@ analysisOutputImage = langmuir_multi_2d_analysis.png
 [Langmuir_multi_2d_psatd_current_correction_nodal]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rt
-runtime_params = amr.max_grid_size=128 algo.current_deposition=direct psatd.fftw_plan_measure=0 psatd.periodic_single_box_fft=1 psatd.current_correction=1 warpx.do_nodal=1 diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot =Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.7071067811865475
+runtime_params = algo.maxwell_solver=psatd amr.max_grid_size=128 algo.current_deposition=direct psatd.fftw_plan_measure=0 psatd.periodic_single_box_fft=1 psatd.current_correction=1 warpx.do_nodal=1 diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot =Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.7071067811865475
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -607,7 +607,7 @@ analysisOutputImage = langmuir_multi_2d_analysis.png
 [Langmuir_multi_2d_psatd_Vay_deposition]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rt
-runtime_params = amr.max_grid_size=128 psatd.fftw_plan_measure=0 psatd.periodic_single_box_fft=1 algo.current_deposition=vay diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot = Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.7071067811865475
+runtime_params = algo.maxwell_solver=psatd amr.max_grid_size=128 psatd.fftw_plan_measure=0 psatd.periodic_single_box_fft=1 algo.current_deposition=vay diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot = Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.7071067811865475
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -625,7 +625,7 @@ analysisOutputImage = langmuir_multi_2d_analysis.png
 [Langmuir_multi_2d_psatd_Vay_deposition_nodal]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rt
-runtime_params = amr.max_grid_size=128 warpx.do_nodal=1 psatd.fftw_plan_measure=0 psatd.periodic_single_box_fft=1 algo.current_deposition=vay diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot = Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.7071067811865475
+runtime_params = algo.maxwell_solver=psatd amr.max_grid_size=128 warpx.do_nodal=1 psatd.fftw_plan_measure=0 psatd.periodic_single_box_fft=1 algo.current_deposition=vay diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot = Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.7071067811865475
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -643,7 +643,7 @@ analysisOutputImage = langmuir_multi_2d_analysis.png
 [Langmuir_multi_2d_psatd_nodal]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rt
-runtime_params = psatd.fftw_plan_measure=0 warpx.do_nodal=1 algo.current_deposition=direct diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell warpx.cfl = 0.7071067811865475
+runtime_params = algo.maxwell_solver=psatd psatd.fftw_plan_measure=0 warpx.do_nodal=1 algo.current_deposition=direct diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell warpx.cfl = 0.7071067811865475
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -681,7 +681,7 @@ tolerance = 1.e-14
 [Langmuir_multi_rz_psatd]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rz_rt
-runtime_params = diag1.electrons.variables=w ux uy uz diag1.ions.variables=w ux uy uz diag1.dump_rz_modes=0 algo.current_deposition=direct warpx.do_dive_cleaning=0
+runtime_params = algo.maxwell_solver=psatd diag1.electrons.variables=w ux uy uz diag1.ions.variables=w ux uy uz diag1.dump_rz_modes=0 algo.current_deposition=direct warpx.do_dive_cleaning=0
 dim = 2
 addToCompileString = USE_RZ=TRUE USE_PSATD=TRUE BLAS_LIB=-lblas LAPACK_LIB=-llapack
 restartTest = 0
@@ -700,7 +700,7 @@ tolerance = 1.e-14
 [Langmuir_multi_rz_psatd_current_correction]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rz_rt
-runtime_params = diag1.electrons.variables=w ux uy uz diag1.ions.variables=w ux uy uz diag1.dump_rz_modes=0 algo.current_deposition=direct warpx.do_dive_cleaning=0 amr.max_grid_size=128 psatd.periodic_single_box_fft=1 psatd.current_correction=1 diag1.fields_to_plot=jx jz Ex Ez By rho divE
+runtime_params = algo.maxwell_solver=psatd diag1.electrons.variables=w ux uy uz diag1.ions.variables=w ux uy uz diag1.dump_rz_modes=0 algo.current_deposition=direct warpx.do_dive_cleaning=0 amr.max_grid_size=128 psatd.periodic_single_box_fft=1 psatd.current_correction=1 diag1.fields_to_plot=jx jz Ex Ez By rho divE
 dim = 2
 addToCompileString = USE_RZ=TRUE USE_PSATD=TRUE BLAS_LIB=-lblas LAPACK_LIB=-llapack
 restartTest = 0
@@ -1515,7 +1515,7 @@ tolerance = 1.e-14
 [Maxwell_Hybrid_QED_solver]
 buildDir = .
 inputFile = Examples/Tests/Maxwell_Hybrid_QED/inputs_2d
-runtime_params = warpx.cfl = 0.7071067811865475
+runtime_params = warpx.cfl=0.7071067811865475
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -1582,6 +1582,7 @@ tolerance = 1e-12
 [galilean_2d_psatd]
 buildDir = .
 inputFile = Examples/Tests/galilean/inputs_2d
+runtime_params = warpx.do_nodal=1 algo.current_deposition=direct
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -1592,7 +1593,6 @@ numthreads = 1
 compileTest = 0
 doVis = 0
 compareParticles = 1
-runtime_params = warpx.do_nodal=1 algo.current_deposition=direct
 particleTypes = electrons ions
 analysisRoutine = Examples/Tests/galilean/analysis_2d.py
 tolerance = 1.e-14
@@ -1690,6 +1690,7 @@ tolerance = 1.e-14
 [galilean_3d_psatd]
 buildDir = .
 inputFile = Examples/Tests/galilean/inputs_3d
+runtime_params = warpx.do_nodal=1 algo.current_deposition=direct
 dim = 3
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -1700,7 +1701,6 @@ numthreads = 1
 compileTest = 0
 doVis = 0
 compareParticles = 1
-runtime_params = warpx.do_nodal=1 algo.current_deposition=direct
 particleTypes = electrons ions
 analysisRoutine = Examples/Tests/galilean/analysis_3d.py
 tolerance = 1.e-14
@@ -1726,6 +1726,7 @@ tolerance = 1.e-14
 [averaged_galilean_2d_psatd]
 buildDir = .
 inputFile = Examples/Tests/averaged_galilean/inputs_avg_2d
+runtime_params =
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -1743,6 +1744,7 @@ tolerance = 1e-6
 [averaged_galilean_3d_psatd]
 buildDir = .
 inputFile = Examples/Tests/averaged_galilean/inputs_avg_3d
+runtime_params =
 dim = 3
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H
index e01d89aff..3130f84b2 100644
--- a/Source/BoundaryConditions/PML.H
+++ b/Source/BoundaryConditions/PML.H
@@ -101,9 +101,8 @@ public:
     PML (const amrex::BoxArray& ba, const amrex::DistributionMapping& dm,
          const amrex::Geometry* geom, const amrex::Geometry* cgeom,
          int ncell, int delta, int ref_ratio,
-#ifdef WARPX_USE_PSATD
+         // PSATD
          amrex::Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal,
-#endif
          int do_dive_cleaning, int do_moving_window,
          int pml_has_particles, int do_pml_in_domain,
          const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(),
diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp
index 7bf11a227..2ede4dbf2 100644
--- a/Source/BoundaryConditions/PML.cpp
+++ b/Source/BoundaryConditions/PML.cpp
@@ -427,9 +427,8 @@ MultiSigmaBox::ComputePMLFactorsE (const Real* dx, Real dt)
 PML::PML (const BoxArray& grid_ba, const DistributionMapping& /*grid_dm*/,
           const Geometry* geom, const Geometry* cgeom,
           int ncell, int delta, int ref_ratio,
-#ifdef WARPX_USE_PSATD
+          // PSATD
           Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal,
-#endif
           int do_dive_cleaning, int do_moving_window,
           int /*pml_has_particles*/, int do_pml_in_domain,
           const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi)
@@ -478,24 +477,25 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& /*grid_dm*/,
     int ngf_int = (do_moving_window) ? 2 : 0;
     if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::CKC) ngf_int = std::max( ngf_int, 1 );
     IntVect ngf = IntVect(AMREX_D_DECL(ngf_int, ngf_int, ngf_int));
-#ifdef WARPX_USE_PSATD
-    // Increase the number of guard cells, in order to fit the extent
-    // of the stencil for the spectral solver
-    IntVect ngFFT;
-    if (do_nodal) {
-        ngFFT = IntVect(AMREX_D_DECL(nox_fft, noy_fft, noz_fft));
-    } else {
-        ngFFT = IntVect(AMREX_D_DECL(nox_fft/2, noy_fft/2, noz_fft/2));
+
+    if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+        // Increase the number of guard cells, in order to fit the extent
+        // of the stencil for the spectral solver
+        IntVect ngFFT;
+        if (do_nodal) {
+            ngFFT = IntVect(AMREX_D_DECL(nox_fft, noy_fft, noz_fft));
+        } else {
+            ngFFT = IntVect(AMREX_D_DECL(nox_fft / 2, noy_fft / 2, noz_fft / 2));
+        }
+        // Set the number of guard cells to the maximum of each field
+        // (all fields should have the same number of guard cells)
+        ngFFT = ngFFT.max(nge);
+        ngFFT = ngFFT.max(ngb);
+        ngFFT = ngFFT.max(ngf);
+        nge = ngFFT;
+        ngb = ngFFT;
+        ngf = ngFFT;
     }
-    // Set the number of guard cells to the maximum of each field
-    // (all fields should have the same number of guard cells)
-    ngFFT = ngFFT.max(nge);
-    ngFFT = ngFFT.max(ngb);
-    ngFFT = ngFFT.max(ngf);
-    nge = ngFFT;
-    ngb = ngFFT;
-    ngf = ngFFT;
-#endif
 
     pml_E_fp[0] = std::make_unique<MultiFab>(amrex::convert( ba,
         WarpX::GetInstance().getEfield_fp(0,0).ixType().toIntVect() ), dm, 3, nge );
@@ -541,27 +541,29 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& /*grid_dm*/,
         sigba_fp = std::make_unique<MultiSigmaBox>(ba, dm, grid_ba, geom->CellSize(), ncell, delta);
     }
 
-
-#ifdef WARPX_USE_PSATD
-    AMREX_ALWAYS_ASSERT_WITH_MESSAGE( do_pml_in_domain==false,
-        "PSATD solver cannot be used with `do_pml_in_domain`.");
-    const bool in_pml = true; // Tells spectral solver to use split-PML equations
-    const RealVect dx{AMREX_D_DECL(geom->CellSize(0), geom->CellSize(1), geom->CellSize(2))};
-    // Get the cell-centered box, with guard cells
-    BoxArray realspace_ba = ba;  // Copy box
-    Array<Real,3> v_galilean_zero = {0., 0., 0.};
-    Array<Real,3> v_comoving_zero = {0., 0., 0.};
-    realspace_ba.enclosedCells().grow(nge); // cell-centered + guard cells
-    spectral_solver_fp = std::make_unique<SpectralSolver>(realspace_ba, dm,
-        nox_fft, noy_fft, noz_fft, do_nodal, v_galilean_zero, v_comoving_zero, dx, dt, in_pml );
+    if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+#ifndef WARPX_USE_PSATD
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false,
+                                         "PML: PSATD solver selected but not built.");
+#else
+        const bool in_pml = true; // Tells spectral solver to use split-PML equations
+        const RealVect dx{AMREX_D_DECL(geom->CellSize(0), geom->CellSize(1), geom->CellSize(2))};
+        // Get the cell-centered box, with guard cells
+        BoxArray realspace_ba = ba;  // Copy box
+        Array<Real,3> const v_galilean_zero = {0., 0., 0.};
+        Array<Real,3> const v_comoving_zero = {0., 0., 0.};
+        realspace_ba.enclosedCells().grow(nge); // cell-centered + guard cells
+        spectral_solver_fp = std::make_unique<SpectralSolver>(realspace_ba, dm,
+            nox_fft, noy_fft, noz_fft, do_nodal, v_galilean_zero, v_comoving_zero, dx, dt, in_pml );
 #endif
+    }
 
     if (cgeom)
     {
-#ifndef WARPX_USE_PSATD
-        nge = IntVect(AMREX_D_DECL(1, 1, 1));
-        ngb = IntVect(AMREX_D_DECL(1, 1, 1));
-#endif
+        if (WarpX::maxwell_solver_id != MaxwellSolverAlgo::PSATD) {
+            nge = IntVect(AMREX_D_DECL(1, 1, 1));
+            ngb = IntVect(AMREX_D_DECL(1, 1, 1));
+        }
 
         BoxArray grid_cba = grid_ba;
         grid_cba.coarsen(ref_ratio);
@@ -615,16 +617,23 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& /*grid_dm*/,
             sigba_cp = std::make_unique<MultiSigmaBox>(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta);
         }
 
-#ifdef WARPX_USE_PSATD
-        const RealVect cdx{AMREX_D_DECL(cgeom->CellSize(0), cgeom->CellSize(1), cgeom->CellSize(2))};
-        // Get the cell-centered box, with guard cells
-        BoxArray realspace_cba = cba;  // Copy box
-        // const bool in_pml = true; // Tells spectral solver to use split-PML equations
-
-        realspace_cba.enclosedCells().grow(nge); // cell-centered + guard cells
-        spectral_solver_cp = std::make_unique<SpectralSolver>(realspace_cba, cdm,
-            nox_fft, noy_fft, noz_fft, do_nodal, v_galilean_zero, v_comoving_zero, cdx, dt, in_pml );
+        if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+#ifndef WARPX_USE_PSATD
+            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false,
+                                             "PML: PSATD solver selected but not built.");
+#else
+            const RealVect cdx{AMREX_D_DECL(cgeom->CellSize(0), cgeom->CellSize(1), cgeom->CellSize(2))};
+            // Get the cell-centered box, with guard cells
+            BoxArray realspace_cba = cba;  // Copy box
+            Array<Real,3> const v_galilean_zero = {0., 0., 0.};
+            Array<Real,3> const v_comoving_zero = {0., 0., 0.};
+            const bool in_pml = true; // Tells spectral solver to use split-PML equations
+
+            realspace_cba.enclosedCells().grow(nge); // cell-centered + guard cells
+            spectral_solver_cp = std::make_unique<SpectralSolver>(realspace_cba, cdm,
+                nox_fft, noy_fft, noz_fft, do_nodal, v_galilean_zero, v_comoving_zero, cdx, dt, in_pml );
 #endif
+        }
     }
 }
 
diff --git a/Source/Diagnostics/ComputeDiagFunctors/DivEFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/DivEFunctor.cpp
index 26104bdf3..6abdb8938 100644
--- a/Source/Diagnostics/ComputeDiagFunctors/DivEFunctor.cpp
+++ b/Source/Diagnostics/ComputeDiagFunctors/DivEFunctor.cpp
@@ -20,13 +20,13 @@ DivEFunctor::operator()(amrex::MultiFab& mf_dst, const int dcomp, const int /*i_
     // output Multifab, mf_dst, the guard-cell data is not needed especially considering
     // the operations performend in the CoarsenAndInterpolate function.
     constexpr int ng = 1;
-#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD)
-    // For RZ spectral, all quantities are cell centered.
-    amrex::IntVect cell_type = amrex::IntVect::TheCellVector();
-#else
     // For staggered and nodal calculations, divE is computed on the nodes.
     // The temporary divE MultiFab is generated to comply with the location of divE.
     amrex::IntVect cell_type = amrex::IntVect::TheNodeVector();
+#ifdef WARPX_DIM_RZ
+    // For RZ spectral, all quantities are cell centered.
+    if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD)
+        cell_type = amrex::IntVect::TheCellVector();
 #endif
     const amrex::BoxArray& ba = amrex::convert(warpx.boxArray(m_lev), cell_type);
     amrex::MultiFab divE(ba, warpx.DistributionMap(m_lev), 2*warpx.n_rz_azimuthal_modes-1, ng );
diff --git a/Source/Diagnostics/FullDiagnostics.cpp b/Source/Diagnostics/FullDiagnostics.cpp
index 6f214b4d7..d00720a54 100644
--- a/Source/Diagnostics/FullDiagnostics.cpp
+++ b/Source/Diagnostics/FullDiagnostics.cpp
@@ -393,13 +393,15 @@ FullDiagnostics::InitializeFieldFunctors (int lev)
             m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_current_fp(lev, 2), lev, m_crse_ratio);
         } else if ( m_varnames[comp] == "rho" ){
             if ( WarpX::do_back_transformed_diagnostics ) {
-#ifdef WARPX_USE_PSATD
-                // rho_new is stored in component 1 of rho_fp when using PSATD
-                amrex::MultiFab* rho_new = new amrex::MultiFab(*warpx.get_pointer_rho_fp(lev), amrex::make_alias, 1, 1);
-                m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(rho_new, lev, m_crse_ratio);
-#else
-                m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_rho_fp(lev), lev, m_crse_ratio);
-#endif
+                if ( WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD ) {
+                    // rho_new is stored in component 1 of rho_fp when using PSATD
+                    amrex::MultiFab *rho_new = new amrex::MultiFab(*warpx.get_pointer_rho_fp(lev), amrex::make_alias,
+                                                                   1, 1);
+                    m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(rho_new, lev, m_crse_ratio);
+                } else {
+                    m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_rho_fp(lev),
+                                                                                          lev, m_crse_ratio);
+                }
             }
             else {
                 // Initialize rho functor to dump total rho
diff --git a/Source/Diagnostics/WarpXOpenPMD.cpp b/Source/Diagnostics/WarpXOpenPMD.cpp
index d7072ac56..a82d5836a 100644
--- a/Source/Diagnostics/WarpXOpenPMD.cpp
+++ b/Source/Diagnostics/WarpXOpenPMD.cpp
@@ -719,15 +719,12 @@ WarpXOpenPMDPlot::WriteOpenPMDFields( //const std::string& filename,
 
   auto meshes = series_iteration.meshes;
   meshes.setAttribute( "fieldSolver", [](){
-#ifdef WARPX_USE_PSATD
-      return "PSATD";
-#else
-      switch( WarpX::particle_pusher_algo ) {
+      switch( WarpX::maxwell_solver_id ) {
           case MaxwellSolverAlgo::Yee : return "Yee";
           case MaxwellSolverAlgo::CKC : return "CK";
+          case MaxwellSolverAlgo::PSATD : return "PSATD";
           default: return "other";
       }
-#endif
   }() );
   meshes.setAttribute( "fieldBoundary", fieldBoundary );
   meshes.setAttribute( "particleBoundary", particleBoundary );
diff --git a/Source/Evolve/WarpXComputeDt.cpp b/Source/Evolve/WarpXComputeDt.cpp
index a58fa5e7b..6237f85c4 100644
--- a/Source/Evolve/WarpXComputeDt.cpp
+++ b/Source/Evolve/WarpXComputeDt.cpp
@@ -23,42 +23,37 @@ WarpX::ComputeDt ()
     const amrex::Real* dx = geom[max_level].CellSize();
     amrex::Real deltat = 0.;
 
-#if WARPX_USE_PSATD
-    // Computation of dt for spectral algorithm
-
-#   if (defined WARPX_DIM_RZ)
-    // - In RZ geometry: dz/c
-    deltat = cfl * dx[1]/PhysConst::c;
-#   elif (defined WARPX_DIM_XZ)
-    // - In Cartesian 2D geometry: determined by the minimum cell size in all direction
-    deltat = cfl * std::min( dx[0], dx[1] )/PhysConst::c;
-#   else
-    // - In Cartesian 3D geometry: determined by the minimum cell size in all direction
-    deltat = cfl * std::min( dx[0], std::min( dx[1], dx[2] ) )/PhysConst::c;
-#   endif
-
-
+    if (maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+        // Computation of dt for spectral algorithm
+#if (defined WARPX_DIM_RZ)
+        // - In RZ geometry: dz/c
+        deltat = cfl * dx[1]/PhysConst::c;
+#elif (defined WARPX_DIM_XZ)
+        // - In Cartesian 2D geometry: determined by the minimum cell size in all direction
+        deltat = cfl * std::min( dx[0], dx[1] )/PhysConst::c;
 #else
-    // Computation of dt for FDTD algorithm
-
-#   ifdef WARPX_DIM_RZ
-    // - In RZ geometry
-    if (maxwell_solver_id == MaxwellSolverAlgo::Yee) {
-        deltat = cfl * CylindricalYeeAlgorithm::ComputeMaxDt(dx,  n_rz_azimuthal_modes);
-#   else
-    // - In Cartesian geometry
-    if (do_nodal) {
-        deltat = cfl * CartesianNodalAlgorithm::ComputeMaxDt( dx );
-    } else if (maxwell_solver_id == MaxwellSolverAlgo::Yee) {
-        deltat = cfl * CartesianYeeAlgorithm::ComputeMaxDt( dx );
-    } else if (maxwell_solver_id == MaxwellSolverAlgo::CKC) {
-        deltat = cfl * CartesianCKCAlgorithm::ComputeMaxDt( dx );
-#   endif
+        // - In Cartesian 3D geometry: determined by the minimum cell size in all direction
+        deltat = cfl * std::min(dx[0], std::min(dx[1], dx[2])) / PhysConst::c;
+#endif
     } else {
-        amrex::Abort("Unknown algorithm");
-    }
-
+        // Computation of dt for FDTD algorithm
+#ifdef WARPX_DIM_RZ
+        // - In RZ geometry
+        if (maxwell_solver_id == MaxwellSolverAlgo::Yee) {
+            deltat = cfl * CylindricalYeeAlgorithm::ComputeMaxDt(dx,  n_rz_azimuthal_modes);
+#else
+        // - In Cartesian geometry
+        if (do_nodal) {
+            deltat = cfl * CartesianNodalAlgorithm::ComputeMaxDt(dx);
+        } else if (maxwell_solver_id == MaxwellSolverAlgo::Yee) {
+            deltat = cfl * CartesianYeeAlgorithm::ComputeMaxDt(dx);
+        } else if (maxwell_solver_id == MaxwellSolverAlgo::CKC) {
+            deltat = cfl * CartesianCKCAlgorithm::ComputeMaxDt(dx);
 #endif
+        } else {
+            amrex::Abort("ComputeDt: Unknown algorithm");
+        }
+    }
 
     dt.resize(0);
     dt.resize(max_level+1,deltat);
diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp
index bf25e4263..a037c973c 100644
--- a/Source/Evolve/WarpXEvolve.cpp
+++ b/Source/Evolve/WarpXEvolve.cpp
@@ -56,9 +56,8 @@ WarpX::Evolve (int numsteps)
 
         amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(0);
         if (cost) {
-#ifdef WARPX_USE_PSATD
-            amrex::Abort("LoadBalance for PSATD: TODO");
-#endif
+            if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD)
+                amrex::Abort("LoadBalance for PSATD: TODO");
             if (step > 0 && load_balance_intervals.contains(step+1))
             {
                 LoadBalance();
@@ -114,10 +113,9 @@ WarpX::Evolve (int numsteps)
                 FillBoundaryE_avg(guard_cells.ng_FieldGather, guard_cells.ng_Extra);
                 FillBoundaryB_avg(guard_cells.ng_FieldGather, guard_cells.ng_Extra);
             }
-#ifndef WARPX_USE_PSATD
             // TODO Remove call to FillBoundaryAux before UpdateAuxilaryData?
-            FillBoundaryAux(guard_cells.ng_UpdateAux);
-#endif
+            if (WarpX::maxwell_solver_id != MaxwellSolverAlgo::PSATD)
+                FillBoundaryAux(guard_cells.ng_UpdateAux);
             UpdateAuxilaryData();
             FillBoundaryAux(guard_cells.ng_UpdateAux);
         }
@@ -288,14 +286,16 @@ WarpX::OneStep_nosub (Real cur_time)
 // TODO
 // Apply current correction in Fourier space: for domain decomposition with local
 // FFTs over guard cells, apply this before calling SyncCurrent
-#ifdef WARPX_USE_PSATD
-    if ( !fft_periodic_single_box && current_correction )
-        amrex::Abort("\nCurrent correction does not guarantee charge conservation with local FFTs over guard cells:\n"
-                     "set psatd.periodic_single_box_fft=1 too, in order to guarantee charge conservation");
-    if ( !fft_periodic_single_box && (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) )
-        amrex::Abort("\nVay current deposition does not guarantee charge conservation with local FFTs over guard cells:\n"
-                     "set psatd.periodic_single_box_fft=1 too, in order to guarantee charge conservation");
-#endif
+    if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+        if (!fft_periodic_single_box && current_correction)
+            amrex::Abort(
+                    "\nCurrent correction does not guarantee charge conservation with local FFTs over guard cells:\n"
+                    "set psatd.periodic_single_box_fft=1 too, in order to guarantee charge conservation");
+        if (!fft_periodic_single_box && (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay))
+            amrex::Abort(
+                    "\nVay current deposition does not guarantee charge conservation with local FFTs over guard cells:\n"
+                    "set psatd.periodic_single_box_fft=1 too, in order to guarantee charge conservation");
+    }
 
 #ifdef WARPX_QED
     doQEDEvents();
@@ -309,13 +309,13 @@ WarpX::OneStep_nosub (Real cur_time)
     SyncCurrent();
     SyncRho();
 
-// Apply current correction in Fourier space: for periodic single-box global FFTs
-// without guard cells, apply this after calling SyncCurrent
-#ifdef WARPX_USE_PSATD
-    if ( fft_periodic_single_box && current_correction ) CurrentCorrection();
-    if ( fft_periodic_single_box && (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) )
-        VayDeposition();
-#endif
+    // Apply current correction in Fourier space: for periodic single-box global FFTs
+    // without guard cells, apply this after calling SyncCurrent
+    if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+        if (fft_periodic_single_box && current_correction) CurrentCorrection();
+        if (fft_periodic_single_box && (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay))
+            VayDeposition();
+    }
 
 
     // At this point, J is up-to-date inside the domain, and E and B are
@@ -330,7 +330,7 @@ WarpX::OneStep_nosub (Real cur_time)
     // Electromagnetic solver:
     // Push E and B from {n} to {n+1}
     // (And update guard cells immediately afterwards)
-#ifdef WARPX_USE_PSATD
+    if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
         if (use_hybrid_QED)
         {
             WarpX::Hybrid_QED_Push(dt);
@@ -344,13 +344,12 @@ WarpX::OneStep_nosub (Real cur_time)
         {
             WarpX::Hybrid_QED_Push(dt);
             FillBoundaryE(guard_cells.ng_alloc_EB, guard_cells.ng_Extra);
-
         }
         if (do_pml) DampPML();
-#else
-        EvolveF(0.5*dt[0], DtType::FirstHalf);
+    } else {
+        EvolveF(0.5 * dt[0], DtType::FirstHalf);
         FillBoundaryF(guard_cells.ng_FieldSolverF);
-        EvolveB(0.5*dt[0]); // We now have B^{n+1/2}
+        EvolveB(0.5 * dt[0]); // We now have B^{n+1/2}
 
         FillBoundaryB(guard_cells.ng_FieldSolver, IntVect::TheZeroVector());
         if (WarpX::em_solver_medium == MediumForEM::Vacuum) {
@@ -364,8 +363,8 @@ WarpX::OneStep_nosub (Real cur_time)
         }
 
         FillBoundaryE(guard_cells.ng_FieldSolver, IntVect::TheZeroVector());
-        EvolveF(0.5*dt[0], DtType::SecondHalf);
-        EvolveB(0.5*dt[0]); // We now have B^{n+1}
+        EvolveF(0.5 * dt[0], DtType::SecondHalf);
+        EvolveB(0.5 * dt[0]); // We now have B^{n+1}
         if (do_pml) {
             FillBoundaryF(guard_cells.ng_alloc_F);
             DampPML();
@@ -375,10 +374,10 @@ WarpX::OneStep_nosub (Real cur_time)
         }
         // E and B are up-to-date in the domain, but all guard cells are
         // outdated.
-        if ( safe_guard_cells )
+        if (safe_guard_cells)
             FillBoundaryB(guard_cells.ng_alloc_EB, guard_cells.ng_Extra);
-#endif
-    }
+    } // !PSATD
+    } // !do_electrostatic
 }
 
 /* /brief Perform one PIC iteration, with subcycling
@@ -735,27 +734,45 @@ WarpX::applyMirrors(Real time){
 }
 
 // Apply current correction in Fourier space
-#ifdef WARPX_USE_PSATD
 void
 WarpX::CurrentCorrection ()
 {
-    for ( int lev = 0; lev <= finest_level; ++lev )
+#ifdef WARPX_USE_PSATD
+    if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD)
     {
-        spectral_solver_fp[lev]->CurrentCorrection( current_fp[lev], rho_fp[lev] );
-        if ( spectral_solver_cp[lev] ) spectral_solver_cp[lev]->CurrentCorrection( current_cp[lev], rho_cp[lev] );
+        for ( int lev = 0; lev <= finest_level; ++lev )
+        {
+            spectral_solver_fp[lev]->CurrentCorrection( current_fp[lev], rho_fp[lev] );
+            if ( spectral_solver_cp[lev] ) spectral_solver_cp[lev]->CurrentCorrection( current_cp[lev], rho_cp[lev] );
+        }
+    } else {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE( false,
+            "WarpX::CurrentCorrection: only implemented for spectral solver.");
     }
-}
+#else
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE( false,
+    "WarpX::CurrentCorrection: requires WarpX build with spectral solver support.");
 #endif
+}
 
 // Compute current from Vay deposition in Fourier space
-#ifdef WARPX_USE_PSATD
 void
 WarpX::VayDeposition ()
 {
-    for (int lev = 0; lev <= finest_level; ++lev)
+#ifdef WARPX_USE_PSATD
+    if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD)
     {
-        spectral_solver_fp[lev]->VayDeposition(current_fp[lev]);
-        if (spectral_solver_cp[lev]) spectral_solver_cp[lev]->VayDeposition(current_cp[lev]);
+        for (int lev = 0; lev <= finest_level; ++lev)
+        {
+            spectral_solver_fp[lev]->VayDeposition(current_fp[lev]);
+            if (spectral_solver_cp[lev]) spectral_solver_cp[lev]->VayDeposition(current_cp[lev]);
+        }
+    } else {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE( false,
+            "WarpX::VayDeposition: only implemented for spectral solver.");
     }
-}
+#else
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE( false,
+    "WarpX::CurrentCorrection: requires WarpX build with spectral solver support.");
 #endif
+}
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/ComputeDivE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/ComputeDivE.cpp
index a6376e03c..baadd0cfc 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/ComputeDivE.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/ComputeDivE.cpp
@@ -49,7 +49,7 @@ void FiniteDifferenceSolver::ComputeDivE (
 
 #endif
     } else {
-        amrex::Abort("Unknown algorithm");
+        amrex::Abort("ComputeDivE: Unknown algorithm");
     }
 
 }
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
index 9afec55e0..0ca576da2 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
@@ -48,7 +48,7 @@ void FiniteDifferenceSolver::EvolveB (
 
 #endif
     } else {
-        amrex::Abort("Unknown algorithm");
+        amrex::Abort("EvolveB: Unknown algorithm");
     }
 
 }
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp
index d449333c1..8346eb067 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp
@@ -47,7 +47,7 @@ void FiniteDifferenceSolver::EvolveBPML (
         EvolveBPMLCartesian <CartesianCKCAlgorithm> ( Bfield, Efield, dt );
 
     } else {
-        amrex::Abort("Unknown algorithm");
+        amrex::Abort("EvolveBPML: Unknown algorithm");
     }
 #endif
 }
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp
index fc43fe209..db1d1fe1b 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp
@@ -52,7 +52,7 @@ void FiniteDifferenceSolver::EvolveE (
 
 #endif
     } else {
-        amrex::Abort("Unknown algorithm");
+        amrex::Abort("EvolveE: Unknown algorithm");
     }
 
 }
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp
index a06d98f9b..d5ef1179b 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveEPML.cpp
@@ -56,7 +56,7 @@ void FiniteDifferenceSolver::EvolveEPML (
             Efield, Bfield, Jfield, Ffield, sigba, dt, pml_has_particles );
 
     } else {
-        amrex::Abort("Unknown algorithm");
+        amrex::Abort("EvolveEPML: Unknown algorithm");
     }
 #endif
 }
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveF.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveF.cpp
index 6627d2ef6..74499661e 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveF.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveF.cpp
@@ -52,7 +52,7 @@ void FiniteDifferenceSolver::EvolveF (
 
 #endif
     } else {
-        amrex::Abort("Unknown algorithm");
+        amrex::Abort("EvolveF: Unknown algorithm");
     }
 
 }
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveFPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveFPML.cpp
index ba7a90483..cdcab750d 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveFPML.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveFPML.cpp
@@ -47,7 +47,7 @@ void FiniteDifferenceSolver::EvolveFPML (
         EvolveFPMLCartesian <CartesianCKCAlgorithm> ( Ffield, Efield, dt );
 
     } else {
-        amrex::Abort("Unknown algorithm");
+        amrex::Abort("EvolveFPML: Unknown algorithm");
     }
 #endif
 }
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp
index fcab20510..3a752ca4b 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp
@@ -26,6 +26,10 @@ FiniteDifferenceSolver::FiniteDifferenceSolver (
     m_fdtd_algo = fdtd_algo;
     m_do_nodal = do_nodal;
 
+    // return if not FDTD
+    if (fdtd_algo == MaxwellSolverAlgo::PSATD)
+        return;
+
     // Calculate coefficients of finite-difference stencil
 #ifdef WARPX_DIM_RZ
     m_dr = cell_size[0];
@@ -46,7 +50,7 @@ FiniteDifferenceSolver::FiniteDifferenceSolver (
                               m_stencil_coefs_z.begin());
         amrex::Gpu::synchronize();
     } else {
-        amrex::Abort("Unknown algorithm");
+        amrex::Abort("FiniteDifferenceSolver: Unknown algorithm");
     }
 #else
     amrex::Vector<amrex::Real> stencil_coefs_x, stencil_coefs_y, stencil_coefs_z;
@@ -67,7 +71,7 @@ FiniteDifferenceSolver::FiniteDifferenceSolver (
             stencil_coefs_x, stencil_coefs_y, stencil_coefs_z );
 
     } else {
-        amrex::Abort("Unknown algorithm");
+        amrex::Abort("FiniteDifferenceSolver: Unknown algorithm");
     }
 
     m_stencil_coefs_x.resize(stencil_coefs_x.size());
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp
index 8ed0f5927..4a2d940aa 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/MacroscopicEvolveE.cpp
@@ -63,7 +63,7 @@ void FiniteDifferenceSolver::MacroscopicEvolveE (
         }
 
     } else {
-        amrex::Abort("Unknown algorithm");
+        amrex::Abort("MacroscopicEvolveE: Unknown algorithm");
     }
 #endif
 
diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp
index 7670a277b..168384413 100644
--- a/Source/FieldSolver/WarpXPushFieldsEM.cpp
+++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp
@@ -114,10 +114,15 @@ namespace {
 #endif
     }
 }
+#endif
 
 void
 WarpX::PushPSATD (amrex::Real a_dt)
 {
+#ifndef WARPX_USE_PSATD
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false,
+                                     "PushFieldsEM: PSATD solver selected but not built.");
+#else
     for (int lev = 0; lev <= finest_level; ++lev) {
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(dt[lev] == a_dt, "dt must be consistent");
         PushPSATD(lev, a_dt);
@@ -127,11 +132,19 @@ WarpX::PushPSATD (amrex::Real a_dt)
             pml[lev]->PushPSATD();
         }
     }
+#endif
 }
 
 void
-WarpX::PushPSATD (int lev, amrex::Real /* dt */)
-{
+WarpX::PushPSATD (int lev, amrex::Real /* dt */) {
+#ifndef WARPX_USE_PSATD
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false,
+                                     "PushFieldsEM: PSATD solver selected but not built.");
+#else
+    if (WarpX::maxwell_solver_id != MaxwellSolverAlgo::PSATD) {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false,
+                                         "WarpX::PushPSATD: only supported for PSATD solver.");
+    }
     // Update the fields on the fine and coarse patch
     PushPSATDSinglePatch( *spectral_solver_fp[lev],
         Efield_fp[lev], Bfield_fp[lev], Efield_avg_fp[lev], Bfield_avg_fp[lev], current_fp[lev], rho_fp[lev] );
@@ -142,8 +155,8 @@ WarpX::PushPSATD (int lev, amrex::Real /* dt */)
     if (use_damp_fields_in_z_guard) {
         DampFieldsInGuards( Efield_fp[lev], Bfield_fp[lev] );
     }
-}
 #endif
+}
 
 void
 WarpX::EvolveB (amrex::Real a_dt)
diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp
index 14acc5647..571848942 100644
--- a/Source/Initialization/WarpXInitData.cpp
+++ b/Source/Initialization/WarpXInitData.cpp
@@ -171,9 +171,8 @@ WarpX::InitPML ()
 #endif
         pml[0] = std::make_unique<PML>(boxArray(0), DistributionMap(0), &Geom(0), nullptr,
                              pml_ncell, pml_delta, 0,
-#ifdef WARPX_USE_PSATD
+                             // PSATD
                              dt[0], nox_fft, noy_fft, noz_fft, do_nodal,
-#endif
                              do_dive_cleaning, do_moving_window,
                              pml_has_particles, do_pml_in_domain,
                              do_pml_Lo_corrected, do_pml_Hi);
@@ -189,9 +188,8 @@ WarpX::InitPML ()
             pml[lev] = std::make_unique<PML>(boxArray(lev), DistributionMap(lev),
                                    &Geom(lev), &Geom(lev-1),
                                    pml_ncell, pml_delta, refRatio(lev-1)[0],
-#ifdef WARPX_USE_PSATD
+                                   // PSATD
                                    dt[lev], nox_fft, noy_fft, noz_fft, do_nodal,
-#endif
                                    do_dive_cleaning, do_moving_window,
                                    pml_has_particles, do_pml_in_domain,
                                    do_pml_Lo_MR, amrex::IntVect::TheUnitVector());
@@ -254,9 +252,9 @@ WarpX::InitFilter (){
 void
 WarpX::PostRestart ()
 {
-#ifdef WARPX_USE_PSATD
-    amrex::Abort("WarpX::PostRestart: TODO for PSATD");
-#endif
+    if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+        amrex::Abort("WarpX::PostRestart: TODO for PSATD");
+    }
     mypc->PostRestart();
 }
 
diff --git a/Source/Parallelization/GuardCellManager.cpp b/Source/Parallelization/GuardCellManager.cpp
index 85490f2e0..74e3bff24 100644
--- a/Source/Parallelization/GuardCellManager.cpp
+++ b/Source/Parallelization/GuardCellManager.cpp
@@ -8,14 +8,14 @@
 #include "Filter/NCIGodfreyFilter.H"
 #include "Utils/WarpXAlgorithmSelection.H"
 
-#include <AMReX_Print.H>
 #include <AMReX_ParmParse.H>
 #include <AMReX.H>
 
+
 using namespace amrex;
 
 void
-guardCellManager::Init(
+guardCellManager::Init (
     const bool do_subcycling,
     const bool do_fdtd_nci_corr,
     const bool do_nodal,
@@ -105,23 +105,23 @@ guardCellManager::Init(
     if (maxwell_solver_id == MaxwellSolverAlgo::CKC) ng_alloc_F_int = std::max( ng_alloc_F_int, 1 );
     ng_alloc_F = IntVect(AMREX_D_DECL(ng_alloc_F_int, ng_alloc_F_int, ng_alloc_F_int));
 
-#ifdef WARPX_USE_PSATD
-    // All boxes should have the same number of guard cells
-    // (to avoid temporary parallel copies)
-    // Thus take the max of the required number of guards for each field
-    // Also: the number of guard cell should be enough to contain
-    // the stencil of the FFT solver. Here, this number (`ngFFT`)
-    // is determined *empirically* to be the order of the solver
-    // for nodal, and half the order of the solver for staggered.
+    if (maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+        // All boxes should have the same number of guard cells
+        // (to avoid temporary parallel copies)
+        // Thus take the max of the required number of guards for each field
+        // Also: the number of guard cell should be enough to contain
+        // the stencil of the FFT solver. Here, this number (`ngFFT`)
+        // is determined *empirically* to be the order of the solver
+        // for nodal, and half the order of the solver for staggered.
 
-    int ngFFt_x = do_nodal ? nox_fft : nox_fft/2;
-    int ngFFt_y = do_nodal ? noy_fft : noy_fft/2;
-    int ngFFt_z = do_nodal ? noz_fft : noz_fft/2;
+        int ngFFt_x = do_nodal ? nox_fft : nox_fft / 2;
+        int ngFFt_y = do_nodal ? noy_fft : noy_fft / 2;
+        int ngFFt_z = do_nodal ? noz_fft : noz_fft / 2;
 
-    ParmParse pp("psatd");
-    pp.query("nx_guard", ngFFt_x);
-    pp.query("ny_guard", ngFFt_y);
-    pp.query("nz_guard", ngFFt_z);
+        ParmParse pp("psatd");
+        pp.query("nx_guard", ngFFt_x);
+        pp.query("ny_guard", ngFFt_y);
+        pp.query("nz_guard", ngFFt_z);
 
 #if (AMREX_SPACEDIM == 3)
         IntVect ngFFT = IntVect(ngFFt_x, ngFFt_y, ngFFt_z);
@@ -129,35 +129,33 @@ guardCellManager::Init(
         IntVect ngFFT = IntVect(ngFFt_x, ngFFt_z);
 #endif
 
-    for (int i_dim=0; i_dim<AMREX_SPACEDIM; i_dim++ ){
-        int ng_required = ngFFT[i_dim];
-        // Get the max
-        ng_required = std::max( ng_required, ng_alloc_EB[i_dim] );
-        ng_required = std::max( ng_required, ng_alloc_J[i_dim] );
-        ng_required = std::max( ng_required, ng_alloc_Rho[i_dim] );
-        ng_required = std::max( ng_required, ng_alloc_F[i_dim] );
-        // Set the guard cells to this max
-        ng_alloc_EB[i_dim] = ng_required;
-        ng_alloc_J[i_dim] = ng_required;
-        ng_alloc_F[i_dim] = ng_required;
-        ng_alloc_Rho[i_dim] = ng_required;
-        ng_alloc_F_int = ng_required;
+        for (int i_dim = 0; i_dim < AMREX_SPACEDIM; i_dim++) {
+            int ng_required = ngFFT[i_dim];
+            // Get the max
+            ng_required = std::max(ng_required, ng_alloc_EB[i_dim]);
+            ng_required = std::max(ng_required, ng_alloc_J[i_dim]);
+            ng_required = std::max(ng_required, ng_alloc_Rho[i_dim]);
+            ng_required = std::max(ng_required, ng_alloc_F[i_dim]);
+            // Set the guard cells to this max
+            ng_alloc_EB[i_dim] = ng_required;
+            ng_alloc_J[i_dim] = ng_required;
+            ng_alloc_F[i_dim] = ng_required;
+            ng_alloc_Rho[i_dim] = ng_required;
+            ng_alloc_F_int = ng_required;
+        }
+        ng_alloc_F = IntVect(AMREX_D_DECL(ng_alloc_F_int, ng_alloc_F_int, ng_alloc_F_int));
     }
-    ng_alloc_F = IntVect(AMREX_D_DECL(ng_alloc_F_int, ng_alloc_F_int, ng_alloc_F_int));
-#else
-    ignore_unused(nox_fft, noy_fft, noz_fft);
-#endif
 
     ng_Extra = IntVect(static_cast<int>(aux_is_nodal and !do_nodal));
 
     // Compute number of cells required for Field Solver
-#ifdef WARPX_USE_PSATD
-    ng_FieldSolver = ng_alloc_EB;
-    ng_FieldSolverF = ng_alloc_EB;
-#else
-    ng_FieldSolver = IntVect(AMREX_D_DECL(1,1,1));
-    ng_FieldSolverF = IntVect(AMREX_D_DECL(1,1,1));
-#endif
+    if (maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+        ng_FieldSolver = ng_alloc_EB;
+        ng_FieldSolverF = ng_alloc_EB;
+    } else {
+        ng_FieldSolver = IntVect(AMREX_D_DECL(1, 1, 1));
+        ng_FieldSolverF = IntVect(AMREX_D_DECL(1, 1, 1));
+    }
 
     if (safe_guard_cells){
         // Run in safe mode: exchange all allocated guard cells at each
diff --git a/Source/Parallelization/WarpXSumGuardCells.H b/Source/Parallelization/WarpXSumGuardCells.H
index 972c1cd2d..b9789b45c 100644
--- a/Source/Parallelization/WarpXSumGuardCells.H
+++ b/Source/Parallelization/WarpXSumGuardCells.H
@@ -16,22 +16,23 @@
  * This is typically called for the sources of the Maxwell equations (J/rho)
  * after deposition from the macroparticles.
  *
- *  - When WarpX is compiled with a finite-difference scheme: this only
+ *  - When WarpX is used with a finite-difference scheme: this only
  *    updates the *valid* cells of `mf`
- *  - When WarpX is compiled with a spectral scheme (WARPX_USE_PSATD): this
+ *  - When WarpX is used with a spectral scheme (PSATD): this
  *    updates both the *valid* cells and *guard* cells. (This is because a
  *    spectral solver requires the value of the sources over a large stencil.)
  */
 inline void
 WarpXSumGuardCells(amrex::MultiFab& mf, const amrex::Periodicity& period,
-                   const int icomp=0, const int ncomp=1){
-#ifdef WARPX_USE_PSATD
-   // Update both valid cells and guard cells
-   const amrex::IntVect n_updated_guards = mf.nGrowVect();
-#else
-   // Update only the valid cells
-   const amrex::IntVect n_updated_guards = amrex::IntVect::TheZeroVector();
-#endif
+                   const int icomp=0, const int ncomp=1)
+{
+    amrex::IntVect n_updated_guards;
+
+    // Update both valid cells and guard cells
+    if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD)
+        n_updated_guards = mf.nGrowVect();
+    else  // Update only the valid cells
+        n_updated_guards = amrex::IntVect::TheZeroVector();
     mf.SumBoundary(icomp, ncomp, n_updated_guards, period);
 }
 
@@ -41,9 +42,9 @@ WarpXSumGuardCells(amrex::MultiFab& mf, const amrex::Periodicity& period,
  * This is typically called for the sources of the Maxwell equations (J/rho)
  * after deposition from the macroparticles + filtering.
  *
- *  - When WarpX is compiled with a finite-difference scheme: this only
+ *  - When WarpX is used with a finite-difference scheme: this only
  *    updates the *valid* cells of `dst`
- *  - When WarpX is compiled with a spectral scheme (WARPX_USE_PSATD): this
+ *  - When WarpX is used with a spectral scheme (PSATD): this
  *    updates both the *valid* cells and *guard* cells. (This is because a
  *    spectral solver requires the value of the sources over a large stencil.)
  *
@@ -53,14 +54,16 @@ WarpXSumGuardCells(amrex::MultiFab& mf, const amrex::Periodicity& period,
 inline void
 WarpXSumGuardCells(amrex::MultiFab& dst, amrex::MultiFab& src,
                    const amrex::Periodicity& period,
-                   const int icomp=0, const int ncomp=1){
-#ifdef WARPX_USE_PSATD
+                   const int icomp=0, const int ncomp=1)
+{
+    amrex::IntVect n_updated_guards;
+
     // Update both valid cells and guard cells
-    const amrex::IntVect n_updated_guards = dst.nGrowVect();
-#else
-    // Update only the valid cells
-    const amrex::IntVect n_updated_guards = amrex::IntVect::TheZeroVector();
-#endif
+    if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD)
+        n_updated_guards = dst.nGrowVect();
+    else  // Update only the valid cells
+        n_updated_guards = amrex::IntVect::TheZeroVector();
+
     src.SumBoundary(0, ncomp, n_updated_guards, period);
     amrex::Copy( dst, src, 0, icomp, ncomp, n_updated_guards );
 }
diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp
index df39bd250..4c5bd469f 100644
--- a/Source/Particles/WarpXParticleContainer.cpp
+++ b/Source/Particles/WarpXParticleContainer.cpp
@@ -630,9 +630,14 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local)
     const auto& ba = m_gdb->ParticleBoxArray(lev);
     const auto& dm = m_gdb->DistributionMap(lev);
     BoxArray nba = ba;
-#if (!defined WARPX_DIM_RZ) || (!defined WARPX_USE_PSATD)
-    nba.surroundingNodes();
+
+    bool is_PSATD_RZ = false;
+#ifdef WARPX_DIM_RZ
+    if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD)
+        is_PSATD_RZ = true;
 #endif
+    if( !is_PSATD_RZ )
+        nba.surroundingNodes();
 
     // Number of guard cells for local deposition of rho
     WarpX& warpx = WarpX::GetInstance();
diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H
index ed3acfeba..c4d2be38a 100644
--- a/Source/Utils/WarpXAlgorithmSelection.H
+++ b/Source/Utils/WarpXAlgorithmSelection.H
@@ -38,7 +38,8 @@ struct MacroscopicSolverAlgo {
 struct MaxwellSolverAlgo {
     enum {
         Yee = 0,
-        CKC = 1
+        CKC = 1,
+        PSATD = 2
     };
 };
 
diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp
index de7944a7f..869f3b462 100644
--- a/Source/Utils/WarpXAlgorithmSelection.cpp
+++ b/Source/Utils/WarpXAlgorithmSelection.cpp
@@ -6,6 +6,7 @@
  *
  * License: BSD-3-Clause-LBNL
  */
+#include "WarpX.H"
 #include "WarpXAlgorithmSelection.H"
 
 #include <algorithm>
@@ -18,9 +19,8 @@
 
 const std::map<std::string, int> maxwell_solver_algo_to_int = {
     {"yee",     MaxwellSolverAlgo::Yee },
-#ifndef WARPX_DIM_RZ // Not available in RZ
     {"ckc",     MaxwellSolverAlgo::CKC },
-#endif
+    {"psatd",   MaxwellSolverAlgo::PSATD },
     {"default", MaxwellSolverAlgo::Yee }
 };
 
@@ -42,11 +42,7 @@ const std::map<std::string, int> current_deposition_algo_to_int = {
     {"esirkepov", CurrentDepositionAlgo::Esirkepov },
     {"direct",    CurrentDepositionAlgo::Direct },
     {"vay",       CurrentDepositionAlgo::Vay },
-#ifdef WARPX_USE_PSATD
-    {"default",   CurrentDepositionAlgo::Direct }
-#else
-    {"default",   CurrentDepositionAlgo::Esirkepov }
-#endif
+    {"default",   CurrentDepositionAlgo::Esirkepov } // NOTE: overwritten for PSATD below
 };
 
 const std::map<std::string, int> charge_deposition_algo_to_int = {
@@ -97,6 +93,8 @@ GetAlgorithmInteger( amrex::ParmParse& pp, const char* pp_search_key ){
         algo_to_int = particle_pusher_algo_to_int;
     } else if (0 == std::strcmp(pp_search_key, "current_deposition")) {
         algo_to_int = current_deposition_algo_to_int;
+        if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD)
+            algo_to_int["default"] = CurrentDepositionAlgo::Direct;
     } else if (0 == std::strcmp(pp_search_key, "charge_deposition")) {
         algo_to_int = charge_deposition_algo_to_int;
     } else if (0 == std::strcmp(pp_search_key, "field_gathering")) {
diff --git a/Source/Utils/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp
index 7b66b416e..b7a046eb4 100644
--- a/Source/Utils/WarpXUtil.cpp
+++ b/Source/Utils/WarpXUtil.cpp
@@ -5,9 +5,10 @@
  *
  * License: BSD-3-Clause-LBNL
  */
-#include "WarpXUtil.H"
-#include "WarpXConst.H"
 #include "WarpX.H"
+#include "WarpXAlgorithmSelection.H"
+#include "WarpXConst.H"
+#include "WarpXUtil.H"
 
 #include <AMReX_ParmParse.H>
 
@@ -270,7 +271,17 @@ getWithParser (const amrex::ParmParse& a_pp, char const * const str, amrex::Real
  */
 void CheckGriddingForRZSpectral ()
 {
-#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD)
+#ifndef WARPX_DIM_RZ
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false,
+        "CheckGriddingForRZSpectral: WarpX was not built with RZ geometry.");
+#endif
+
+    ParmParse pp("algo");
+    int maxwell_solver_id = GetAlgorithmInteger(pp, "maxwell_solver");
+
+    // only check for PSATD in RZ
+    if (maxwell_solver_id != MaxwellSolverAlgo::PSATD)
+        return;
 
     int max_level;
     Vector<int> n_cell(AMREX_SPACEDIM, -1);
@@ -336,8 +347,6 @@ void CheckGriddingForRZSpectral ()
         mg[0] /= 2;
     }
     pp_amr.addarr("max_grid_size_y", mg);
-
-#endif
 }
 
 namespace WarpXUtilMsg{
diff --git a/Source/WarpX.H b/Source/WarpX.H
index 949247ac9..3036dcf62 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -137,17 +137,13 @@ public:
     static int em_solver_medium;
     static int macroscopic_solver_algo;
 
-#ifdef WARPX_USE_PSATD
-    // If true (overwritten by the user in the input file), the current correction
+    // PSATD: If true (overwritten by the user in the input file), the current correction
     // defined in equation (19) of https://doi.org/10.1016/j.jcp.2013.03.010 is applied
     bool current_correction = false;
-#endif
 
-#ifdef WARPX_USE_PSATD
-    // If true, the update equation for E contains both J and rho (at times n and n+1):
+    // PSATD: If true, the update equation for E contains both J and rho (at times n and n+1):
     // default is false for standard PSATD and true for Galilean PSATD (set in WarpX.cpp)
-    bool update_with_rho;
-#endif
+    bool update_with_rho = false;
 
     // div E cleaning
     static int do_dive_cleaning;
@@ -881,14 +877,14 @@ private:
     // Domain decomposition on Level 0
     amrex::IntVect numprocs{0};
 
-#ifdef WARPX_USE_PSATD
 private:
-    void EvolvePSATD (int numsteps);
+    // void EvolvePSATD (int numsteps);
     void PushPSATD (amrex::Real dt);
     void PushPSATD (int lev, amrex::Real dt);
 
-    int fftw_plan_measure = 1;
+    int fftw_plan_measure = 1; // used with PSATD
 
+#ifdef WARPX_USE_PSATD
 #   ifdef WARPX_DIM_RZ
         amrex::Vector<std::unique_ptr<SpectralSolverRZ>> spectral_solver_fp;
         amrex::Vector<std::unique_ptr<SpectralSolverRZ>> spectral_solver_cp;
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 9b9304a28..6e9c05427 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -250,39 +250,39 @@ WarpX::WarpX ()
         && WarpX::load_balance_costs_update_algo==LoadBalanceCostsUpdateAlgo::Heuristic)
     {
 #ifdef AMREX_USE_GPU
-#ifdef WARPX_USE_PSATD
-        switch (WarpX::nox)
-        {
-            case 1:
-                costs_heuristic_cells_wt = 0.575;
-                costs_heuristic_particles_wt = 0.425;
-                break;
-            case 2:
-                costs_heuristic_cells_wt = 0.405;
-                costs_heuristic_particles_wt = 0.595;
-                break;
-            case 3:
-                costs_heuristic_cells_wt = 0.250;
-                costs_heuristic_particles_wt = 0.750;
-                break;
-        }
-#else // FDTD
-        switch (WarpX::nox)
-        {
-            case 1:
-                costs_heuristic_cells_wt = 0.401;
-                costs_heuristic_particles_wt = 0.599;
-                break;
-            case 2:
-                costs_heuristic_cells_wt = 0.268;
-                costs_heuristic_particles_wt = 0.732;
-                break;
-            case 3:
-                costs_heuristic_cells_wt = 0.145;
-                costs_heuristic_particles_wt = 0.855;
-                break;
+        if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+            switch (WarpX::nox)
+            {
+                case 1:
+                    costs_heuristic_cells_wt = 0.575;
+                    costs_heuristic_particles_wt = 0.425;
+                    break;
+                case 2:
+                    costs_heuristic_cells_wt = 0.405;
+                    costs_heuristic_particles_wt = 0.595;
+                    break;
+                case 3:
+                    costs_heuristic_cells_wt = 0.250;
+                    costs_heuristic_particles_wt = 0.750;
+                    break;
+            }
+        } else { // FDTD
+            switch (WarpX::nox)
+            {
+                case 1:
+                    costs_heuristic_cells_wt = 0.401;
+                    costs_heuristic_particles_wt = 0.599;
+                    break;
+                case 2:
+                    costs_heuristic_cells_wt = 0.268;
+                    costs_heuristic_particles_wt = 0.732;
+                    break;
+                case 3:
+                    costs_heuristic_cells_wt = 0.145;
+                    costs_heuristic_particles_wt = 0.855;
+                    break;
+            }
         }
-#endif // WARPX_USE_PSATD
 #else // CPU
         costs_heuristic_cells_wt = 0.1;
         costs_heuristic_particles_wt = 0.9;
@@ -291,11 +291,15 @@ WarpX::WarpX ()
 
     // Allocate field solver objects
 #ifdef WARPX_USE_PSATD
-    spectral_solver_fp.resize(nlevs_max);
-    spectral_solver_cp.resize(nlevs_max);
+    if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+        spectral_solver_fp.resize(nlevs_max);
+        spectral_solver_cp.resize(nlevs_max);
+    }
 #endif
-    m_fdtd_solver_fp.resize(nlevs_max);
-    m_fdtd_solver_cp.resize(nlevs_max);
+    if (WarpX::maxwell_solver_id != MaxwellSolverAlgo::PSATD) {
+        m_fdtd_solver_fp.resize(nlevs_max);
+        m_fdtd_solver_cp.resize(nlevs_max);
+    }
 
     // NCI Godfrey filters can have different stencils
     // at different levels (the stencil depends on c*dt/dz)
@@ -305,10 +309,10 @@ WarpX::WarpX ()
     // Sanity checks. Must be done after calling the MultiParticleContainer
     // constructor, as it reads additional parameters
     // (e.g., use_fdtd_nci_corr)
-#ifdef WARPX_USE_PSATD
-    AMREX_ALWAYS_ASSERT(use_fdtd_nci_corr == 0);
-    AMREX_ALWAYS_ASSERT(do_subcycling == 0);
-#endif
+    if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+        AMREX_ALWAYS_ASSERT(use_fdtd_nci_corr == 0);
+        AMREX_ALWAYS_ASSERT(do_subcycling == 0);
+    }
 }
 
 WarpX::~WarpX ()
@@ -610,25 +614,42 @@ WarpX::ReadParameters ()
 
         // Only needs to be set with WARPX_DIM_RZ, otherwise defaults to 1
         pp.query("n_rz_azimuthal_modes", n_rz_azimuthal_modes);
-
-#if defined WARPX_DIM_RZ
-        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(Geom(0).isPeriodic(0) == 0,
-                   "The problem must not be periodic in the radial direction");
-#endif
-#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD)
-        // Force do_nodal=true (that is, not staggered) and
-        // use same shape factors in all directions, for gathering
-        do_nodal = true;
-        galerkin_interpolation = false;
-#endif
     }
 
     {
         ParmParse pp("algo");
+        maxwell_solver_id = GetAlgorithmInteger(pp, "maxwell_solver");
+#ifdef WARPX_DIM_RZ
+        if (maxwell_solver_id == MaxwellSolverAlgo::CKC) {
+            AMREX_ALWAYS_ASSERT_WITH_MESSAGE( false,
+                "algo.maxwell_solver = ckc is not (yet) available for RZ geometry");
+        }
+#endif
+#ifndef WARPX_USE_PSATD
+        if (maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+            AMREX_ALWAYS_ASSERT_WITH_MESSAGE( false,
+                                              "algo.maxwell_solver = psatd is not supported because WarpX was built without spectral solvers");
+        }
+#endif
+
+#ifdef WARPX_DIM_RZ
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(Geom(0).isPeriodic(0) == 0,
+            "The problem must not be periodic in the radial direction");
+
+        if (maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+            // Force do_nodal=true (that is, not staggered) and
+            // use same shape factors in all directions, for gathering
+            do_nodal = true;
+            galerkin_interpolation = false;
+        }
+#endif
+
+        // note: current_deposition must be set after maxwell_solver is already determined,
+        //       because its default depends on the solver selection
         current_deposition_algo = GetAlgorithmInteger(pp, "current_deposition");
         charge_deposition_algo = GetAlgorithmInteger(pp, "charge_deposition");
         particle_pusher_algo = GetAlgorithmInteger(pp, "particle_pusher");
-        maxwell_solver_id = GetAlgorithmInteger(pp, "maxwell_solver");
+
         field_gathering_algo = GetAlgorithmInteger(pp, "field_gathering");
         if (field_gathering_algo == GatheringAlgo::MomentumConserving) {
             // Use same shape factors in all directions, for gathering
@@ -642,7 +663,6 @@ WarpX::ReadParameters ()
         queryWithParser(pp, "costs_heuristic_cells_wt", costs_heuristic_cells_wt);
         queryWithParser(pp, "costs_heuristic_particles_wt", costs_heuristic_particles_wt);
     }
-
     {
         ParmParse pp("interpolation");
         pp.query("nox", nox);
@@ -672,7 +692,7 @@ WarpX::ReadParameters ()
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE( nox >= 1, "warpx.nox must >= 1");
     }
 
-#ifdef WARPX_USE_PSATD
+    if (maxwell_solver_id == MaxwellSolverAlgo::PSATD)
     {
         ParmParse pp("psatd");
         pp.query("periodic_single_box_fft", fft_periodic_single_box);
@@ -756,7 +776,6 @@ WarpX::ReadParameters ()
 #   endif
 
     }
-#endif
 
     // for slice generation //
     {
@@ -996,21 +1015,21 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
         jz_nodal_flag  = IntVect::TheNodeVector();
         rho_nodal_flag = IntVect::TheNodeVector();
     }
-#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD)
-    // Force cell-centered IndexType in r and z
-    Ex_nodal_flag  = IntVect::TheCellVector();
-    Ey_nodal_flag  = IntVect::TheCellVector();
-    Ez_nodal_flag  = IntVect::TheCellVector();
-    Bx_nodal_flag  = IntVect::TheCellVector();
-    By_nodal_flag  = IntVect::TheCellVector();
-    Bz_nodal_flag  = IntVect::TheCellVector();
-    jx_nodal_flag  = IntVect::TheCellVector();
-    jy_nodal_flag  = IntVect::TheCellVector();
-    jz_nodal_flag  = IntVect::TheCellVector();
-    rho_nodal_flag = IntVect::TheCellVector();
-#endif
+#ifdef WARPX_DIM_RZ
+    if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) {
+        // Force cell-centered IndexType in r and z
+        Ex_nodal_flag  = IntVect::TheCellVector();
+        Ey_nodal_flag  = IntVect::TheCellVector();
+        Ez_nodal_flag  = IntVect::TheCellVector();
+        Bx_nodal_flag  = IntVect::TheCellVector();
+        By_nodal_flag  = IntVect::TheCellVector();
+        Bz_nodal_flag  = IntVect::TheCellVector();
+        jx_nodal_flag  = IntVect::TheCellVector();
+        jy_nodal_flag  = IntVect::TheCellVector();
+        jz_nodal_flag  = IntVect::TheCellVector();
+        rho_nodal_flag = IntVect::TheCellVector();
+    }
 
-#if defined WARPX_DIM_RZ
     // With RZ multimode, there is a real and imaginary component
     // for each mode, except mode 0 which is purely real
     // Component 0 is mode 0.
@@ -1073,40 +1092,49 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
     {
         F_fp[lev] = std::make_unique<MultiFab>(amrex::convert(ba,IntVect::TheUnitVector()),dm,ncomps, ngF.max());
     }
-#ifdef WARPX_USE_PSATD
-    // Allocate and initialize the spectral solver
+    bool const pml_flag_false = false;
+    if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD)
+    {
+        // Allocate and initialize the spectral solver
+#ifndef WARPX_USE_PSATD
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE( false,
+            "WarpX::AllocLevelMFs: PSATD solver requires WarpX build with spectral solver support.");
+#else
+        if (!do_dive_cleaning)
+            rho_fp[lev] = std::make_unique<MultiFab>(amrex::convert(ba,rho_nodal_flag),dm,2*ncomps,ngRho);
+
 #   if (AMREX_SPACEDIM == 3)
-    RealVect dx_vect(dx[0], dx[1], dx[2]);
+        RealVect dx_vect(dx[0], dx[1], dx[2]);
 #   elif (AMREX_SPACEDIM == 2)
-    RealVect dx_vect(dx[0], dx[2]);
+        RealVect dx_vect(dx[0], dx[2]);
 #   endif
-    // Check whether the option periodic, single box is valid here
-    if (fft_periodic_single_box) {
+        // Check whether the option periodic, single box is valid here
+        if (fft_periodic_single_box) {
 #   ifdef WARPX_DIM_RZ
-        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
-            geom[0].isPeriodic(1)          // domain is periodic in z
-            && ba.size() == 1 && lev == 0, // domain is decomposed in a single box
-            "The option `psatd.periodic_single_box_fft` can only be used for a periodic domain, decomposed in a single box");
+            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+                geom[0].isPeriodic(1)          // domain is periodic in z
+                && ba.size() == 1 && lev == 0, // domain is decomposed in a single box
+                "The option `psatd.periodic_single_box_fft` can only be used for a periodic domain, decomposed in a single box");
 #   else
-        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
-            geom[0].isAllPeriodic()        // domain is periodic in all directions
-            && ba.size() == 1 && lev == 0, // domain is decomposed in a single box
-            "The option `psatd.periodic_single_box_fft` can only be used for a periodic domain, decomposed in a single box");
+            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+                geom[0].isAllPeriodic()        // domain is periodic in all directions
+                && ba.size() == 1 && lev == 0, // domain is decomposed in a single box
+                "The option `psatd.periodic_single_box_fft` can only be used for a periodic domain, decomposed in a single box");
 #   endif
-    }
-    // Get the cell-centered box
-    BoxArray realspace_ba = ba;  // Copy box
-    realspace_ba.enclosedCells(); // Make it cell-centered
-    // Define spectral solver
+        }
+        // Get the cell-centered box
+        BoxArray realspace_ba = ba;  // Copy box
+        realspace_ba.enclosedCells(); // Make it cell-centered
+        // Define spectral solver
 #   ifdef WARPX_DIM_RZ
-    if ( fft_periodic_single_box == false ) {
-        realspace_ba.grow(1, ngE[1]); // add guard cells only in z
-    }
-    spectral_solver_fp[lev] = std::make_unique<SpectralSolverRZ>( realspace_ba, dm,
-        n_rz_azimuthal_modes, noz_fft, do_nodal, m_v_galilean, dx_vect, dt[lev], lev );
-    if (use_kspace_filter) {
-        spectral_solver_fp[lev]->InitFilter(filter_npass_each_dir, use_filter_compensation);
-    }
+        if ( fft_periodic_single_box == false ) {
+            realspace_ba.grow(1, ngE[1]); // add guard cells only in z
+        }
+        spectral_solver_fp[lev] = std::make_unique<SpectralSolverRZ>( realspace_ba, dm,
+            n_rz_azimuthal_modes, noz_fft, do_nodal, m_v_galilean, dx_vect, dt[lev], lev );
+        if (use_kspace_filter) {
+            spectral_solver_fp[lev]->InitFilter(filter_npass_each_dir, use_filter_compensation);
+        }
 #   else
     if ( fft_periodic_single_box == false ) {
         realspace_ba.grow(ngE); // add guard cells
@@ -1117,8 +1145,11 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
         pml_flag_false, fft_periodic_single_box, update_with_rho, fft_do_time_averaging );
 #   endif
 #endif
-    m_fdtd_solver_fp[lev] = std::make_unique<FiniteDifferenceSolver>(
-        maxwell_solver_id, dx, do_nodal);
+    } // MaxwellSolverAlgo::PSATD
+    else {
+        m_fdtd_solver_fp[lev] = std::make_unique<FiniteDifferenceSolver>(maxwell_solver_id, dx, do_nodal);
+    }
+
     //
     // The Aux patch (i.e., the full solution)
     //
@@ -1206,28 +1237,32 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
         {
             F_cp[lev] = std::make_unique<MultiFab>(amrex::convert(cba,IntVect::TheUnitVector()),dm,ncomps, ngF.max());
         }
-#ifdef WARPX_USE_PSATD
-        else
+        if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD)
         {
-            rho_cp[lev] = std::make_unique<MultiFab>(amrex::convert(cba,rho_nodal_flag),dm,2*ncomps,ngRho);
-        }
-        // Allocate and initialize the spectral solver
+            // Allocate and initialize the spectral solver
+#ifndef WARPX_USE_PSATD
+            AMREX_ALWAYS_ASSERT_WITH_MESSAGE( false,
+                "WarpX::AllocLevelMFs: PSATD solver requires WarpX build with spectral solver support.");
+#else
+            if (!do_dive_cleaning)
+                rho_cp[lev] = std::make_unique<MultiFab>( amrex::convert(cba,rho_nodal_flag),dm,2*ncomps,ngRho );
+
 #   if (AMREX_SPACEDIM == 3)
-        RealVect cdx_vect(cdx[0], cdx[1], cdx[2]);
+            RealVect cdx_vect(cdx[0], cdx[1], cdx[2]);
 #   elif (AMREX_SPACEDIM == 2)
-        RealVect cdx_vect(cdx[0], cdx[2]);
+            RealVect cdx_vect(cdx[0], cdx[2]);
 #   endif
-        // Get the cell-centered box, with guard cells
-        BoxArray c_realspace_ba = cba;// Copy box
-        c_realspace_ba.enclosedCells(); // Make it cell-centered
-        // Define spectral solver
+            // Get the cell-centered box, with guard cells
+            BoxArray c_realspace_ba = cba;// Copy box
+            c_realspace_ba.enclosedCells(); // Make it cell-centered
+            // Define spectral solver
 #   ifdef WARPX_DIM_RZ
-        c_realspace_ba.grow(1, ngE[1]); // add guard cells only in z
-        spectral_solver_cp[lev] = std::make_unique<SpectralSolverRZ>( c_realspace_ba, dm,
-            n_rz_azimuthal_modes, noz_fft, do_nodal, m_v_galilean, cdx_vect, dt[lev], lev );
-        if (use_kspace_filter) {
-            spectral_solver_cp[lev]->InitFilter(filter_npass_each_dir, use_filter_compensation);
-        }
+            c_realspace_ba.grow(1, ngE[1]); // add guard cells only in z
+            spectral_solver_cp[lev] = std::make_unique<SpectralSolverRZ>( c_realspace_ba, dm,
+                n_rz_azimuthal_modes, noz_fft, do_nodal, m_v_galilean, cdx_vect, dt[lev], lev );
+            if (use_kspace_filter) {
+                spectral_solver_cp[lev]->InitFilter(filter_npass_each_dir, use_filter_compensation);
+            }
 #   else
         c_realspace_ba.grow(ngE); // add guard cells
         spectral_solver_cp[lev] = std::make_unique<SpectralSolver>( c_realspace_ba, dm,
@@ -1235,8 +1270,11 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
             pml_flag_false, fft_periodic_single_box, update_with_rho, fft_do_time_averaging );
 #   endif
 #endif
-        m_fdtd_solver_cp[lev] = std::make_unique<FiniteDifferenceSolver>(
-            maxwell_solver_id, cdx, do_nodal);
+        } // MaxwellSolverAlgo::PSATD
+        else {
+            m_fdtd_solver_cp[lev] = std::make_unique<FiniteDifferenceSolver>(maxwell_solver_id, cdx,
+                                                                             do_nodal);
+        }
     }
 
     //
@@ -1436,11 +1474,15 @@ WarpX::ComputeDivB (amrex::MultiFab& divB, int const dcomp,
 void
 WarpX::ComputeDivE(amrex::MultiFab& divE, const int lev)
 {
+    if ( WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD ) {
 #ifdef WARPX_USE_PSATD
-    spectral_solver_fp[lev]->ComputeSpectralDivE( Efield_aux[lev], divE );
+        spectral_solver_fp[lev]->ComputeSpectralDivE( Efield_aux[lev], divE );
 #else
-    m_fdtd_solver_fp[lev]->ComputeDivE( Efield_aux[lev], divE );
+        amrex::Abort("ComputeDivE: PSATD requested but not compiled");
 #endif
+    } else {
+        m_fdtd_solver_fp[lev]->ComputeDivE( Efield_aux[lev], divE );
+    }
 }
 
 PML*
diff --git a/Source/main.cpp b/Source/main.cpp
index 0ac1e9358..7f4a7cc46 100644
--- a/Source/main.cpp
+++ b/Source/main.cpp
@@ -36,7 +36,9 @@ int main(int argc, char* argv[])
 
     ConvertLabParamsToBoost();
 
+#ifdef WARPX_DIM_RZ
     CheckGriddingForRZSpectral();
+#endif
 
     WARPX_PROFILE_VAR("main()", pmain);