diff --git a/examples/fokker-planck-1D2V/fokker-planck-1D2V-even_nz-shorttest-nstep200.toml b/examples/fokker-planck-1D2V/fokker-planck-1D2V-even_nz-shorttest-nstep200.toml
index eb7d60b18..7422ef0b0 100644
--- a/examples/fokker-planck-1D2V/fokker-planck-1D2V-even_nz-shorttest-nstep200.toml
+++ b/examples/fokker-planck-1D2V/fokker-planck-1D2V-even_nz-shorttest-nstep200.toml
@@ -9,7 +9,6 @@ evolve_moments_parallel_flow = false
 evolve_moments_parallel_pressure = false
 evolve_moments_conservation = false
 #force_Er_zero_at_wall = false #true
-#Er_constant = 0.0
 #epsilon_offset = 0.1
 #use_vpabar_in_mms_dfni = true
 T_e = 1.0
@@ -95,4 +94,4 @@ vperp_discretization = "gausslegendre_pseudospectral"
 [fokker_planck_collisions]
 use_fokker_planck = true
 nuii = 1.0
-frequency_option = "manual"
\ No newline at end of file
+frequency_option = "manual"
diff --git a/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas.toml b/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas.toml
new file mode 100644
index 000000000..e6b75e29c
--- /dev/null
+++ b/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas.toml
@@ -0,0 +1,92 @@
+# cheap input file for a 0D2V relaxation to a collisional Maxwellian distribution with self-ion collisions and collisions with fixed Maxwellian background of cold ions and electrons.
+n_ion_species = 1
+n_neutral_species = 0
+electron_physics = "boltzmann_electron_response"
+evolve_moments_density = false
+evolve_moments_parallel_flow = false
+evolve_moments_parallel_pressure = false
+evolve_moments_conservation = false
+T_e = 1.0
+T_wall = 1.0
+initial_density1 = 1.0
+initial_temperature1 = 1.0
+z_IC_option1 = "sinusoid"
+z_IC_density_amplitude1 = 0.0
+z_IC_density_phase1 = 0.0
+z_IC_upar_amplitude1 = 0.0
+z_IC_upar_phase1 = 0.0
+z_IC_temperature_amplitude1 = 0.0
+z_IC_temperature_phase1 = 0.0
+vpa_IC_option1 = "isotropic-beam"
+vpa_IC_v01 = 1.0
+vpa_IC_vth01 = 0.1
+#vpa_IC_option1 = "directed-beam"
+#vpa_IC_vpa01 = -1.5
+#vpa_IC_vperp01 = 0.0
+charge_exchange_frequency = 0.0
+ionization_frequency = 0.0
+constant_ionization_rate = false
+
+z_ngrid = 1
+z_nelement = 1
+z_nelement_local = 1
+z_bc = "wall"
+z_discretization = "chebyshev_pseudospectral"
+r_ngrid = 1
+r_nelement = 1
+r_nelement_local = 1
+r_bc = "periodic"
+r_discretization = "chebyshev_pseudospectral"
+vpa_ngrid = 5
+vpa_nelement = 32
+vpa_L = 3.0
+vpa_bc = "zero"
+vpa_discretization = "gausslegendre_pseudospectral"
+vperp_ngrid = 5
+vperp_nelement = 16
+vperp_L = 1.5
+vperp_discretization = "gausslegendre_pseudospectral"
+vperp_bc = "zero"
+# Fokker-Planck operator requires the "gausslegendre_pseudospectral
+# options for the vpa and vperp grids
+
+[fokker_planck_collisions]
+# nuii sets the normalised input C[F,F] Fokker-Planck collision frequency
+nuii = 1.876334222e-3 #(1/nu_alphae, as computed from input diagnostic)
+Zi = 2.0
+self_collisions = true
+slowing_down_test = true
+frequency_option = "manual"
+use_fokker_planck = true
+use_conserving_corrections = true
+sd_density = 1.0
+sd_temp = 0.025 # TD/Ealpha
+sd_mi = 0.5 # mD/malpha
+sd_me = 0.000013616 # 0.25/1836.0 me/malpha
+sd_q = 1.0
+
+[ion_source]
+active=false
+source_strength=0.0
+source_T=0.005
+source_n=1.0
+r_profile="constant"
+r_width=1.0
+r_relative_minimum=0.0
+z_profile="gaussian"
+z_width=0.1
+z_relative_minimum=0.0
+source_v0 = 1.0
+#source_type="alphas"
+source_type="alphas-with-losses"
+#source_type="beam-with-losses"
+#source_vpa0 = 1.0
+#source_vperp0 = 1.0
+sink_strength = 0.0
+sink_vth=0.07071067811865475                     
+
+[timestepping]
+nstep = 50000
+dt = 2.5e-4
+nwrite = 100
+nwrite_dfns = 100
diff --git a/examples/fokker-planck/fokker-planck-relaxation.toml b/examples/fokker-planck/fokker-planck-relaxation.toml
index f2d5651a1..0f633c12b 100644
--- a/examples/fokker-planck/fokker-planck-relaxation.toml
+++ b/examples/fokker-planck/fokker-planck-relaxation.toml
@@ -13,7 +13,7 @@ initial_temperature1 = 1.0
 initial_density2 = 0.5
 initial_temperature2 = 1.0
 z_IC_option1 = "sinusoid"
-z_IC_density_amplitude1 = 0.001
+z_IC_density_amplitude1 = 0.0
 z_IC_density_phase1 = 0.0
 z_IC_upar_amplitude1 = 0.0
 z_IC_upar_phase1 = 0.0
diff --git a/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas-no-self-collisions.toml b/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas-no-self-collisions.toml
new file mode 100644
index 000000000..1a95fc210
--- /dev/null
+++ b/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas-no-self-collisions.toml
@@ -0,0 +1,91 @@
+# cheap input file for a 0D2V relaxation to a collisional Maxwellian distribution with self-ion collisions and collisions with fixed Maxwellian background of cold ions and electrons.
+n_ion_species = 1
+n_neutral_species = 0
+electron_physics = "boltzmann_electron_response"
+evolve_moments_density = false
+evolve_moments_parallel_flow = false
+evolve_moments_parallel_pressure = false
+evolve_moments_conservation = false
+T_e = 1.0
+T_wall = 1.0
+initial_density1 = 1.0
+initial_temperature1 = 1.0
+z_IC_option1 = "sinusoid"
+z_IC_density_amplitude1 = 0.0
+z_IC_density_phase1 = 0.0
+z_IC_upar_amplitude1 = 0.0
+z_IC_upar_phase1 = 0.0
+z_IC_temperature_amplitude1 = 0.0
+z_IC_temperature_phase1 = 0.0
+vpa_IC_option1 = "isotropic-beam"
+vpa_IC_v01 = 1.0
+vpa_IC_vth01 = 0.1
+#vpa_IC_option1 = "directed-beam"
+#vpa_IC_vpa01 = -1.5
+#vpa_IC_vperp01 = 0.0
+charge_exchange_frequency = 0.0
+ionization_frequency = 0.0
+constant_ionization_rate = false
+
+z_ngrid = 1
+z_nelement = 1
+z_nelement_local = 1
+z_bc = "wall"
+z_discretization = "chebyshev_pseudospectral"
+r_ngrid = 1
+r_nelement = 1
+r_nelement_local = 1
+r_bc = "periodic"
+r_discretization = "chebyshev_pseudospectral"
+vpa_ngrid = 5
+vpa_nelement = 32
+vpa_L = 3.0
+vpa_bc = "zero"
+vpa_discretization = "gausslegendre_pseudospectral"
+vperp_ngrid = 5
+vperp_nelement = 16
+vperp_L = 1.5
+vperp_discretization = "gausslegendre_pseudospectral"
+vperp_bc = "zero"
+# Fokker-Planck operator requires the "gausslegendre_pseudospectral
+# options for the vpa and vperp grids
+
+[fokker_planck_collisions]
+# nuii sets the normalised input C[F,F] Fokker-Planck collision frequency
+nuii = 1.876334222e-3 #(1/nu_alphae, as computed from input diagnostic)
+Zi = 2.0
+self_collisions = false
+slowing_down_test = true
+frequency_option = "manual"
+use_fokker_planck = true
+sd_density = 1.0
+sd_temp = 0.0025 # TD/Ealpha
+sd_mi = 0.5 # mD/malpha
+sd_me = 0.000013616 # 0.25/1836.0 me/malpha
+sd_q = 1.0
+
+[ion_source]
+active=true
+source_strength=1.0
+source_T=0.005
+source_n=1.0
+r_profile="constant"
+r_width=1.0
+r_relative_minimum=0.0
+z_profile="gaussian"
+z_width=0.1
+z_relative_minimum=0.0
+source_v0 = 1.0
+#source_type="alphas"
+source_type="alphas-with-losses"
+#source_type="beam-with-losses"
+#source_vpa0 = 1.0
+#source_vperp0 = 1.0
+sink_strength = 1.0
+sink_vth=0.07071067811865475                     
+
+[timestepping]
+nstep = 250000
+dt = 1.0e-4
+nwrite = 500
+nwrite_dfns = 500
diff --git a/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas-source-sink.toml b/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas.toml
similarity index 97%
rename from examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas-source-sink.toml
rename to examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas.toml
index 098dc787a..cd49c82e7 100644
--- a/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas-source-sink.toml
+++ b/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas.toml
@@ -85,7 +85,7 @@ sink_strength = 1.0
 sink_vth=0.07071067811865475                     
 
 [timestepping]
-nstep = 50000
-dt = 5.0e-4
-nwrite = 100
-nwrite_dfns = 100
+nstep = 250000
+dt = 1.0e-4
+nwrite = 500
+nwrite_dfns = 500
diff --git a/makie_post_processing/makie_post_processing/src/shared_utils.jl b/makie_post_processing/makie_post_processing/src/shared_utils.jl
index 031bcda75..32acaf1ec 100644
--- a/makie_post_processing/makie_post_processing/src/shared_utils.jl
+++ b/makie_post_processing/makie_post_processing/src/shared_utils.jl
@@ -5,7 +5,6 @@ export calculate_and_write_frequencies, get_geometry, get_composition
 using moment_kinetics.analysis: fit_delta_phi_mode
 using moment_kinetics.array_allocation: allocate_float
 using moment_kinetics.coordinates: define_coordinate
-using moment_kinetics.geo: init_magnetic_geometry
 using moment_kinetics.input_structs: boltzmann_electron_response,
                                      boltzmann_electron_response_with_simple_sheath,
                                      grid_input, geometry_input, species_composition
@@ -13,7 +12,7 @@ using moment_kinetics.moment_kinetics_input: get_default_rhostar, setup_referenc
 using moment_kinetics.type_definitions: mk_float, mk_int
 using moment_kinetics.reference_parameters: setup_reference_parameters
 using moment_kinetics.moment_kinetics_input: get_default_rhostar
-using moment_kinetics.geo: init_magnetic_geometry
+using moment_kinetics.geo: init_magnetic_geometry, setup_geometry_input
 using MPI
 
 """
@@ -90,7 +89,6 @@ function get_composition(scan_input)
     use_test_neutral_wall_pdf = get(scan_input, "use_test_neutral_wall_pdf", false)
     gyrokinetic_ions = get(scan_input, "gyrokinetic_ions", false)
     # constant to be used to test nonzero Er in wall boundary condition
-    Er_constant = get(scan_input, "Er_constant", 0.0)
     recycling_fraction = get(scan_input, "recycling_fraction", 1.0)
     # constant to be used to control Ez divergences
     epsilon_offset = get(scan_input, "epsilon_offset", 0.001)
@@ -106,7 +104,7 @@ function get_composition(scan_input)
     # ratio of the electron particle mass to the ion particle mass
     me_over_mi = 1.0/1836.0
     composition = species_composition(n_species, n_ion_species, n_neutral_species,
-        electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, Er_constant,
+        electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall,
         mn_over_mi, me_over_mi, recycling_fraction, gyrokinetic_ions, allocate_float(n_species))
     return composition
 
@@ -114,18 +112,10 @@ end
 
 function get_geometry(scan_input,z,r)
     reference_params = setup_reference_parameters(scan_input)
-    # set geometry_input
-    # MRH need to get this in way that does not duplicate code
-    # MRH from moment_kinetics_input.jl
-    option = get(scan_input, "geometry_option", "constant-helical") #"1D-mirror"
-    pitch = get(scan_input, "pitch", 1.0)
-    rhostar = get(scan_input, "rhostar", get_default_rhostar(reference_params))
-    DeltaB = get(scan_input, "DeltaB", 1.0)
-    geo_in = geometry_input(rhostar,option,pitch,DeltaB)
+    reference_rhostar = get_default_rhostar(reference_params)
+    geo_in = setup_geometry_input(scan_input, reference_rhostar)
     geometry = init_magnetic_geometry(geo_in,z,r)
-    
     return geometry
-
 end
 
 end # shared_utils.jl
diff --git a/moment_kinetics/debug_test/runtests.jl b/moment_kinetics/debug_test/runtests.jl
index 94b5dbbc7..7317ee974 100644
--- a/moment_kinetics/debug_test/runtests.jl
+++ b/moment_kinetics/debug_test/runtests.jl
@@ -8,7 +8,7 @@ function runtests()
         include(joinpath(@__DIR__, "fokker_planck_collisions_tests.jl"))
         include(joinpath(@__DIR__, "wall_bc_tests.jl"))
         include(joinpath(@__DIR__, "harrisonthompson.jl"))
-        include(joinpath(@__DIR__, "mms_tests.jl"))
+        #include(joinpath(@__DIR__, "mms_tests.jl"))
         include(joinpath(@__DIR__, "restart_interpolation_tests.jl"))
         include(joinpath(@__DIR__, "recycling_fraction_tests.jl"))
         include(joinpath(@__DIR__, "gyroaverage_tests.jl"))
diff --git a/moment_kinetics/ext/manufactured_solns_ext.jl b/moment_kinetics/ext/manufactured_solns_ext.jl
index 0707af732..e65c964db 100644
--- a/moment_kinetics/ext/manufactured_solns_ext.jl
+++ b/moment_kinetics/ext/manufactured_solns_ext.jl
@@ -43,6 +43,8 @@ using IfElse
     # does not appear in a particular geometric coefficient, because 
     # that coefficient is a constant. 
     struct geometric_coefficients_sym{}
+        Er_constant::mk_float
+        Ez_constant::mk_float
         rhostar::mk_float
         Bzed::Union{mk_float,Num}
         Bzeta::Union{mk_float,Num}
@@ -62,6 +64,8 @@ using IfElse
         option = geometry_input_data.option
         rhostar = geometry_input_data.rhostar
         pitch = geometry_input_data.pitch
+        Er_constant = geometry_input_data.Er_constant
+        Ez_constant = geometry_input_data.Ez_constant
         if option == "constant-helical" || option == "default"
             bzed = pitch
             bzeta = sqrt(1 - bzed^2)
@@ -90,10 +94,20 @@ using IfElse
             end
             dBdz = expand_derivatives(Dz(Bmag))
             jacobian = 1.0
+        elseif option == "0D-Spitzer-test"
+            Bmag = 1.0
+            dBdz = geometry_input_data.dBdz_constant
+            dBdr = geometry_input_data.dBdr_constant
+            bzed = pitch
+            bzeta = sqrt(1 - bzed^2)
+            Bzed = Bmag*bzed
+            Bzeta = Bmag*bzeta
+            jacobian = 1.0
         else
             input_option_error("$option", option)
         end
-        geo_sym = geometric_coefficients_sym(rhostar,Bzed,Bzeta,Bmag,bzed,bzeta,dBdz,dBdr,jacobian)
+        geo_sym = geometric_coefficients_sym(Er_constant,Ez_constant,
+          rhostar,Bzed,Bzeta,Bmag,bzed,bzeta,dBdz,dBdr,jacobian)
         return geo_sym
     end
 
@@ -272,7 +286,7 @@ using IfElse
             upari = 0.0
         elseif z_bc == "wall"
             densi = densi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species)
-            Er, Ez, phi = electric_fields(Lr,Lz,r_bc,z_bc,composition,nr,manufactured_solns_input,species)
+            Er, Ez, phi = electric_fields(Lr,Lz,r_bc,z_bc,composition,geometry,nr,manufactured_solns_input,species)
             Bzeta = geometry.Bzeta
             Bmag = geometry.Bmag
             rhostar = geometry.rhostar
@@ -357,7 +371,7 @@ using IfElse
 
         if manufactured_solns_input.type == "default"
             # calculate the electric fields and the potential
-            Er, Ez, phi = electric_fields(Lr, Lz, r_bc, z_bc, composition, nr,
+            Er, Ez, phi = electric_fields(Lr, Lz, r_bc, z_bc, composition, geometry, nr,
                                           manufactured_solns_input, species)
 
             # get geometric/composition data
@@ -402,7 +416,7 @@ using IfElse
         return dfni
     end
 
-    function electric_fields(Lr, Lz, r_bc, z_bc, composition, nr,
+    function electric_fields(Lr, Lz, r_bc, z_bc, composition, geometry, nr,
                              manufactured_solns_input, species)
        
         # define derivative operators
@@ -412,7 +426,7 @@ using IfElse
         # get N_e factor for boltzmann response
         if composition.electron_physics == boltzmann_electron_response_with_simple_sheath && nr == 1 
             # so 1D MMS test with 3V neutrals where ion current can be calculated prior to knowing Er
-            jpari_into_LHS_wall = jpari_into_LHS_wall_sym(Lr, Lz, r_bc, z_bc, composition,
+            jpari_into_LHS_wall = jpari_into_LHS_wall_sym(Lr, Lz, r_bc, z_bc,
                                                           manufactured_solns_input)
             N_e = -2.0*sqrt(pi*composition.me_over_mi)*exp(-composition.phi_wall/composition.T_e)*jpari_into_LHS_wall
         elseif composition.electron_physics == boltzmann_electron_response_with_simple_sheath && nr > 1 
@@ -437,8 +451,8 @@ using IfElse
         # calculate the electric fields
         dense = densi # get the electron density via quasineutrality with Zi = 1
         phi = composition.T_e*log(dense/N_e) # use the adiabatic response of electrons for me/mi -> 0
-        Er = -Dr(phi)*rfac + composition.Er_constant
-        Ez = -Dz(phi)
+        Er = -Dr(phi)*rfac + geometry.Er_constant
+        Ez = -Dz(phi)      + geometry.Ez_constant
         
         Er_expanded = expand_derivatives(Er)
         Ez_expanded = expand_derivatives(Ez)
@@ -500,11 +514,13 @@ using IfElse
         return manufactured_solns_list
     end 
     
-    function manufactured_electric_fields(Lr, Lz, r_bc, z_bc, composition, nr,
+    function manufactured_electric_fields(Lr, Lz, r_bc, z_bc, composition, geometry_input_data::geometry_input, nr,
                                           manufactured_solns_input, species)
         
+        # calculate the geometry symbolically
+        geometry = geometry_sym(geometry_input_data,Lz,Lr,nr)
         # calculate the electric fields and the potential
-        Er, Ez, phi = electric_fields(Lr, Lz, r_bc, z_bc, composition, nr,
+        Er, Ez, phi = electric_fields(Lr, Lz, r_bc, z_bc, composition, geometry, nr,
                                       manufactured_solns_input, species)
         
         Er_func = build_function(Er, z, r, t, expression=Val{false})
@@ -603,7 +619,7 @@ using IfElse
         
         # calculate the electric fields and the potential
         Er, Ez, phi = electric_fields(r_coord.L, z_coord.L, r_coord.bc, z_coord.bc,
-                                      composition, r_coord.n, manufactured_solns_input,
+                                      composition, geometry, r_coord.n, manufactured_solns_input,
                                       ion_species)
 
         # the adiabatic invariant (for compactness)
diff --git a/moment_kinetics/src/boundary_conditions.jl b/moment_kinetics/src/boundary_conditions.jl
index 61dff1a08..b014abe99 100644
--- a/moment_kinetics/src/boundary_conditions.jl
+++ b/moment_kinetics/src/boundary_conditions.jl
@@ -954,6 +954,8 @@ function enforce_v_boundary_condition_local!(f, bc, speed, v_diffusion, v, v_spe
     elseif bc == "periodic"
         f[1] = 0.5*(f[1]+f[end])
         f[end] = f[1]
+    elseif bc == "none"
+        # Do nothing
     else
         error("Unsupported boundary condition option '$bc' for $(v.name)")
     end
@@ -972,7 +974,7 @@ function enforce_vperp_boundary_condition!(f::AbstractArray{mk_float,5}, bc, vpe
 end
 
 function enforce_vperp_boundary_condition!(f::AbstractArray{mk_float,4}, bc, vperp, vperp_spectral, vperp_advect, diffusion)
-    if bc == "zero"
+    if bc == "zero" || bc == "zero-impose-regularity"
         nvperp = vperp.n
         ngrid = vperp.ngrid
         # set zero boundary condition
@@ -982,7 +984,7 @@ function enforce_vperp_boundary_condition!(f::AbstractArray{mk_float,4}, bc, vpe
             end
         end
         # set regularity condition d F / d vperp = 0 at vperp = 0
-        if vperp.discretization == "gausslegendre_pseudospectral" || vperp.discretization == "chebyshev_pseudospectral"
+        if bc == "zero-impose-regularity" && (vperp.discretization == "gausslegendre_pseudospectral" || vperp.discretization == "chebyshev_pseudospectral")
             D0 = vperp_spectral.radau.D0
             buffer = @view vperp.scratch[1:ngrid-1]
             @loop_r_z_vpa ir iz ivpa begin
@@ -992,6 +994,8 @@ function enforce_vperp_boundary_condition!(f::AbstractArray{mk_float,4}, bc, vpe
                     f[ivpa,1,iz,ir] = -sum(buffer)/D0[1]
                 end
             end
+        elseif bc == "zero"
+            # do nothing
         else
             println("vperp.bc=\"$bc\" not supported by discretization "
                     * "$(vperp.discretization)")
diff --git a/moment_kinetics/src/calculus.jl b/moment_kinetics/src/calculus.jl
index d2c16a478..a7bda13e8 100644
--- a/moment_kinetics/src/calculus.jl
+++ b/moment_kinetics/src/calculus.jl
@@ -131,6 +131,70 @@ function second_derivative!(d2f, f, coord, spectral; handle_periodic=true)
     return nothing
 end
 
+function laplacian_derivative!(d2f, f, coord, spectral; handle_periodic=true)
+    # computes (1/coord) d / coord ( coord d f / d(coord)) for vperp coordinate
+    # For spectral element methods, calculate second derivative by applying first
+    # derivative twice, with special treatment for element boundaries
+
+    # First derivative
+    elementwise_derivative!(coord, f, spectral)
+    derivative_elements_to_full_grid!(coord.scratch3, coord.scratch_2d, coord)
+    if coord.name == "vperp"
+        # include the Jacobian
+        @. coord.scratch3 *= coord.grid
+    end
+    # MPI reconcile code here if used with z or r coords
+
+    # Save elementwise first derivative result
+    coord.scratch2_2d .= coord.scratch_2d
+
+    # Second derivative for element interiors
+    elementwise_derivative!(coord, coord.scratch3, spectral)
+    derivative_elements_to_full_grid!(d2f, coord.scratch_2d, coord)
+    if coord.name == "vperp"
+        # include the Jacobian
+        @. d2f /= coord.grid
+    end
+    # MPI reconcile code here if used with z or r coords
+
+    # Add contribution to penalise discontinuous first derivatives at element
+    # boundaries. For smooth functions this would do nothing so should not affect
+    # convergence of the second derivative. Aims to stabilise numerical instability when
+    # spike develops at an element boundary. The coefficient is an arbitrary choice, it
+    # should probably be large enough for stability but as small as possible.
+    #
+    # Arbitrary numerical coefficient
+    C = 1.0
+    function penalise_discontinuous_first_derivative!(d2f, imin, imax, df)
+        # Left element boundary
+        d2f[imin] += C * df[1]
+
+        # Right element boundary
+        d2f[imax] -= C * df[end]
+
+        return nothing
+    end
+    @views penalise_discontinuous_first_derivative!(d2f, 1, coord.imax[1],
+                                                    coord.scratch2_2d[:,1])
+    for ielement ∈ 2:coord.nelement_local
+        @views penalise_discontinuous_first_derivative!(d2f, coord.imin[ielement]-1,
+                                                        coord.imax[ielement],
+                                                        coord.scratch2_2d[:,ielement])
+    end
+
+    if coord.bc ∈ ("zero", "both_zero", "zero-no-regularity")
+        # For stability don't contribute to evolution at boundaries, in case these
+        # points are not set by a boundary condition.
+        # Full grid may be across processes and bc only applied to extreme ends of the
+        # domain.
+        if coord.irank == coord.nrank - 1
+            d2f[end] = 0.0
+        end
+    else
+        error("Unsupported bc '$(coord.bc)'")
+    end
+    return nothing
+end
 """
     mass_matrix_solve!(f, b, spectral::weak_discretization_info)
 
@@ -176,9 +240,11 @@ function laplacian_derivative!(d2f, f, coord, spectral::weak_discretization_info
     # for all other coord.name, do exactly the same as second_derivative! above.
     mul!(coord.scratch, spectral.L_matrix, f)
 
-    if handle_periodic && coord.bc == "periodic"
+    if coord.bc == "periodic" && coord.name == "vperp"
+        error("laplacian_derivative!() cannot handle periodic boundaries for vperp")
+    elseif coord.bc == "periodic"
         if coord.nrank > 1
-            error("second_derivative!() cannot handle periodic boundaries for a "
+            error("laplacian_derivative!() cannot handle periodic boundaries for a "
                   * "distributed coordinate")
         end
 
diff --git a/moment_kinetics/src/coordinates.jl b/moment_kinetics/src/coordinates.jl
index d034a2dd0..b761ea1fa 100644
--- a/moment_kinetics/src/coordinates.jl
+++ b/moment_kinetics/src/coordinates.jl
@@ -252,8 +252,7 @@ function define_coordinate(input, parallel_io::Bool=false; run_directory=nothing
     elseif input.discretization == "gausslegendre_pseudospectral"
         # create arrays needed for explicit GaussLegendre pseudospectral treatment in this
         # coordinate and create the matrices for differentiation
-        spectral = setup_gausslegendre_pseudospectral(coord, collision_operator_dim=collision_operator_dim,
-                                                      dirichlet_bc=occursin("zero", coord.bc))
+        spectral = setup_gausslegendre_pseudospectral(coord, collision_operator_dim=collision_operator_dim)
         # obtain the local derivatives of the uniform grid with respect to the used grid
         derivative!(coord.duniform_dgrid, coord.uniform_grid, coord, spectral)
     else
diff --git a/moment_kinetics/src/em_fields.jl b/moment_kinetics/src/em_fields.jl
index 6591e9062..67534d3c3 100644
--- a/moment_kinetics/src/em_fields.jl
+++ b/moment_kinetics/src/em_fields.jl
@@ -31,7 +31,7 @@ end
 """
 update_phi updates the electrostatic potential, phi
 """
-function update_phi!(fields, fvec, vperp, z, r, composition, z_spectral, r_spectral, scratch_dummy, gyroavs::gyro_operators)
+function update_phi!(fields, fvec, vperp, z, r, composition, geometry, z_spectral, r_spectral, scratch_dummy, gyroavs::gyro_operators)
     n_ion_species = composition.n_ion_species
     @boundscheck size(fields.phi,1) == z.n || throw(BoundsError(fields.phi))
     @boundscheck size(fields.phi,2) == r.n || throw(BoundsError(fields.phi))
@@ -125,8 +125,8 @@ function update_phi!(fields, fvec, vperp, z, r, composition, z_spectral, r_spect
         end
     else
         @loop_r_z ir iz begin
-            fields.Er[iz,ir] = composition.Er_constant
-            # Er_constant defaults to 0.0 in moment_kinetics_input.jl
+            fields.Er[iz,ir] = geometry.input.Er_constant
+            # Er_constant defaults to 0.0 in geo.jl
         end
     end
     #Ez = - d phi / dz
@@ -136,8 +136,9 @@ function update_phi!(fields, fvec, vperp, z, r, composition, z_spectral, r_spect
                 scratch_dummy.buffer_r_3, scratch_dummy.buffer_r_4,
                 z_spectral,z)
     else
-        @serial_region begin
-            fields.Ez[:,:] .= 0.0
+        @loop_r_z ir iz begin
+            fields.Ez[iz,ir] = geometry.input.Ez_constant
+            # Ez_constant defaults to 0.0 in geo.jl
         end
     end
     
diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl
index 1f520e9f6..fb9c6ef19 100644
--- a/moment_kinetics/src/fokker_planck.jl
+++ b/moment_kinetics/src/fokker_planck.jl
@@ -90,6 +90,7 @@ function setup_fkpl_collisions_input(toml_input::Dict, reference_params)
        nuii = -1.0,
        frequency_option = "reference_parameters",
        self_collisions = true,
+       use_conserving_corrections = true,
        slowing_down_test = false,
        sd_density = 1.0,
        sd_temp = 0.01,
@@ -257,6 +258,7 @@ function explicit_fp_collisions_weak_form_Maxwellian_cross_species!(pdf_out,pdf_
     @boundscheck r.n == size(dSdt,2) || throw(BoundsError(dSdt))
     @boundscheck n_ion_species == size(dSdt,3) || throw(BoundsError(dSdt))
     
+    use_conserving_corrections = collisions.fkpl.use_conserving_corrections
     fkin = collisions.fkpl
     # masses charge numbers and collision frequencies
     mref = 1.0 # generalise if multiple ions evolved
@@ -282,8 +284,10 @@ function explicit_fp_collisions_weak_form_Maxwellian_cross_species!(pdf_out,pdf_
         # enforce the boundary conditions on CC before it is used for timestepping
         enforce_vpavperp_BCs!(fkpl_arrays.CC,vpa,vperp,vpa_spectral,vperp_spectral)
         # make sure that the cross-species terms conserve density
-        density_conserving_correction!(fkpl_arrays.CC, pdf_in[:,:,iz,ir,is], vpa, vperp,
+        if use_conserving_corrections
+            density_conserving_correction!(fkpl_arrays.CC, pdf_in[:,:,iz,ir,is], vpa, vperp,
                                 fkpl_arrays.S_dummy)
+        end
         # advance this part of s,r,z with the resulting sum_s' C[Fs,Fs']
         begin_anyv_vperp_vpa_region()
         CC = fkpl_arrays.CC
@@ -326,6 +330,7 @@ function explicit_fokker_planck_collisions_weak_form!(pdf_out,pdf_in,dSdt,compos
     nuref = collisions.fkpl.nuii # generalise!
     Zi = collisions.fkpl.Zi # generalise!
     nussp = nuref*(Zi^4) # include charge number factor for self collisions
+    use_conserving_corrections = collisions.fkpl.use_conserving_corrections
     # N.B. parallelisation using special 'anyv' region
     begin_s_r_z_anyv_region()
     @loop_s_r_z is ir iz begin
@@ -336,9 +341,10 @@ function explicit_fokker_planck_collisions_weak_form!(pdf_out,pdf_in,dSdt,compos
         # enforce the boundary conditions on CC before it is used for timestepping
         enforce_vpavperp_BCs!(fkpl_arrays.CC,vpa,vperp,vpa_spectral,vperp_spectral)
         # make ad-hoc conserving corrections
-        conserving_corrections!(fkpl_arrays.CC, pdf_in[:,:,iz,ir,is], vpa, vperp,
+        if use_conserving_corrections
+            conserving_corrections!(fkpl_arrays.CC, pdf_in[:,:,iz,ir,is], vpa, vperp,
                                 fkpl_arrays.S_dummy)
-        
+        end
         # advance this part of s,r,z with the resulting C[Fs,Fs]
         begin_anyv_vperp_vpa_region()
         CC = fkpl_arrays.CC
diff --git a/moment_kinetics/src/fokker_planck_calculus.jl b/moment_kinetics/src/fokker_planck_calculus.jl
index 300e0ca6e..f11dcf43c 100644
--- a/moment_kinetics/src/fokker_planck_calculus.jl
+++ b/moment_kinetics/src/fokker_planck_calculus.jl
@@ -2287,21 +2287,29 @@ function enforce_vpavperp_BCs!(pdf,vpa,vperp,vpa_spectral,vperp_spectral)
     D0 = vperp_spectral.radau.D0
     # vpa boundary conditions
     # zero at infinity
-    begin_anyv_vperp_region()
-    @loop_vperp ivperp begin
-        pdf[1,ivperp] = 0.0
-        pdf[nvpa,ivperp] = 0.0
+    if vpa.bc == "zero"
+        begin_anyv_vperp_region()
+        @loop_vperp ivperp begin
+            pdf[1,ivperp] = 0.0
+            pdf[nvpa,ivperp] = 0.0
+        end
     end
     # vperp boundary conditions
     # zero boundary condition at infinity
     # set regularity condition d F / d vperp = 0 at vperp = 0
     # adjust F(vperp = 0) so that d F / d vperp = 0 at vperp = 0
     begin_anyv_vpa_region()
-    buffer = @view vperp.scratch[1:ngrid_vperp-1]
-    @loop_vpa ivpa begin
-        pdf[ivpa,nvperp] = 0.0
-        @views @. buffer = D0[2:ngrid_vperp] * pdf[ivpa,2:ngrid_vperp]
-        pdf[ivpa,1] = -sum(buffer)/D0[1]
+    if vperp.bc in ("zero", "zero-impose-regularity")
+        @loop_vpa ivpa begin
+            pdf[ivpa,nvperp] = 0.0
+        end
+    end
+    if vperp.bc == "zero-impose-regularity"
+        buffer = @view vperp.scratch[1:ngrid_vperp-1]
+        @loop_vpa ivpa begin
+            @views @. buffer = D0[2:ngrid_vperp] * pdf[ivpa,2:ngrid_vperp]
+            pdf[ivpa,1] = -sum(buffer)/D0[1]
+        end
     end
 end
 
diff --git a/moment_kinetics/src/gauss_legendre.jl b/moment_kinetics/src/gauss_legendre.jl
index 539a5fd86..9a0afb6b3 100644
--- a/moment_kinetics/src/gauss_legendre.jl
+++ b/moment_kinetics/src/gauss_legendre.jl
@@ -100,7 +100,7 @@ struct gausslegendre_info{TSparse, TLU} <: weak_discretization_info
     Qmat::Array{mk_float,2}
 end
 
-function setup_gausslegendre_pseudospectral(coord; collision_operator_dim=true, dirichlet_bc=true)
+function setup_gausslegendre_pseudospectral(coord; collision_operator_dim=true, dirichlet_bc=false)
     lobatto = setup_gausslegendre_pseudospectral_lobatto(coord,collision_operator_dim=collision_operator_dim)
     radau = setup_gausslegendre_pseudospectral_radau(coord,collision_operator_dim=collision_operator_dim)
 
@@ -887,9 +887,12 @@ function setup_global_weak_form_matrix!(QQ_global::Array{mk_float,2},
     if dirichlet_bc
         # Make matrix diagonal for first/last grid points so it does not change the values
         # there
-        if coord.irank == 0
-            QQ_global[1,:] .= 0.0
-            QQ_global[1,1] = 1.0
+        if !(coord.name == "vperp") 
+            # modify lower endpoint if not a radial/cylindrical coordinate
+            if coord.irank == 0
+                QQ_global[1,:] .= 0.0
+                QQ_global[1,1] = 1.0
+            end
         end
         if coord.irank == coord.nrank - 1
             QQ_global[end,:] .= 0.0
diff --git a/moment_kinetics/src/geo.jl b/moment_kinetics/src/geo.jl
index 9ad6d4885..3635230dc 100644
--- a/moment_kinetics/src/geo.jl
+++ b/moment_kinetics/src/geo.jl
@@ -8,7 +8,7 @@ module geo
 export init_magnetic_geometry
 export setup_geometry_input
 
-using ..input_structs: geometry_input
+using ..input_structs: geometry_input, set_defaults_and_check_section!
 using ..file_io: input_option_error
 using ..array_allocation: allocate_float
 using ..type_definitions: mk_float, mk_int
@@ -73,11 +73,26 @@ option = ""
 
 """
 function setup_geometry_input(toml_input::Dict, reference_rhostar)
-    input_section = get(toml_input, "geometry", Dict{String,Any}())
-    if !("rhostar" ∈ keys(input_section))
-        # Set default rhostar with reference value
-        input_section["rhostar"] = get(input_section, "rhostar", reference_rhostar)
-    end
+    # read the input toml and specify a sensible default
+    input_section = set_defaults_and_check_section!(toml_input, "geometry",
+        # begin default inputs (as kwargs)
+        # rhostar ion (ref)
+        rhostar = reference_rhostar, #used to premultiply ExB drift terms
+        # magnetic geometry option
+        option = "constant-helical",# "1D-mirror"
+        # pitch ( = Bzed/Bmag if geometry_option == "constant-helical")
+        pitch = 1.0,
+        # DeltaB ( = (Bzed(z=L/2) - Bzed(0))/Bref if geometry_option == "1D-mirror")
+        DeltaB = 0.0,
+        # constant for testing nonzero Er when nr = 1
+        Er_constant = 0.0,
+        # constant for testing nonzero Ez when nz = 1
+        Ez_constant = 0.0,
+        # constant for testing nonzero dBdz when nz = 1
+        dBdz_constant = 0.0,
+        # constant for testing nonzero dBdr when nr = 1
+        dBdr_constant = 0.0)
+    
     input = Dict(Symbol(k)=>v for (k,v) in input_section)
     #println(input)
     return geometry_input(; input...)
@@ -187,6 +202,31 @@ function init_magnetic_geometry(geometry_input_data::geometry_input,z,r)
                 gbdriftz[iz,ir] = cvdriftz[iz,ir]
             end
         end
+    elseif option == "0D-Spitzer-test"
+     # a 0D configuration with certain geometrical factors
+     # set to be constants to enable testing of velocity
+     # space operators such as mirror or vperp advection terms
+     pitch = geometry_input_data.pitch
+     dBdz_constant = geometry_input_data.dBdz_constant
+     dBdr_constant = geometry_input_data.dBdr_constant
+     B0 = 1.0 # chose reference field strength to be Bzeta at r = 1
+     for ir in 1:nr
+            for iz in 1:nz
+                Bmag[iz,ir] = B0
+                bzed[iz,ir] = pitch
+                bzeta[iz,ir] = sqrt(1 - pitch^2)
+                Bzed[iz,ir] = bzed[iz,ir]*Bmag[iz,ir]
+                Bzeta[iz,ir] = bzeta[iz,ir]*Bmag[iz,ir]
+                dBdz[iz,ir] = dBdz_constant
+                dBdr[iz,ir] = dBdr_constant
+                jacobian[iz,ir] = 1.0
+                
+                cvdriftr[iz,ir] = 0.0
+                cvdriftz[iz,ir] = 0.0
+                gbdriftr[iz,ir] = 0.0
+                gbdriftz[iz,ir] = 0.0
+            end
+        end
     else 
         input_option_error("$option", option)
     end
diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl
index f81d28a30..b740615c9 100644
--- a/moment_kinetics/src/input_structs.jl
+++ b/moment_kinetics/src/input_structs.jl
@@ -318,7 +318,6 @@ mutable struct species_composition
     # wall potential used if electron_physics=boltzmann_electron_response_with_simple_sheath
     phi_wall::mk_float
     # constant for testing nonzero Er
-    Er_constant::mk_float
     # ratio of the neutral particle mass to the ion mass
     mn_over_mi::mk_float
     # ratio of the electron particle mass to the ion mass
@@ -378,6 +377,8 @@ Base.@kwdef struct fkpl_collisions_input
     nuii::mk_float
     # option to determine if self collisions are used (for physics test)
     self_collisions::Bool
+    # option to determine if ad-hoc moment_kinetics-style conserving corrections are used
+    use_conserving_corrections::Bool
     # option to determine if cross-collisions against fixed Maxwellians are used
     slowing_down_test::Bool
     # Setting to switch between different options for Fokker-Planck collision frequency input
@@ -425,6 +426,14 @@ Base.@kwdef struct geometry_input
     pitch::mk_float = 1.0
     # DeltaB ( = (Bzed(z=L/2) - Bzed(0))/Bref if geometry_option == "1D-mirror")
     DeltaB::mk_float = 0.0
+    # constant for testing nonzero Er when nr = 1
+    Er_constant::mk_float
+    # constant for testing nonzero Ez when nz = 1
+    Ez_constant::mk_float
+    # constant for testing nonzero dBdz when nz = 1
+    dBdz_constant::mk_float
+    # constant for testing nonzero dBdr when nr = 1
+    dBdr_constant::mk_float
 end
 
 @enum binary_format_type hdf5 netcdf
@@ -556,6 +565,8 @@ struct pp_input
 	animate_parallel_flow_vs_r_z::Bool
 	# if plot_parallel_pressure_vs_r0_z = true  plot last timestep parallel_pressure[z,ir0]
 	plot_parallel_pressure_vs_r0_z::Bool
+	# if plot_perpendicular_pressure_vs_r0_z = true  plot last timestep perpendicular_pressure[z,ir0]
+	plot_perpendicular_pressure_vs_r0_z::Bool
 	# if plot_wall_parallel_pressure_vs_r = true  plot last timestep parallel_pressure[z_wall,r]
 	plot_wall_parallel_pressure_vs_r::Bool
     # if plot_parallel_pressure_vs_r_z = true plot parallel_pressure vs r z at last timestep 
diff --git a/moment_kinetics/src/moment_kinetics_input.jl b/moment_kinetics/src/moment_kinetics_input.jl
index 82b8b00fb..3bdba79f8 100644
--- a/moment_kinetics/src/moment_kinetics_input.jl
+++ b/moment_kinetics/src/moment_kinetics_input.jl
@@ -103,7 +103,7 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
     # if false use true Knudsen cosine for neutral wall bc
     composition.use_test_neutral_wall_pdf = get(scan_input, "use_test_neutral_wall_pdf", false)
     # constant to be used to test nonzero Er in wall boundary condition
-    composition.Er_constant = get(scan_input, "Er_constant", 0.0)
+    #composition.Er_constant = get(scan_input, "Er_constant", 0.0)
     # The ion flux reaching the wall that is recycled as neutrals is reduced by
     # `recycling_fraction` to account for ions absorbed by the wall.
     composition.recycling_fraction = get(scan_input, "recycling_fraction", 1.0)
@@ -550,7 +550,7 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
             species.ion[is].vpa_IC.upar_amplitude, species.ion[is].vpa_IC.upar_phase,
             species.ion[is].vpa_IC.temperature_amplitude,
             species.ion[is].vpa_IC.temperature_phase, species.ion[is].vpa_IC.monomial_degree,
-            get(scan_input, "vpa_IC_v0_$is", 0.5*sqrt(vperp.L^2 + (0.5*vpa.L)^2)),
+            get(scan_input, "vpa_IC_v0$is", 0.5*sqrt(vperp.L^2 + (0.5*vpa.L)^2)),
             get(scan_input, "vpa_IC_vth0$is", 0.1*sqrt(vperp.L^2 + (0.5*vpa.L)^2)),
             get(scan_input, "vpa_IC_vpa0$is", 0.25*0.5*abs(vpa.L)),
             get(scan_input, "vpa_IC_vperp0$is", 0.5*abs(vperp.L)))
@@ -1016,8 +1016,6 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics)
     T_wall = 1.0
     # wall potential at z = 0
     phi_wall = 0.0
-    # constant to test nonzero Er
-    Er_constant = 0.0
     # ratio of the neutral particle mass to the ion particle mass
     mn_over_mi = 1.0
     # ratio of the electron particle mass to the ion particle mass
@@ -1027,7 +1025,7 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics)
     recycling_fraction = 1.0
     gyrokinetic_ions = false
     composition = species_composition(n_species, n_ion_species, n_neutral_species,
-        electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, Er_constant,
+        electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall,
         mn_over_mi, me_over_mi, recycling_fraction, gyrokinetic_ions, allocate_float(n_species))
     
     species_ion = Array{species_parameters_mutable,1}(undef,n_ion_species)
@@ -1198,6 +1196,8 @@ function check_coordinate_input(coord, coord_name, io)
         println(io,">$coord_name.bc = 'constant'.  enforcing constant incoming BC in $coord_name.")
     elseif coord.bc == "zero"
         println(io,">$coord_name.bc = 'zero'.  enforcing zero incoming BC in $coord_name. Enforcing zero at both boundaries if diffusion operator is present.")
+    elseif coord.bc == "zero-impose-regularity"
+        println(io,">$coord_name.bc = 'zero'.  enforcing zero incoming BC in $coord_name. Enforcing zero at both boundaries if diffusion operator is present. Enforce dF/dcoord = 0 at origin if coord = vperp.")
     elseif coord.bc == "zero_gradient"
         println(io,">$coord_name.bc = 'zero_gradient'.  enforcing zero gradients at both limits of $coord_name domain.")
     elseif coord.bc == "both_zero"
@@ -1216,9 +1216,9 @@ function check_coordinate_input(coord, coord_name, io)
                 coord.nelement_global, " elements across the $coord_name domain [",
                 0.0, ",", coord.L, "].")
 
-        if coord.bc != "zero" && coord.n_global > 1 && global_rank[] == 0
-            println("WARNING: regularity condition (df/dvperp=0 at vperp=0) not being "
-                    * "imposed. Collisions or vperp-diffusion will be unstable.")
+        if coord.bc == "zero-impose-regularity" && coord.n_global > 1 && global_rank[] == 0
+            println("WARNING: regularity condition (df/dvperp=0 at vperp=0) being "
+                    * "imposed explicitly.")
         end
     else
         println(io,">using ", coord.ngrid, " grid points per $coord_name element on ",
diff --git a/moment_kinetics/src/numerical_dissipation.jl b/moment_kinetics/src/numerical_dissipation.jl
index 86e517229..f48a158c6 100644
--- a/moment_kinetics/src/numerical_dissipation.jl
+++ b/moment_kinetics/src/numerical_dissipation.jl
@@ -12,7 +12,7 @@ export setup_numerical_dissipation, vpa_boundary_buffer_decay!,
 using Base.Iterators: flatten
 
 using ..looping
-using ..calculus: derivative!, second_derivative!
+using ..calculus: derivative!, second_derivative!, laplacian_derivative!
 using ..derivatives: derivative_r!, derivative_z!, second_derivative_r!,
                      second_derivative_z!
 using ..type_definitions: mk_float
@@ -369,7 +369,7 @@ function vperp_dissipation!(f_out, f_in, vperp, spectral::T_spectral, dt,
     end
     
     @loop_s_r_z_vpa is ir iz ivpa begin
-        @views second_derivative!(vperp.scratch, f_in[ivpa,:,iz,ir,is], vperp, spectral)
+        @views laplacian_derivative!(vperp.scratch, f_in[ivpa,:,iz,ir,is], vperp, spectral)
         @views @. f_out[ivpa,:,iz,ir,is] += dt * diffusion_coefficient * vperp.scratch
     end
 
diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl
index 2460dba15..6bfad2755 100644
--- a/moment_kinetics/src/time_advance.jl
+++ b/moment_kinetics/src/time_advance.jl
@@ -529,7 +529,7 @@ function setup_time_advance!(pdf, fields, vz, vr, vzeta, vpa, vperp, z, r, gyrop
 
     # initialize the electrostatic potential
     begin_serial_region()
-    update_phi!(fields, scratch[1], vperp, z, r, composition, z_spectral, r_spectral, scratch_dummy, gyroavs)
+    update_phi!(fields, scratch[1], vperp, z, r, composition, geometry, z_spectral, r_spectral, scratch_dummy, gyroavs)
     @serial_region begin
         # save the initial phi(z) for possible use later (e.g., if forcing phi)
         fields.phi0 .= fields.phi
@@ -727,7 +727,7 @@ function setup_time_advance!(pdf, fields, vz, vr, vzeta, vpa, vperp, z, r, gyrop
         end
     end
 
-    update_phi!(fields, scratch[1], vperp, z, r, composition, z_spectral, r_spectral,
+    update_phi!(fields, scratch[1], vperp, z, r, composition, geometry, z_spectral, r_spectral,
                 scratch_dummy, gyroavs)
     calculate_ion_moment_derivatives!(moments, scratch[1], scratch_dummy, z, z_spectral, 
                                       ion_mom_diss_coeff)
@@ -792,8 +792,8 @@ function setup_advance_flags(moments, composition, t_params, collisions,
         # default for non-split operators is to include both vpa and z advection together
         # If using an IMEX scheme and implicit vpa advection has been requested, then vpa
         # advection is not included in the explicit part of the timestep.
-        advance_vpa_advection = vpa.n > 1 && z.n > 1 && !(t_params.implicit_ion_advance || t_params.implicit_vpa_advection)
-        advance_vperp_advection = vperp.n > 1 && z.n > 1 && !t_params.implicit_ion_advance
+        advance_vpa_advection = vpa.n > 1 && !(t_params.implicit_ion_advance || t_params.implicit_vpa_advection)
+        advance_vperp_advection = vperp.n > 1 && !t_params.implicit_ion_advance
         advance_z_advection = z.n > 1 && !t_params.implicit_ion_advance
         advance_r_advection = r.n > 1 && !t_params.implicit_ion_advance
         if collisions.fkpl.nuii > 0.0 && vperp.n > 1 && !t_params.implicit_ion_advance
@@ -1871,7 +1871,7 @@ function apply_all_bcs_constraints_update_moments!(
                                       num_diss_params.ion.moment_dissipation_coefficient)
 
     # update the electrostatic potential phi
-    update_phi!(fields, this_scratch, vperp, z, r, composition, z_spectral, r_spectral,
+    update_phi!(fields, this_scratch, vperp, z, r, composition, geometry, z_spectral, r_spectral,
                 scratch_dummy, gyroavs)
 
     if composition.n_neutral_species > 0
@@ -2476,7 +2476,7 @@ function euler_time_advance!(fvec_out, fvec_in, pdf, fields, moments,
     end
 
     # r advection relies on derivatives in z to get ExB
-    if advance.r_advection && r.n > 1
+    if advance.r_advection
         r_advection!(fvec_out.pdf, fvec_in, moments, fields, r_advect, r, z, vperp, vpa,
                      dt, r_spectral, composition, geometry, scratch_dummy)
     end
@@ -2484,7 +2484,8 @@ function euler_time_advance!(fvec_out, fvec_in, pdf, fields, moments,
     # so call vperp_advection! only after z and r advection routines
     if advance.vperp_advection
         vperp_advection!(fvec_out.pdf, fvec_in, vperp_advect, r, z, vperp, vpa,
-                      dt, vperp_spectral, composition, z_advect, r_advect, geometry)
+                      dt, vperp_spectral, composition, z_advect, r_advect, geometry,
+                      moments, fields, t)
     end
 
     if advance.source_terms
@@ -2497,7 +2498,7 @@ function euler_time_advance!(fvec_out, fvec_in, pdf, fields, moments,
             r, z, vzeta, vr, vz, dt, t, z_spectral, composition, scratch_dummy)
     end
 
-    if advance.neutral_r_advection && r.n > 1
+    if advance.neutral_r_advection
         neutral_advection_r!(fvec_out.pdf_neutral, fvec_in, neutral_r_advect,
             r, z, vzeta, vr, vz, dt, r_spectral, composition, geometry, scratch_dummy)
     end
diff --git a/moment_kinetics/src/vperp_advection.jl b/moment_kinetics/src/vperp_advection.jl
index 16c2697df..2dc6bd5c8 100644
--- a/moment_kinetics/src/vperp_advection.jl
+++ b/moment_kinetics/src/vperp_advection.jl
@@ -6,13 +6,21 @@ export update_speed_vperp!
 using ..advection: advance_f_local!
 using ..chebyshev: chebyshev_info
 using ..looping
+using ..z_advection: update_speed_z!
+using ..r_advection: update_speed_r!
 
 # do a single stage time advance (potentially as part of a multi-stage RK scheme)
 function vperp_advection!(f_out, fvec_in, vperp_advect, r, z, vperp, vpa,
-                      dt, vperp_spectral, composition, z_advect, r_advect, geometry)
+                      dt, vperp_spectral, composition, z_advect, r_advect, geometry,
+                      moments, fields, t)
+    
+    # if appropriate, update z and r speeds
+    update_z_r_speeds!(z_advect, r_advect, fvec_in, moments, fields,
+                            geometry, vpa, vperp, z, r, t)
+    
     begin_s_r_z_vpa_region()
     @loop_s is begin
-        # get the updated speed along the r direction using the current f
+        # get the updated speed along the vperp direction using the current f
         @views update_speed_vperp!(vperp_advect[is], vpa, vperp, z, r, z_advect[is], r_advect[is], geometry)
         @loop_r_z_vpa ir iz ivpa begin
             @views advance_f_local!(f_out[ivpa,:,iz,ir,is], fvec_in.pdf[ivpa,:,iz,ir,is],
@@ -58,4 +66,34 @@ function update_speed_vperp!(vperp_advect, vpa, vperp, z, r, z_advect, r_advect,
     return nothing
 end
 
+function update_z_r_speeds!(z_advect, r_advect, fvec_in, moments, fields,
+                            geometry, vpa, vperp, z, r, t)
+    update_z_speed = (z.n == 1 && geometry.input.option == "0D-Spitzer-test")
+    if update_z_speed
+    # make sure z speed is physical despite
+    # 0D nature of simulation
+        begin_s_r_vperp_vpa_region()
+        @loop_s is begin
+            # get the updated speed along the z direction using the current f
+            @views update_speed_z!(z_advect[is], fvec_in.upar[:,:,is],
+                                   moments.ion.vth[:,:,is], moments.evolve_upar,
+                                   moments.evolve_ppar, fields, vpa, vperp, z, r, t, geometry, is)
+        end
+    end
+    
+    update_r_speed = (r.n == 1 && geometry.input.option == "0D-Spitzer-test")
+    if update_r_speed
+    # make sure r speed is physical despite
+    # 0D nature of simulation
+        begin_s_z_vperp_vpa_region()
+        @loop_s is begin
+            # get the updated speed along the r direction using the current f
+            update_speed_r!(r_advect[is], fvec_in.upar[:,:,is],
+                               moments.ion.vth[:,:,is], fields, moments.evolve_upar,
+                               moments.evolve_ppar, vpa, vperp, z, r, geometry, is)
+        end
+    end
+    return nothing
+end    
+
 end
diff --git a/moment_kinetics/test/calculus_tests.jl b/moment_kinetics/test/calculus_tests.jl
index e439ac4cf..ef13b75bb 100644
--- a/moment_kinetics/test/calculus_tests.jl
+++ b/moment_kinetics/test/calculus_tests.jl
@@ -5,6 +5,7 @@ include("setup.jl")
 using moment_kinetics.input_structs: grid_input, advection_input
 using moment_kinetics.coordinates: define_coordinate
 using moment_kinetics.calculus: derivative!, second_derivative!, integral
+using moment_kinetics.calculus: laplacian_derivative!
 
 using MPI
 using Random
@@ -1355,6 +1356,101 @@ function runtests()
             end
         end
         
+        @testset "Chebyshev pseudospectral cylindrical laplacian derivatives (4 argument), zero" verbose=false begin
+            @testset "$nelement $ngrid" for (nelement, ngrid, rtol) ∈
+                    (
+                     (4, 7, 2.e-1),
+                     (4, 8, 2.e-1),
+                     (4, 9, 4.e-2),
+                     (4, 10, 4.e-2),
+                     (4, 11, 5.e-3),
+                     (4, 12, 5.e-3),
+                     (4, 13, 5.e-3),
+                     (4, 14, 5.e-3),
+                     (4, 15, 5.e-3),
+                     (4, 16, 5.e-3),
+                     (4, 17, 5.e-3),
+                     (4, 18, 5.e-3),
+                     (4, 19, 5.e-3),
+                     (4, 20, 5.e-3),
+                     (4, 21, 5.e-3),
+                     (4, 22, 5.e-3),
+                     (4, 23, 5.e-3),
+                     (4, 24, 4.e-3),
+                     (4, 25, 4.e-3),
+                     (4, 26, 4.e-3),
+                     (4, 27, 4.e-3),
+                     (4, 28, 4.e-3),
+                     (4, 29, 4.e-3),
+                     (4, 30, 4.e-3),
+                     (4, 31, 4.e-3),
+                     (4, 32, 4.e-3),
+                     (4, 33, 4.e-3),
+
+                     (5, 7, 2.e-1),
+                     (5, 8, 8.e-2),
+                     (5, 9, 5.e-2),
+                     (5, 10, 8.e-3),
+                     (5, 11, 8.e-3),
+                     (5, 12, 8.e-3),
+                     (5, 13, 8.e-3),
+                     (5, 14, 8.e-3),
+                     (5, 15, 8.e-3),
+                     (5, 16, 2.e-3),
+                     (5, 17, 2.e-3),
+                     (5, 18, 2.e-3),
+                     (5, 19, 2.e-3),
+                     (5, 20, 2.e-3),
+                     (5, 21, 2.e-3),
+                     (5, 22, 2.e-3),
+                     (5, 23, 2.e-3),
+                     (5, 24, 2.e-3),
+                     (5, 25, 2.e-3),
+                     (5, 26, 2.e-3),
+                     (5, 27, 2.e-3),
+                     (5, 28, 2.e-3),
+                     (5, 29, 2.e-3),
+                     (5, 30, 2.e-3),
+                     (5, 31, 2.e-3),
+                     (5, 32, 2.e-3),
+                     (5, 33, 2.e-3),
+                    ), cheb_option in ("FFT","matrix")
+
+                # define inputs needed for the test
+                L = 6.0
+                bc = "zero"
+                # fd_option and adv_input not actually used so given values unimportant
+                fd_option = ""
+                adv_input = advection_input("default", 1.0, 0.0, 0.0)
+                # create the 'input' struct containing input info needed to create a
+                # coordinate
+                nelement_local = nelement
+				nrank_per_block = 0 # dummy value
+				irank = 0 # dummy value
+				comm = MPI.COMM_NULL # dummy value
+                element_spacing_option = "uniform"
+				input = grid_input("vperp", ngrid, nelement,
+                    nelement_local, nrank_per_block, irank, L,
+                    "chebyshev_pseudospectral", fd_option, cheb_option, bc, adv_input, comm,
+                    element_spacing_option)
+                # create the coordinate struct 'x' and info for derivatives, etc.
+                # This test runs effectively in serial, so use `ignore_MPI=true` to avoid
+                # errors due to communicators not being fully set up.
+                x, spectral = define_coordinate(input; ignore_MPI=true)
+
+                f = @. exp(-x.grid^2)
+                expected_d2f = @. 4.0*(x.grid^2 - 1.0)*exp(-x.grid^2)
+                # create array for the derivative d2f/dx2
+                d2f = similar(f)
+
+                # differentiate f
+                laplacian_derivative!(d2f, f, x, spectral)
+
+                @test isapprox(d2f, expected_d2f, rtol=rtol, atol=1.e-10,
+                               norm=maxabs_norm)
+            end
+        end
+        
         @testset "GaussLegendre pseudospectral second derivatives (4 argument), periodic" verbose=false begin
             @testset "$nelement $ngrid" for (nelement, ngrid, rtol) ∈
                     (
@@ -1460,6 +1556,111 @@ function runtests()
                 # differentiate f
                 second_derivative!(d2f, f, x, spectral)
 
+                @test isapprox(d2f, expected_d2f, rtol=rtol, atol=1.e-10,
+                               norm=maxabs_norm)
+            end
+        end
+        
+    @testset "GaussLegendre pseudospectral cylindrical laplacian derivatives (4 argument), zero" verbose=false begin
+            @testset "$nelement $ngrid" for (nelement, ngrid, rtol) ∈
+                    (
+                     (1, 8, 5.e-2),
+                     (1, 9, 1.e-1),
+                     (1, 10, 2.e-1),
+                     (1, 11, 5.e-2),
+                     (1, 12, 5.e-2),
+                     (1, 13, 5.e-2),
+                     (1, 14, 5.e-2),
+                     (1, 15, 5.e-3),
+                     (1, 16, 5.e-2),
+                     (1, 17, 5.e-3),
+                     
+                     (2, 6, 8.e-2),
+                     (2, 7, 8.e-2),
+                     (2, 8, 5.e-2),
+                     (2, 9, 5.e-2),
+                     (2, 10, 5.e-2),
+                     (2, 11, 5.e-3),
+                     (2, 12, 5.e-3),
+                     (2, 13, 5.e-4),
+                     (2, 14, 5.e-4),
+                     (2, 15, 5.e-4),
+                     (2, 16, 5.e-4),
+                     (2, 17, 5.e-4),
+                     
+                     (3, 6, 5.e-2),
+                     (3, 7, 5.e-3),
+                     (3, 8, 5.e-2),
+                     (3, 9, 5.e-4),
+                     (3, 10, 5.e-3),
+                     (3, 11, 5.e-4),
+                     (3, 12, 5.e-5),
+                     (3, 13, 5.e-5),
+                     (3, 14, 5.e-6),
+                     (3, 15, 5.e-6),
+                     (3, 16, 5.e-6),
+                     (3, 17, 5.e-8),
+                     
+                     (4, 5, 5.e-2),
+                     (4, 6, 2.e-2),
+                     (4, 7, 2.e-2),
+                     (4, 8, 2.e-3),
+                     (4, 9, 1.e-3),
+                     (4, 10, 1.e-4),
+                     (4, 11, 8.e-5),
+                     (4, 12, 8.e-6),
+                     (4, 13, 8.e-6),
+                     (4, 14, 8.e-7),
+                     (4, 15, 8.e-7),
+                     (4, 16, 8.e-8),
+                     (4, 17, 8.e-8),
+                     
+                     (5, 5, 4.e-2),
+                     (5, 6, 8.e-3),
+                     (5, 7, 5.e-3),
+                     (5, 8, 5.e-4),
+                     (5, 9, 8.e-5),
+                     (5, 10, 5.e-6),
+                     (5, 11, 8.e-6),
+                     (5, 12, 4.e-7),
+                     (5, 13, 2.e-7),
+                     (5, 14, 2.e-7),
+                     (5, 15, 8.e-7),
+                     (5, 16, 8.e-10),
+                     (5, 17, 8.e-10),
+                     )
+
+                # define inputs needed for the test
+                L = 6.0
+                bc = "zero"
+                # fd_option and adv_input not actually used so given values unimportant
+                fd_option = ""
+                adv_input = advection_input("default", 1.0, 0.0, 0.0)
+                cheb_option = ""
+                # create the 'input' struct containing input info needed to create a
+                # coordinate
+                nelement_local = nelement
+				nrank_per_block = 0 # dummy value
+				irank = 0 # dummy value
+				comm = MPI.COMM_NULL # dummy value
+                element_spacing_option = "uniform"
+				input = grid_input("vperp", ngrid, nelement,
+                    nelement_local, nrank_per_block, irank, L,
+                    "gausslegendre_pseudospectral", fd_option, cheb_option, bc, adv_input, comm,
+                    element_spacing_option)
+                # create the coordinate struct 'x' and info for derivatives, etc.
+                # This test runs effectively in serial, so use `ignore_MPI=true` to avoid
+                # errors due to communicators not being fully set up.
+                x, spectral = define_coordinate(input; ignore_MPI=true, collision_operator_dim=false)
+
+                f = @. exp(-x.grid^2)
+                expected_d2f = @. 4.0*(x.grid^2 - 1.0)*exp(-x.grid^2)
+                # create array for the derivative d2f/dx2
+                d2f = similar(f)
+
+                # differentiate f
+                laplacian_derivative!(d2f, f, x, spectral)
+
                 @test isapprox(d2f, expected_d2f, rtol=rtol, atol=1.e-10,
                                norm=maxabs_norm)
             end
diff --git a/moment_kinetics/test/fokker_planck_time_evolution_tests.jl b/moment_kinetics/test/fokker_planck_time_evolution_tests.jl
index 08239c0f7..a02fddf70 100644
--- a/moment_kinetics/test/fokker_planck_time_evolution_tests.jl
+++ b/moment_kinetics/test/fokker_planck_time_evolution_tests.jl
@@ -34,7 +34,7 @@ struct expected_data
     f_ion::Array{mk_float, 3} # vpa, vperp, time
 end
 
-const expected =
+const expected_zero_impose_regularity =
   expected_data(
    [-3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0],
    [0.155051025721682, 0.644948974278318, 1.000000000000000, 1.500000000000000, 2.000000000000000, 2.500000000000000, 3.000000000000000],
@@ -81,34 +81,145 @@ const expected =
     0.005877384425022965 0.00466291282312835 0.0028530945920350282 0.0008130106752860425 2.4399432485774198e-5 -9.743480904212476e-5 0.0;
     0.00017108987342944414 7.097261590252536e-5 -7.822316408657793e-5 -0.0001534897546375058 -9.087984332447822e-5 -3.3079379573126077e-5 0.0;
     0.0 0.0 0.0 0.0 0.0 0.0 0.0])
+const expected_zero = 
+expected_data(
+   [-3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0],
+   [0.155051025721682, 0.644948974278318, 1.000000000000000, 1.500000000000000, 2.000000000000000, 2.500000000000000, 3.000000000000000],
+   # Expected phi:
+   [-1.254259088243025, -1.254259088243286],
+   # Expected n_ion:
+   [0.285287142537587, 0.285287142537513],
+   # Expected upar_ion:
+   [0.0, 0.0],
+   # Expected ppar_ion:
+   [0.182220654804438, 0.156448883483764],
+   # Expected pperp_ion
+   [0.143306715174515, 0.156192600834786],
+   # Expected qpar_ion
+   [0.0, 0.0],
+   # Expected v_t_ion
+   [1.046701532502699, 1.046701532502689],
+   # Expected dSdt
+   [0.000000000000000, -0.000000000865997],
+   # Expected f_ion:
+   [0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 ;
+    0.000706724195475 0.000477575434536 0.000266315395816 0.000076300638379 0.000013259062819 0.000001397494940 0.000000000000000 ;
+    0.006705212475828 0.004531109564814 0.002526730124657 0.000723920301085 0.000125798485463 0.000013259062819 0.000000000000000 ;
+    0.038585833650058 0.026074735222541 0.014540327934436 0.004165873701136 0.000723920300963 0.000076300638366 0.000000000000000 ;
+    0.134671678398729 0.091005636629981 0.050748427134097 0.014539667807029 0.002526615411287 0.000266303305116 0.000000000000000 ;
+    0.261709398369233 0.176852554997681 0.098620144126568 0.028255144359305 0.004910007858078 0.000517511020834 0.000000000000000 ;
+    0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 ;
+    0.261709398369233 0.176852554997682 0.098620144126568 0.028255144359305 0.004910007858078 0.000517511020834 0.000000000000000 ;
+    0.134671678398729 0.091005636629981 0.050748427134097 0.014539667807029 0.002526615411287 0.000266303305116 0.000000000000000 ;
+    0.038585833650058 0.026074735222541 0.014540327934436 0.004165873701136 0.000723920300963 0.000076300638366 0.000000000000000 ;
+    0.006705212475828 0.004531109564814 0.002526730124657 0.000723920301085 0.000125798485463 0.000013259062819 0.000000000000000 ;
+    0.000706724195475 0.000477575434536 0.000266315395816 0.000076300638379 0.000013259062819 0.000001397494940 0.000000000000000 ;
+    0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 ;;;
+    0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 ;
+    0.000473874379555 0.000190726186343 -0.000067413540342 -0.000219129923463 -0.000165609965457 -0.000069940808382 0.000000000000000 ;
+    0.007980077110978 0.005701898246888 0.003327438535529 0.000925955494251 -0.000037963012435 -0.000174809935222 0.000000000000000 ;
+    0.034010103035913 0.024541441396039 0.014741892515397 0.004975785949797 0.000869171347288 -0.000277723510529 0.000000000000000 ;
+    0.095320265225675 0.069160937049138 0.041953067849510 0.014711859721633 0.003235165444330 -0.000199356792174 0.000000000000000 ;
+    0.176443375955397 0.128736406466194 0.078839749105573 0.028150317244036 0.006528121645593 0.000119350436884 0.000000000000000 ;
+    0.215552595150568 0.157471466895046 0.096656199821596 0.034660213647059 0.008120924795998 0.000303073103709 0.000000000000000 ;
+    0.176443375955397 0.128736406466194 0.078839749105573 0.028150317244036 0.006528121645593 0.000119350436884 0.000000000000000 ;
+    0.095320265225675 0.069160937049138 0.041953067849510 0.014711859721633 0.003235165444330 -0.000199356792174 0.000000000000000 ;
+    0.034010103035913 0.024541441396039 0.014741892515397 0.004975785949797 0.000869171347288 -0.000277723510529 0.000000000000000 ;
+    0.007980077110978 0.005701898246888 0.003327438535529 0.000925955494251 -0.000037963012435 -0.000174809935222 0.000000000000000 ;
+    0.000473874379555 0.000190726186343 -0.000067413540342 -0.000219129923463 -0.000165609965457 -0.000069940808382 0.000000000000000 ;
+    0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000])
+
+const expected_none_bc = 
+expected_data(
+   [-3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0],
+   [0.155051025721682, 0.644948974278318, 1.000000000000000, 1.500000000000000, 2.000000000000000, 2.500000000000000, 3.000000000000000],
+   # Expected phi:
+   [-1.254104813718096, -1.254104813718360],
+   # Expected n_ion:
+   [0.285331158471152, 0.285331158471077],
+   # Expected upar_ion:
+   [0.0, 0.0],
+   # Expected ppar_ion:
+   [0.182321263834348, 0.154595996906667],
+   # Expected pperp_ion
+   [0.143470130069393, 0.157332763533171],
+   # Expected qpar_ion
+   [0.0, 0.0],
+   # Expected v_t_ion
+   [1.047097792428007, 1.047097792428005],
+   # Expected dSdt
+   [0.000000000000000, 0.000000019115425],
+   # Expected f_ion:
+   [0.000045179366280 0.000030530376095 0.000017024973661 0.000004877736620 0.000000847623528 0.000000089338863 0.000000005711242 ;
+    0.000706724195475 0.000477575434536 0.000266315395816 0.000076300638379 0.000013259062819 0.000001397494940 0.000000089338863 ;
+    0.006705212475828 0.004531109564814 0.002526730124657 0.000723920301085 0.000125798485463 0.000013259062819 0.000000847623528 ;
+    0.038585833650058 0.026074735222541 0.014540327934436 0.004165873701136 0.000723920300963 0.000076300638366 0.000004877736619 ;
+    0.134671678398729 0.091005636629981 0.050748427134097 0.014539667807029 0.002526615411287 0.000266303305116 0.000017024200728 ;
+    0.261709398369233 0.176852554997681 0.098620144126568 0.028255144359305 0.004910007858078 0.000517511020834 0.000033083372713 ;
+    0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 ;
+    0.261709398369233 0.176852554997682 0.098620144126568 0.028255144359305 0.004910007858078 0.000517511020834 0.000033083372713 ;
+    0.134671678398729 0.091005636629981 0.050748427134097 0.014539667807029 0.002526615411287 0.000266303305116 0.000017024200728 ;
+    0.038585833650058 0.026074735222541 0.014540327934436 0.004165873701136 0.000723920300963 0.000076300638366 0.000004877736619 ;
+    0.006705212475828 0.004531109564814 0.002526730124657 0.000723920301085 0.000125798485463 0.000013259062819 0.000000847623528 ;
+    0.000706724195475 0.000477575434536 0.000266315395816 0.000076300638379 0.000013259062819 0.000001397494940 0.000000089338863 ;
+    0.000045179366280 0.000030530376095 0.000017024973661 0.000004877736620 0.000000847623528 0.000000089338863 0.000000005711242 ;;;
+    0.000447615535468 0.000364801746561 0.000270093752229 0.000138454732861 0.000052917434226 0.000014331973588 -0.000000022222325 ;
+    0.000530676740860 0.000337644325418 0.000161654162384 0.000044642657051 0.000030848860156 0.000021605350900 0.000014496785601 ;
+    0.006652401806692 0.004725322877613 0.002773765135978 0.000908975538425 0.000210924069594 0.000037798409136 0.000057517483301 ;
+    0.030674916797523 0.021423584760802 0.012211671337997 0.003838328661082 0.000895825759878 0.000045643481270 0.000159762122747 ;
+    0.097751022951371 0.068620792366922 0.039365625553724 0.012397445282748 0.002818297719268 0.000149127033057 0.000311439525732 ;
+    0.194564361178625 0.137570544970062 0.079900546066037 0.025549350260855 0.005673491663727 0.000408759661293 0.000445141805523 ;
+    0.242684098382075 0.171867841736931 0.100113917461427 0.032148650920148 0.007091277425658 0.000554833974109 0.000500018986081 ;
+    0.194564361178625 0.137570544970062 0.079900546066037 0.025549350260855 0.005673491663727 0.000408759661293 0.000445141805523 ;
+    0.097751022951371 0.068620792366922 0.039365625553724 0.012397445282748 0.002818297719268 0.000149127033057 0.000311439525732 ;
+    0.030674916797523 0.021423584760802 0.012211671337997 0.003838328661082 0.000895825759878 0.000045643481270 0.000159762122747 ;
+    0.006652401806692 0.004725322877613 0.002773765135978 0.000908975538425 0.000210924069594 0.000037798409136 0.000057517483301 ;
+    0.000530676740860 0.000337644325418 0.000161654162384 0.000044642657051 0.000030848860156 0.000021605350900 0.000014496785601 ;
+    0.000447615535468 0.000364801746561 0.000270093752229 0.000138454732861 0.000052917434226 0.000014331973588 -0.000000022222325])
 ###########################################################################################
 # to modify the test, with a new expected f, print the new f using the following commands
 # in an interative Julia REPL. The path is the path to the .dfns file. 
 ########################################################################################## 
 """
+using moment_kinetics.load_data: open_readonly_output_file, load_pdf_data, load_ion_moments_data, load_fields_data
 fid = open_readonly_output_file(path, "dfns")
 f_ion_vpavperpzrst = load_pdf_data(fid)
 f_ion = f_ion_vpavperpzrst[:,:,1,1,1,:]
 ntind = 2
 nvpa = 13  #subject to grid choices
 nvperp = 7 #subject to grid choices
+# pdf
 for k in 1:ntind
-  for j in 1:nvperp-1
-      for i in 1:nvpa-1
-         @printf("%.15f ", f_ion[i,j,k])
-         print("; ")
-      end
-      @printf("%.15f ", f_ion[nvpa,j,k])
-      print(";;\n")
-  end
-  for i in 1:nvpa-1
-    @printf("%.15f ", f_ion[i,nvperp,k])
-    print("; ")
-  end
-  @printf("%.15f ", f_ion[nvpa,nvperp,k])
-  if k < ntind
-      print(";;;\n")
-  end  
+     for i in 1:nvpa-1
+         for j in 1:nvperp-1
+            @printf("%.15f ", f_ion[i,j,k])
+         end
+         @printf("%.15f ", f_ion[i,nvperp,k])
+         print(";\n")
+     end
+     for j in 1:nvperp-1
+       @printf("%.15f ", f_ion[nvpa,j,k])
+     end
+     @printf("%.15f ", f_ion[nvpa,nvperp,k])
+     if k < ntind
+         print(";;;\n")
+     end
+end
+# a moment
+n_ion_zrst, upar_ion_zrst, ppar_ion_zrst, pperp_ion_zrst, qpar_ion_zrst, v_t_ion_zrst, dSdt_zrst = load_ion_moments_data(fid,extended_moments=true)
+for k in 1:ntind
+   @printf("%.15f", n_ion_zrst[1,1,1,k])
+   if k < ntind
+       print(", ")
+   end
+end
+# a field
+phi_zrt, Er_zrt, Ez_zrt = load_fields_data(fid)
+for k in 1:ntind
+   @printf("%.15f", phi_zrt[1,1,k])
+   if k < ntind
+       print(", ")
+   end
 end
 """
 # default inputs for tests
@@ -128,7 +239,7 @@ test_input_gauss_legendre = Dict("run_name" => "gausslegendre_pseudospectral",
                               "vperp_nelement" => 3,
                               "vperp_L" => 3.0,
                               "vperp_discretization" => "gausslegendre_pseudospectral",
-                              "split_operators" => false,
+                              #"split_operators" => false,
                               "ionization_frequency" => 0.0,
                               "charge_exchange_frequency" => 0.0,
                               "constant_ionization_rate" => false,
@@ -175,14 +286,14 @@ test_input_gauss_legendre = Dict("run_name" => "gausslegendre_pseudospectral",
 Run a sound-wave test for a single set of parameters
 """
 # Note 'name' should not be shared by any two tests in this file
-function run_test(test_input, rtol, atol, upar_rtol=nothing; args...)
+function run_test(test_input, expected, rtol, atol, upar_rtol=nothing; args...)
     # by passing keyword arguments to run_test, args becomes a Dict which can be used to
     # update the default inputs
 
     # Make a copy to make sure nothing modifies the input Dicts defined in this test
     # script.
     test_input = deepcopy(test_input)
-
+    
     if upar_rtol === nothing
         upar_rtol = rtol
     end
@@ -204,10 +315,8 @@ function run_test(test_input, rtol, atol, upar_rtol=nothing; args...)
 
     # Update default inputs with values to be changed
     input = merge(test_input, modified_inputs)
-    #input = test_input
 
     input["run_name"] = name
-
     # Suppress console output while running
     quietoutput() do
         # run simulation
@@ -325,7 +434,26 @@ function runtests()
         # GaussLegendre pseudospectral
         # Benchmark data is taken from this run (GaussLegendre)
         @testset "Gauss Legendre base" begin
-            run_test(test_input_gauss_legendre, 1.e-14, 1.0e-14 )
+            run_name = "gausslegendre_pseudospectral"
+            vperp_bc = "zero-impose-regularity"
+            run_test(test_input_gauss_legendre,
+             expected_zero_impose_regularity, 1.0e-14, 1.0e-14;
+             vperp_bc=vperp_bc)
+        end
+        @testset "Gauss Legendre no enforced regularity condition at vperp = 0" begin
+            run_name = "gausslegendre_pseudospectral_no_regularity"
+            vperp_bc = "zero"
+            run_test(test_input_gauss_legendre,
+            expected_zero,
+             1.0e-14, 1.0e-14; vperp_bc=vperp_bc)
+        end
+        @testset "Gauss Legendre no (explicitly) enforced boundary conditions" begin
+            run_name = "gausslegendre_pseudospectral_none_bc"
+            vperp_bc = "none"
+            vpa_bc = "none"
+            run_test(test_input_gauss_legendre,
+            expected_none_bc,
+             1.0e-14, 1.0e-14; vperp_bc=vperp_bc, vpa_bc=vpa_bc)
         end
     end
 end
diff --git a/moment_kinetics/test/gyroaverage_tests.jl b/moment_kinetics/test/gyroaverage_tests.jl
index 38a6c7085..21c43f0b7 100644
--- a/moment_kinetics/test/gyroaverage_tests.jl
+++ b/moment_kinetics/test/gyroaverage_tests.jl
@@ -10,7 +10,7 @@ using SpecialFunctions: besselj0
 import moment_kinetics
 using moment_kinetics.input_structs
 using moment_kinetics.coordinates: define_coordinate
-using moment_kinetics.geo: init_magnetic_geometry
+using moment_kinetics.geo: init_magnetic_geometry, setup_geometry_input
 using moment_kinetics.communication
 using moment_kinetics.looping
 using moment_kinetics.array_allocation: allocate_float, allocate_shared_float
@@ -142,11 +142,9 @@ function gyroaverage_test(absolute_error; rhostar=0.1, pitch=0.5, ngrid=5, kr=2,
         gyrophase, gyrophase_spectral = define_coordinate(gyrophase_input; collision_operator_dim=false, run_directory=test_output_directory)
         
         # create test geometry
-        #rhostar = 0.1 #rhostar of ions for ExB drift
         option = "constant-helical"
-        #pitch = 1.0
-        DeltaB = 1.0
-        geometry_in = geometry_input(rhostar,option,pitch,DeltaB)
+        inputdict = Dict("geometry" => Dict("option" => option, "rhostar" => rhostar, "pitch" => pitch))
+        geometry_in = setup_geometry_input(inputdict, 0.1)
         geometry = init_magnetic_geometry(geometry_in,z,r)
         
         # create test composition
@@ -243,7 +241,7 @@ function create_test_composition()
     # wall potential at z = 0
     phi_wall = 0.0
     # constant to test nonzero Er
-    Er_constant = 0.0
+    #Er_constant = 0.0
     # ratio of the neutral particle mass to the ion particle mass
     mn_over_mi = 1.0
     # ratio of the electron particle mass to the ion particle mass
@@ -253,7 +251,7 @@ function create_test_composition()
     recycling_fraction = 1.0
     gyrokinetic_ions = true
     return composition = species_composition(n_species, n_ion_species, n_neutral_species,
-            electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, Er_constant,
+            electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, #Er_constant,
             mn_over_mi, me_over_mi, recycling_fraction, gyrokinetic_ions, allocate_float(n_species))
 end
 
diff --git a/plots_post_processing/plots_post_processing/src/plot_MMS_sequence.jl b/plots_post_processing/plots_post_processing/src/plot_MMS_sequence.jl
index 84bfcf21b..92d780b27 100644
--- a/plots_post_processing/plots_post_processing/src/plot_MMS_sequence.jl
+++ b/plots_post_processing/plots_post_processing/src/plot_MMS_sequence.jl
@@ -255,7 +255,7 @@ function get_MMS_error_data(path_list,scan_type,scan_name)
         dfnn_func = manufactured_solns_list.dfnn_func
         densn_func = manufactured_solns_list.densn_func
         
-        manufactured_E_fields = manufactured_electric_fields(Lr_in,z.L,r_bc,z_bc,composition,r.n,manufactured_solns_input,species)
+        manufactured_E_fields = manufactured_electric_fields(Lr_in,z.L,r_bc,z_bc,composition,geo_in,r.n,manufactured_solns_input,species)
         Er_func = manufactured_E_fields.Er_func
         Ez_func = manufactured_E_fields.Ez_func
         phi_func = manufactured_E_fields.phi_func
diff --git a/plots_post_processing/plots_post_processing/src/plots_post_processing.jl b/plots_post_processing/plots_post_processing/src/plots_post_processing.jl
index 972b92620..8fdb8fd81 100644
--- a/plots_post_processing/plots_post_processing/src/plots_post_processing.jl
+++ b/plots_post_processing/plots_post_processing/src/plots_post_processing.jl
@@ -1097,7 +1097,7 @@ function analyze_and_plot_data(prefix...; run_index=nothing)
         densn_func = manufactured_solns_list.densn_func
         manufactured_E_fields =
             manufactured_electric_fields(Lr_in, z_global.L, r_global.bc, z_global.bc,
-                                         composition, r_global.n,
+                                         composition, geometry.input, r_global.n,
                                          manufactured_solns_input, species)
         Er_func = manufactured_E_fields.Er_func
         Ez_func = manufactured_E_fields.Ez_func
@@ -2806,13 +2806,13 @@ function plot_fields_rt(phi, delta_phi, time, itime_min, itime_max, nwrite_movie
         outfile = string(run_name, "_delta_phi(r0,z0)_vs_t.pdf")
         trysavefig(outfile)
     end
-    if pp.plot_phi_vs_z_t
+    if pp.plot_phi_vs_z_t && r.n > 1
         # make a heatmap plot of ϕ(r,t)
         heatmap(time, r.grid, phi, xlabel="time", ylabel="r", title="ϕ", c = :deep)
         outfile = string(run_name, "_phi_vs_r_t.pdf")
         trysavefig(outfile)
     end
-    if pp.animate_phi_vs_z
+    if pp.animate_phi_vs_z && r.n > 1
         # make a gif animation of ϕ(r) at different times
         anim = @animate for i ∈ itime_min:nwrite_movie:itime_max
             @views plot(r.grid, phi[:,i], xlabel="r", ylabel="ϕ", ylims = (phimin,phimax))
@@ -2826,10 +2826,11 @@ function plot_fields_rt(phi, delta_phi, time, itime_min, itime_max, nwrite_movie
     # plot!(exp.(-(phi[cld(nz,2),end] .- phi[izmid:end,end])) .* erfi.(sqrt.(abs.(phi[cld(nz,2),end] .- phi[izmid:end,end])))/sqrt(pi)/0.688, phi[izmid:end,end] .- phi[izmid,end], label = "analytical", linewidth=2)
     # outfile = string(run_name, "_harrison_comparison.pdf")
     # trysavefig(outfile)
-    plot(r.grid, phi[:,end], xlabel="r/Lr", ylabel="eϕ/Te", label="", linewidth=2)
-    outfile = string(run_name, "_phi(r)_final.pdf")
-    trysavefig(outfile)
-
+    if r.n > 1
+        plot(r.grid, phi[:,end], xlabel="r/Lr", ylabel="eϕ/Te", label="", linewidth=2)
+        outfile = string(run_name, "_phi(r)_final.pdf")
+        trysavefig(outfile)
+    end
     println("done.")
 end
 
@@ -2858,7 +2859,7 @@ function plot_ion_pdf(run_name, run_name_label, vpa, vperp, z, r, z_local, r_loc
         spec_string = [string("_", spec_type)]
     end
     # make a gif animation of f(vpa,z,t) at a given (vperp,r) location
-    if pp.animate_f_vs_vpa_z
+    if pp.animate_f_vs_vpa_z && (vpa.n > 1 && z_local.n > 1)
         pdf = load_distributed_ion_pdf_slice(run_name, nblocks,
                                                  itime_min:iskip:itime_max, n_species,
                                                  r_local, z_local, vperp, vpa;
@@ -2876,7 +2877,7 @@ function plot_ion_pdf(run_name, run_name_label, vpa, vperp, z, r, z_local, r_loc
         end
     end
     # make a gif animation of f(vpa,r,t) at a given (vperp,z) location
-    if pp.animate_f_vs_vpa_r
+    if pp.animate_f_vs_vpa_r && (vpa.n > 1 && r_local.n > 1)
         pdf = load_distributed_ion_pdf_slice(run_name, nblocks,
                                              itime_min:iskip:itime_max, n_species,
                                              r_local, z_local, vperp, vpa;
@@ -2894,7 +2895,7 @@ function plot_ion_pdf(run_name, run_name_label, vpa, vperp, z, r, z_local, r_loc
         end
     end
     # make a gif animation of f(vperp,z,t) at a given (vpa,r) location
-    if pp.animate_f_vs_vperp_z
+    if pp.animate_f_vs_vperp_z && (vperp.n > 1 && z_local.n > 1)
         pdf = load_distributed_ion_pdf_slice(run_name, nblocks,
                                                  itime_min:iskip:itime_max, n_species,
                                                  r_local, z_local, vperp, vpa; ivpa=ivpa0,
@@ -2908,7 +2909,7 @@ function plot_ion_pdf(run_name, run_name_label, vpa, vperp, z, r, z_local, r_loc
         end
     end
     # make a gif animation of f(vperp,r,t) at a given (vpa,z) location
-    if pp.animate_f_vs_vperp_r
+    if pp.animate_f_vs_vperp_r && (vperp.n > 1 && r_local.n > 1)
         pdf = load_distributed_ion_pdf_slice(run_name, nblocks,
                                                  itime_min:iskip:itime_max, n_species,
                                                  r_local, z_local, vperp, vpa; ivpa=ivpa0,
@@ -2970,7 +2971,7 @@ function plot_ion_pdf(run_name, run_name_label, vpa, vperp, z, r, z_local, r_loc
         end
     end
     # make a gif animation of f(z,r,t) at a given (vpa,vperp) location
-    if pp.animate_f_vs_r_z
+    if pp.animate_f_vs_r_z && (z_local.n > 1 && r_local.n > 1)
         pdf = load_distributed_ion_pdf_slice(run_name, nblocks,
                                              itime_min:iskip:itime_max, n_species,
                                              r_local, z_local, vperp, vpa; ivpa=ivpa0,
@@ -3067,15 +3068,16 @@ end
 function plot_fields_2D(phi, Ez, Er, time, z, r, iz0, ir0,
     itime_min, itime_max, nwrite_movie, run_name, pp, description)
     nr = size(r,1)
+    nz = size(z,1)
     print("Plotting fields data...")
     phimin = minimum(phi)
     phimax = maximum(phi)
-    if pp.plot_phi_vs_r0_z # plot last timestep phi[z,ir0]
+    if pp.plot_phi_vs_r0_z && nz > 1 # plot last timestep phi[z,ir0]
         @views plot(z, phi[:,ir0,end], xlabel=L"z/L_z", ylabel=L"\phi")
         outfile = string(run_name, "_phi"*description*"(r0,z)_vs_z.pdf")
         trysavefig(outfile)
     end
-    if pp.animate_phi_vs_r_z && nr > 1
+    if pp.animate_phi_vs_r_z && nr > 1 && nz > 1
         # make a gif animation of ϕ(z) at different times
         anim = @animate for i ∈ itime_min:nwrite_movie:itime_max
             @views heatmap(r, z, phi[:,:,i], xlabel="r", ylabel="z", c = :deep, interpolation = :cubic)
@@ -3087,7 +3089,7 @@ function plot_fields_2D(phi, Ez, Er, time, z, r, iz0, ir0,
         end
         outfile = string(run_name, "_phi(zwall-)_vs_r.gif")
         trygif(anim, outfile, fps=5)
-    elseif pp.animate_phi_vs_r_z && nr == 1 # make a gif animation of ϕ(z) at different times
+    elseif pp.animate_phi_vs_r_z && nr == 1 && nz > 1 # make a gif animation of ϕ(z) at different times
         anim = @animate for i ∈ itime_min:nwrite_movie:itime_max
             @views plot(z, phi[:,1,i], xlabel="z", ylabel=L"\widetilde{\phi}", ylims = (phimin,phimax))
         end
@@ -3096,17 +3098,17 @@ function plot_fields_2D(phi, Ez, Er, time, z, r, iz0, ir0,
     end
     Ezmin = minimum(Ez)
     Ezmax = maximum(Ez)
-    if pp.plot_Ez_vs_r0_z # plot last timestep Ez[z,ir0]
+    if pp.plot_Ez_vs_r0_z && nz > 1 # plot last timestep Ez[z,ir0]
         @views plot(z, Ez[:,ir0,end], xlabel=L"z/L_z", ylabel=L"E_z")
         outfile = string(run_name, "_Ez"*description*"(r0,z)_vs_z.pdf")
         trysavefig(outfile)
     end
-    if pp.plot_wall_Ez_vs_r && nr > 1 # plot last timestep Ez[z_wall,r]
+    if pp.plot_wall_Ez_vs_r && nr > 1 && nz > 1 # plot last timestep Ez[z_wall,r]
         @views plot(r, Ez[1,:,end], xlabel=L"r/L_r", ylabel=L"E_z")
         outfile = string(run_name, "_Ez"*description*"(r,z_wall-)_vs_r.pdf")
         trysavefig(outfile)
     end
-    if pp.animate_Ez_vs_r_z && nr > 1
+    if pp.animate_Ez_vs_r_z && nr > 1 && nz > 1
         # make a gif animation of ϕ(z) at different times
         anim = @animate for i ∈ itime_min:nwrite_movie:itime_max
             @views heatmap(r, z, Ez[:,:,i], xlabel="r", ylabel="z", c = :deep, interpolation = :cubic)
@@ -3118,7 +3120,7 @@ function plot_fields_2D(phi, Ez, Er, time, z, r, iz0, ir0,
         end
         outfile = string(run_name, "_Ez(zwall-)_vs_r.gif")
         trygif(anim, outfile, fps=5)
-    elseif pp.animate_Ez_vs_r_z && nr == 1
+    elseif pp.animate_Ez_vs_r_z && nr == 1 && nz > 1
         anim = @animate for i ∈ itime_min:nwrite_movie:itime_max
             @views plot(z, Ez[:,1,i], xlabel="z", ylabel=L"\widetilde{E}_z", ylims = (Ezmin,Ezmax))
         end
@@ -3127,12 +3129,12 @@ function plot_fields_2D(phi, Ez, Er, time, z, r, iz0, ir0,
     end
     Ermin = minimum(Er)
     Ermax = maximum(Er)
-    if pp.plot_Er_vs_r0_z # plot last timestep Er[z,ir0]
+    if pp.plot_Er_vs_r0_z && nz > 1 # plot last timestep Er[z,ir0]
         @views plot(z, Er[:,ir0,end], xlabel=L"z/L_z", ylabel=L"E_r")
         outfile = string(run_name, "_Er"*description*"(r0,z)_vs_z.pdf")
         trysavefig(outfile)
     end
-    if pp.plot_wall_Er_vs_r && nr > 1 # plot last timestep Er[z_wall,r]
+    if pp.plot_wall_Er_vs_r && nr > 1 && nz > 1 # plot last timestep Er[z_wall,r]
         @views plot(r, Er[1,:,end], xlabel=L"r/L_r", ylabel=L"E_r")
         outfile = string(run_name, "_Er"*description*"(r,z_wall-)_vs_r.pdf")
         trysavefig(outfile)
@@ -3142,7 +3144,7 @@ function plot_fields_2D(phi, Ez, Er, time, z, r, iz0, ir0,
         outfile = string(run_name, "_Er(zwall-)_vs_r.gif")
         trygif(anim, outfile, fps=5)
     end
-    if pp.animate_Er_vs_r_z && nr > 1
+    if pp.animate_Er_vs_r_z && nr > 1 && nz > 1
         # make a gif animation of ϕ(z) at different times
         anim = @animate for i ∈ itime_min:nwrite_movie:itime_max
             @views heatmap(r, z, Er[:,:,i], xlabel="r", ylabel="z", c = :deep, interpolation = :cubic)
@@ -3159,6 +3161,7 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure,
     time, z, r, iz0, ir0, n_ion_species,
     itime_min, itime_max, nwrite_movie, run_name, pp)
     nr = size(r,1)
+    nz = size(z,1)
     ntime = size(time,1)
     print("Plotting ion moments data...")
     for is in 1:n_ion_species
@@ -3166,7 +3169,7 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure,
 		# the density
 		densitymin = minimum(density[:,:,is,:])
 		densitymax = maximum(density)
-		if pp.plot_density_vs_r0_z # plot last timestep density[z,ir0]
+		if pp.plot_density_vs_r0_z && nz > 1 # plot last timestep density[z,ir0]
 			@views plot(z, density[:,ir0,is,end], xlabel=L"z/L_z", ylabel=L"n_i")
 			outfile = string(run_name, "_density"*description*"(r0,z)_vs_z.pdf")
 			trysavefig(outfile)
@@ -3176,7 +3179,7 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure,
 			outfile = string(run_name, "_density"*description*"(r,z_wall)_vs_r.pdf")
 			trysavefig(outfile)
 		end
-		if pp.animate_density_vs_r_z && nr > 1
+		if pp.animate_density_vs_r_z && nr > 1 && nz > 1
 			# make a gif animation of ϕ(z) at different times
 			anim = @animate for i ∈ itime_min:nwrite_movie:itime_max
 				@views heatmap(r, z, density[:,:,is,i], xlabel="r", ylabel="z", c = :deep, interpolation = :cubic)
@@ -3184,7 +3187,7 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure,
 			outfile = string(run_name, "_density"*description*"_vs_r_z.gif")
 			trygif(anim, outfile, fps=5)
 		end
-		if pp.plot_density_vs_r_z && nr > 1
+		if pp.plot_density_vs_r_z && nr > 1 && nz > 1
 			@views heatmap(r, z, density[:,:,is,end], xlabel=L"r", ylabel=L"z", c = :deep, interpolation = :cubic,
 			windowsize = (360,240), margin = 15pt)
 			outfile = string(run_name, "_density"*description*"_vs_r_z.pdf")
@@ -3201,17 +3204,17 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure,
 		# the parallel flow
 		parallel_flowmin = minimum(parallel_flow[:,:,is,:])
 		parallel_flowmax = maximum(parallel_flow)
-		if pp.plot_parallel_flow_vs_r0_z # plot last timestep parallel_flow[z,ir0]
+		if pp.plot_parallel_flow_vs_r0_z && nz > 1 # plot last timestep parallel_flow[z,ir0]
 			@views plot(z, parallel_flow[:,ir0,is,end], xlabel=L"z/L_z", ylabel=L"u_{i\|\|}")
 			outfile = string(run_name, "_parallel_flow"*description*"(r0,z)_vs_z.pdf")
 			trysavefig(outfile)
 		end
-		if pp.plot_wall_parallel_flow_vs_r && nr > 1 # plot last timestep parallel_flow[z_wall,r]
+		if pp.plot_wall_parallel_flow_vs_r && nr > 1 && nz > 1# plot last timestep parallel_flow[z_wall,r]
 			@views plot(r, parallel_flow[end,:,is,end], xlabel=L"r/L_r", ylabel=L"u_{i\|\|}")
 			outfile = string(run_name, "_parallel_flow"*description*"(r,z_wall)_vs_r.pdf")
 			trysavefig(outfile)
 		end
-		if pp.animate_parallel_flow_vs_r_z && nr > 1
+		if pp.animate_parallel_flow_vs_r_z && nr > 1 && nz > 1
 			# make a gif animation of ϕ(z) at different times
 			anim = @animate for i ∈ itime_min:nwrite_movie:itime_max
 				@views heatmap(r, z, parallel_flow[:,:,is,i], xlabel="r", ylabel="z", c = :deep, interpolation = :cubic)
@@ -3219,7 +3222,7 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure,
 			outfile = string(run_name, "_parallel_flow"*description*"_vs_r_z.gif")
 			trygif(anim, outfile, fps=5)
 		end
-		if pp.plot_parallel_flow_vs_r_z && nr > 1
+		if pp.plot_parallel_flow_vs_r_z && nr > 1 && nz > 1
 			@views heatmap(r, z, parallel_flow[:,:,is,end], xlabel=L"r", ylabel=L"z", c = :deep, interpolation = :cubic,
 			windowsize = (360,240), margin = 15pt)
 			outfile = string(run_name, "_parallel_flow"*description*"_vs_r_z.pdf")
@@ -3236,7 +3239,7 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure,
 		# the parallel pressure
 		parallel_pressuremin = minimum(parallel_pressure[:,:,is,:])
 		parallel_pressuremax = maximum(parallel_pressure)
-		if pp.plot_parallel_pressure_vs_r0_z # plot last timestep parallel_pressure[z,ir0]
+		if pp.plot_parallel_pressure_vs_r0_z && nz > 1 # plot last timestep parallel_pressure[z,ir0]
 			@views plot(z, parallel_pressure[:,ir0,is,end], xlabel=L"z/L_z", ylabel=L"p_{i\|\|}")
 			outfile = string(run_name, "_parallel_pressure"*description*"(r0,z)_vs_z.pdf")
 			trysavefig(outfile)
@@ -3246,7 +3249,7 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure,
 			outfile = string(run_name, "_parallel_pressure"*description*"(r,z_wall)_vs_r.pdf")
 			trysavefig(outfile)
 		end
-		if pp.animate_parallel_pressure_vs_r_z && nr > 1
+		if pp.animate_parallel_pressure_vs_r_z && nr > 1 && nz > 1
 			# make a gif animation of ϕ(z) at different times
 			anim = @animate for i ∈ itime_min:nwrite_movie:itime_max
 				@views heatmap(r, z, parallel_pressure[:,:,is,i], xlabel="r", ylabel="z", c = :deep, interpolation = :cubic)
@@ -3254,7 +3257,7 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure,
 			outfile = string(run_name, "_parallel_pressure"*description*"_vs_r_z.gif")
 			trygif(anim, outfile, fps=5)
 		end
-		if pp.plot_parallel_pressure_vs_r_z && nr > 1
+		if pp.plot_parallel_pressure_vs_r_z && nr > 1 && nz > 1
 			@views heatmap(r, z, parallel_pressure[:,:,is,end], xlabel=L"r", ylabel=L"z", c = :deep, interpolation = :cubic,
 			windowsize = (360,240), margin = 15pt)
 			outfile = string(run_name, "_parallel_pressure"*description*"_vs_r_z.pdf")
@@ -3264,17 +3267,17 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure,
                 # Note factor of 2 here because currently temperatures are normalised by
                 # Tref, while pressures are normalised by m*nref*c_ref^2=2*nref*Tref
                 temperature = 2.0 * parallel_pressure ./ density
-                if pp.plot_parallel_temperature_vs_r0_z # plot last timestep parallel_temperature[z,ir0]
+                if pp.plot_parallel_temperature_vs_r0_z  && nz > 1# plot last timestep parallel_temperature[z,ir0]
                     @views plot(z, temperature[:,ir0,is,end], xlabel=L"z/L_z", ylabel=L"T_i")
                     outfile = string(run_name, "_temperature"*description*"(r0,z)_vs_z.pdf")
                     trysavefig(outfile)
                 end
-                if pp.plot_wall_parallel_temperature_vs_r && nr > 1 # plot last timestep parallel_temperature[z_wall,r]
+                if pp.plot_wall_parallel_temperature_vs_r && nr > 1 && nz > 1 # plot last timestep parallel_temperature[z_wall,r]
                     @views plot(r, temperature[end,:,is,end], xlabel=L"r/L_r", ylabel=L"T_i")
                     outfile = string(run_name, "_temperature"*description*"(r,z_wall)_vs_r.pdf")
                     trysavefig(outfile)
                 end
-                if pp.animate_parallel_temperature_vs_r_z && nr > 1
+                if pp.animate_parallel_temperature_vs_r_z && nr > 1 && nz > 1
                     # make a gif animation of T_||(z) at different times
                     anim = @animate for i ∈ itime_min:nwrite_movie:itime_max
                         @views heatmap(r, z, temperature[:,:,is,i], xlabel="r", ylabel="z", c = :deep, interpolation = :cubic)
@@ -3282,7 +3285,7 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure,
                     outfile = string(run_name, "_temperature"*description*"_vs_r_z.gif")
                     trygif(anim, outfile, fps=5)
                 end
-                if pp.plot_parallel_temperature_vs_r_z && nr > 1
+                if pp.plot_parallel_temperature_vs_r_z && nr > 1 && nz > 1
                     @views heatmap(r, z, temperature[:,:,is,end], xlabel=L"r", ylabel=L"z", c = :deep, interpolation = :cubic,
                                   windowsize = (360,240), margin = 15pt)
                     outfile = string(run_name, "_temperature"*description*"_vs_r_z.pdf")
@@ -3305,33 +3308,47 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure,
 			outfile = string(run_name, "_delta_perpendicular_pressure"*description*"(iz0,ir0)_vs_t.pdf")
 			trysavefig(outfile)
         end
+        if pp.plot_perpendicular_pressure_vs_r0_z && nz > 1 # plot last timestep perpendicular_pressure[z,ir0]
+           @views plot(z, perpendicular_pressure[:,ir0,is,end], xlabel=L"z/L_z", ylabel=L"p_{i\perp}",label="")
+           outfile = string(run_name, "_perpendicular_pressure"*description*"(r0,z)_vs_z.pdf")
+           trysavefig(outfile)
+        end
         # the total pressure
+        total_pressure = copy(parallel_pressure)
+        @. total_pressure = (2.0/3.0)*perpendicular_pressure + (1.0/3.0)*parallel_pressure
         if pp.plot_ppar0_vs_t && pp.plot_pperp0_vs_t
+            
             @views plot([time, time, time] , 
             [parallel_pressure[iz0,ir0,is,:], perpendicular_pressure[iz0,ir0,is,:], 
-            (2.0/3.0).*perpendicular_pressure[iz0,ir0,is,:] .+ (1.0/3.0).*parallel_pressure[iz0,ir0,is,:]],
-            xlabel=L"t/ (L_{ref}/c_{ref})", ylabel="", label = [L"p_{i\|\|}(t)" L"p_{i\perp}(t)" L"p_{i}(t)"])
+            total_pressure[iz0,ir0,is,:]],
+            xlabel=L"t/ (L_{ref}/c_{ref})", ylabel="", label = [L"p_{i\|\|}(t)" L"p_{i\perp}(t)" L"p_{i}(t)"],
+            foreground_color_legend = nothing, background_color_legend = nothing)
 			outfile = string(run_name, "_pressures"*description*"(iz0,ir0)_vs_t.pdf")
 			trysavefig(outfile)
             @views plot([time, time, time] , 
             [parallel_pressure[iz0,ir0,is,:] .- parallel_pressure[iz0,ir0,is,1], perpendicular_pressure[iz0,ir0,is,:] .- perpendicular_pressure[iz0,ir0,is,1], 
-            (2.0/3.0).*(perpendicular_pressure[iz0,ir0,is,:] .- perpendicular_pressure[iz0,ir0,is,1]).+
-            (1.0/3.0).*(parallel_pressure[iz0,ir0,is,:] .- parallel_pressure[iz0,ir0,is,1])],
-            xlabel=L"t/ (L_{ref}/c_{ref})", ylabel="", label = [L"p_{i\|\|}(t) - p_{i\|\|}(0)" L"p_{i\perp}(t) - p_{i\perp}(0)" L"p_{i}(t) - p_{i}(0)"])
+             total_pressure[iz0,ir0,is,:] .- total_pressure[iz0,ir0,is,1]],
+             xlabel=L"t/ (L_{ref}/c_{ref})", ylabel="", label = [L"p_{i\|\|}(t) - p_{i\|\|}(0)" L"p_{i\perp}(t) - p_{i\perp}(0)" L"p_{i}(t) - p_{i}(0)"],
+             foreground_color_legend = nothing, background_color_legend = nothing)
 			outfile = string(run_name, "_delta_pressures"*description*"(iz0,ir0)_vs_t.pdf")
 			trysavefig(outfile)
             @views plot([time] , 
-            [(2.0/3.0).*perpendicular_pressure[iz0,ir0,is,:] .+ (1.0/3.0).*parallel_pressure[iz0,ir0,is,:]],
+            [total_pressure[iz0,ir0,is,:]],
             xlabel=L"t/ (L_{ref}/c_{ref})", ylabel=L"p_{i}(t)", label = "")
 			outfile = string(run_name, "_pressure"*description*"(iz0,ir0)_vs_t.pdf")
 			trysavefig(outfile)
             @views plot([time] , 
-            [(2.0/3.0).*(perpendicular_pressure[iz0,ir0,is,:] .- perpendicular_pressure[iz0,ir0,is,1]).+
-            (1.0/3.0).*(parallel_pressure[iz0,ir0,is,:] .- parallel_pressure[iz0,ir0,is,1])],
+            [total_pressure[iz0,ir0,is,:] .- total_pressure[iz0,ir0,is,1]],
             xlabel=L"t/ (L_{ref}/c_{ref})", ylabel=L"p_{i}(t) - p_{i}(0)", label = "")
 			outfile = string(run_name, "_delta_pressure"*description*"(iz0,ir0)_vs_t.pdf")
 			trysavefig(outfile)
         end
+        if pp.plot_perpendicular_pressure_vs_r0_z && pp.plot_parallel_pressure_vs_r0_z && nz > 1 # plot last timestep perpendicular_pressure[z,ir0]
+           @views plot([z,z,z], [parallel_pressure[:,ir0,is,end],perpendicular_pressure[:,ir0,is,end],total_pressure[:,ir0,is,end]], xlabel=L"z/L_z", ylabel="",
+                       label = [L"p_{i\|\|}" L"p_{i\perp}" L"p_{i}"], foreground_color_legend = nothing, background_color_legend = nothing)
+           outfile = string(run_name, "_all_pressures"*description*"(r0,z)_vs_z.pdf")
+           trysavefig(outfile)
+        end
         # the thermal speed
         if pp.plot_vth0_vs_t
             @views plot(time, thermal_speed[iz0,ir0,is,:], xlabel=L"t / (L_{ref}/c_{ref})", ylabel=L"v_{i,th}(t)", label = "")
diff --git a/plots_post_processing/plots_post_processing/src/post_processing_input.jl b/plots_post_processing/plots_post_processing/src/post_processing_input.jl
index 5c3b6f58d..08cad55f8 100644
--- a/plots_post_processing/plots_post_processing/src/post_processing_input.jl
+++ b/plots_post_processing/plots_post_processing/src/post_processing_input.jl
@@ -102,6 +102,7 @@ const plot_wall_density_vs_r = true # plot last timestep density[z_wall,r]
 const plot_density_vs_r_z = true
 const animate_density_vs_r_z = true
 const plot_parallel_flow_vs_r0_z = true # plot last timestep parallel_flow[z,ir0]
+const plot_perpendicular_flow_vs_r0_z = true # plot last timestep perpendicular_flow[z,ir0]
 const plot_wall_parallel_flow_vs_r = true # plot last timestep parallel_flow[z_wall,r]
 const plot_parallel_flow_vs_r_z = true
 const animate_parallel_flow_vs_r_z = true
@@ -171,7 +172,8 @@ pp = pp_input(calculate_frequencies, plot_phi0_vs_t, plot_phi_vs_z_t, animate_ph
     animate_Er_vs_r_z, animate_Ez_vs_r_z, animate_phi_vs_r_z, plot_phi_vs_r0_z,
     plot_Ez_vs_r0_z, plot_wall_Ez_vs_r, plot_Er_vs_r0_z, plot_wall_Er_vs_r,
     plot_density_vs_r0_z, plot_wall_density_vs_r, plot_density_vs_r_z,
-    animate_density_vs_r_z, plot_parallel_flow_vs_r0_z, plot_wall_parallel_flow_vs_r,
+    animate_density_vs_r_z, plot_parallel_flow_vs_r0_z, plot_perpendicular_flow_vs_r0_z, 
+    plot_wall_parallel_flow_vs_r,
     plot_parallel_flow_vs_r_z, animate_parallel_flow_vs_r_z,
     plot_parallel_pressure_vs_r0_z, plot_wall_parallel_pressure_vs_r,
     plot_parallel_pressure_vs_r_z, animate_parallel_pressure_vs_r_z,
diff --git a/runs/1D-mirror_MMS_new_nel_r_1_z_16_vpa_16_vperp_8_diss.toml b/runs/1D-mirror_MMS_new_nel_r_1_z_16_vpa_16_vperp_8_diss.toml
deleted file mode 100644
index 9a475346b..000000000
--- a/runs/1D-mirror_MMS_new_nel_r_1_z_16_vpa_16_vperp_8_diss.toml
+++ /dev/null
@@ -1,98 +0,0 @@
-n_ion_species = 1
-n_neutral_species = 0
-electron_physics = "boltzmann_electron_response"
-#electron_physics = "boltzmann_electron_response_with_simple_sheath"
-evolve_moments_density = false
-evolve_moments_parallel_flow = false
-evolve_moments_parallel_pressure = false
-evolve_moments_conservation = false
-force_Er_zero_at_wall = false #true
-geometry_option="1D-mirror"
-DeltaB=10.0
-#geometry_option="constant-helical"
-#pitch=1.0
-Er_constant = 0.0
-T_e = 1.0
-T_wall = 1.0
-rhostar = 1.0
-initial_density1 = 0.5
-initial_temperature1 = 1.0
-initial_density2 = 0.5
-initial_temperature2 = 1.0
-z_IC_option1 = "sinusoid"
-z_IC_density_amplitude1 = 0.001
-z_IC_density_phase1 = 0.0
-z_IC_upar_amplitude1 = 0.0
-z_IC_upar_phase1 = 0.0
-z_IC_temperature_amplitude1 = 0.0
-z_IC_temperature_phase1 = 0.0
-z_IC_option2 = "sinusoid"
-z_IC_density_amplitude2 = 0.001
-z_IC_density_phase2 = 0.0
-z_IC_upar_amplitude2 = 0.0
-z_IC_upar_phase2 = 0.0
-z_IC_temperature_amplitude2 = 0.0
-z_IC_temperature_phase2 = 0.0
-charge_exchange_frequency = 0.0
-ionization_frequency = 0.0
-nstep = 2000
-dt = 0.0005
-nwrite = 200
-nwrite_dfns = 200
-use_semi_lagrange = false
-n_rk_stages = 4
-split_operators = false
-z_ngrid = 17
-z_nelement = 16
-z_nelement_local = 16
-z_bc = "wall"
-z_discretization = "chebyshev_pseudospectral"
-r_ngrid = 1
-r_nelement = 1
-r_nelement_local = 1
-r_bc = "periodic"
-r_discretization = "chebyshev_pseudospectral"
-vpa_ngrid = 17
-vpa_nelement = 16
-vpa_L = 12.0
-vpa_bc = "zero"
-vpa_discretization = "chebyshev_pseudospectral"
-vperp_ngrid = 17
-vperp_nelement = 8
-vperp_L = 6.0
-#vperp_discretization = "finite_difference"
-vperp_discretization = "chebyshev_pseudospectral"
-vperp_bc = "zero"
-
-vz_ngrid = 17
-vz_nelement = 4
-vz_L = 12.0
-vz_bc = "periodic"
-vz_discretization = "chebyshev_pseudospectral"
-
-vr_ngrid = 17
-vr_nelement = 4
-vr_L = 12.0
-vr_bc = "periodic"
-vr_discretization = "chebyshev_pseudospectral"
-
-vzeta_ngrid = 17
-vzeta_nelement = 4
-vzeta_L = 12.0
-vzeta_bc = "periodic"
-vzeta_discretization = "chebyshev_pseudospectral"
-
-[manufactured_solns]
- use_for_advance=true
- use_for_init=true
- # constant to be used to control Ez divergence in MMS tests
- epsilon_offset=0.1
- # bool to control if dfni is a function of vpa or vpabar in MMS test
- use_vpabar_in_mms_dfni=true
- alpha_switch=1.0
- type="default"
-[ion_numerical_dissipation]
-vpa_dissipation_coefficient = 0.001
-vperp_dissipation_coefficient = 0.001
-#z_dissipation_coefficient = 0.1
-r_dissipation_coefficient = 0.0
diff --git a/runs/1D-mirror_MMS_new_nel_r_1_z_6_vpa_6_vperp_3_diss.toml b/runs/1D-mirror_MMS_new_nel_r_1_z_6_vpa_6_vperp_3_diss.toml
deleted file mode 100644
index 347dd2b04..000000000
--- a/runs/1D-mirror_MMS_new_nel_r_1_z_6_vpa_6_vperp_3_diss.toml
+++ /dev/null
@@ -1,99 +0,0 @@
-n_ion_species = 1
-n_neutral_species = 0
-electron_physics = "boltzmann_electron_response"
-#electron_physics = "boltzmann_electron_response_with_simple_sheath"
-evolve_moments_density = false
-evolve_moments_parallel_flow = false
-evolve_moments_parallel_pressure = false
-evolve_moments_conservation = false
-force_Er_zero_at_wall = false #true
-geometry_option="1D-mirror"
-DeltaB=0.5
-#geometry_option="constant-helical"
-#pitch=1.0
-Er_constant = 0.0
-T_e = 1.0
-T_wall = 1.0
-rhostar = 1.0
-initial_density1 = 0.5
-initial_temperature1 = 1.0
-initial_density2 = 0.5
-initial_temperature2 = 1.0
-z_IC_option1 = "sinusoid"
-z_IC_density_amplitude1 = 0.001
-z_IC_density_phase1 = 0.0
-z_IC_upar_amplitude1 = 0.0
-z_IC_upar_phase1 = 0.0
-z_IC_temperature_amplitude1 = 0.0
-z_IC_temperature_phase1 = 0.0
-z_IC_option2 = "sinusoid"
-z_IC_density_amplitude2 = 0.001
-z_IC_density_phase2 = 0.0
-z_IC_upar_amplitude2 = 0.0
-z_IC_upar_phase2 = 0.0
-z_IC_temperature_amplitude2 = 0.0
-z_IC_temperature_phase2 = 0.0
-charge_exchange_frequency = 0.0
-ionization_frequency = 0.0
-nstep = 2000
-dt = 0.0005
-nwrite = 200
-nwrite_dfns = 200
-use_semi_lagrange = false
-n_rk_stages = 4
-split_operators = false
-z_ngrid = 17
-z_nelement = 6
-z_nelement_local = 6
-z_bc = "wall"
-z_discretization = "chebyshev_pseudospectral"
-r_ngrid = 1
-r_nelement = 1
-r_nelement_local = 1
-r_bc = "periodic"
-r_discretization = "chebyshev_pseudospectral"
-vpa_ngrid = 17
-vpa_nelement = 6
-vpa_L = 12.0
-vpa_bc = "zero"
-vpa_discretization = "chebyshev_pseudospectral"
-vperp_ngrid = 17
-vperp_nelement = 3
-vperp_L = 6.0
-#vperp_discretization = "finite_difference"
-#vperp_discretization = "chebyshev_pseudospectral"
-vperp_discretization = "gausslegendre_pseudospectral"
-vperp_bc = "zero"
-
-vz_ngrid = 17
-vz_nelement = 4
-vz_L = 12.0
-vz_bc = "periodic"
-vz_discretization = "chebyshev_pseudospectral"
-
-vr_ngrid = 17
-vr_nelement = 4
-vr_L = 12.0
-vr_bc = "periodic"
-vr_discretization = "chebyshev_pseudospectral"
-
-vzeta_ngrid = 17
-vzeta_nelement = 4
-vzeta_L = 12.0
-vzeta_bc = "periodic"
-vzeta_discretization = "chebyshev_pseudospectral"
-
-[manufactured_solns]
- use_for_advance=true
- use_for_init=true
- # constant to be used to control Ez divergence in MMS tests
- epsilon_offset=0.1
- # bool to control if dfni is a function of vpa or vpabar in MMS test
- use_vpabar_in_mms_dfni=true
- alpha_switch=1.0
- type="default"
-[ion_numerical_dissipation]
-vpa_dissipation_coefficient = 0.001
-vperp_dissipation_coefficient = 0.001
-#z_dissipation_coefficient = 0.1
-r_dissipation_coefficient = 0.0
diff --git a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_16_vpa_16_vperp_8_diss.toml b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_16_vpa_16_vperp_8_diss.toml
index 883a73275..f22e52894 100644
--- a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_16_vpa_16_vperp_8_diss.toml
+++ b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_16_vpa_16_vperp_8_diss.toml
@@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false
 evolve_moments_parallel_pressure = false
 evolve_moments_conservation = false
 force_Er_zero_at_wall = false #true
-Er_constant = 0.0
 T_e = 1.0
 T_wall = 1.0
 initial_density1 = 0.5
@@ -99,3 +98,4 @@ DeltaB=0.5
 #option="constant-helical"
 pitch=1.0
 rhostar = 1.0
+Er_constant = 0.0
diff --git a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_32_vpa_32_vperp_16_diss.toml b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_32_vpa_32_vperp_16_diss.toml
index 50e878cb3..fefea7811 100644
--- a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_32_vpa_32_vperp_16_diss.toml
+++ b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_32_vpa_32_vperp_16_diss.toml
@@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false
 evolve_moments_parallel_pressure = false
 evolve_moments_conservation = false
 force_Er_zero_at_wall = false #true
-Er_constant = 0.0
 T_e = 1.0
 T_wall = 1.0
 initial_density1 = 0.5
@@ -99,3 +98,4 @@ DeltaB=0.5
 #option="constant-helical"
 pitch=1.0
 rhostar = 1.0
+Er_constant = 0.0
diff --git a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_4_vpa_4_vperp_2_diss.toml b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_4_vpa_4_vperp_2_diss.toml
index 1624ca129..58051541c 100644
--- a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_4_vpa_4_vperp_2_diss.toml
+++ b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_4_vpa_4_vperp_2_diss.toml
@@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false
 evolve_moments_parallel_pressure = false
 evolve_moments_conservation = false
 force_Er_zero_at_wall = false #true
-Er_constant = 0.0
 T_e = 1.0
 T_wall = 1.0
 initial_density1 = 0.5
@@ -99,3 +98,4 @@ DeltaB=0.5
 #option="constant-helical"
 pitch=1.0
 rhostar = 1.0
+Er_constant = 0.0
diff --git a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_64_vpa_64_vperp_32_diss.toml b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_64_vpa_64_vperp_32_diss.toml
index 74e3af388..0cc6b6ef1 100644
--- a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_64_vpa_64_vperp_32_diss.toml
+++ b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_64_vpa_64_vperp_32_diss.toml
@@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false
 evolve_moments_parallel_pressure = false
 evolve_moments_conservation = false
 force_Er_zero_at_wall = false #true
-Er_constant = 0.0
 T_e = 1.0
 T_wall = 1.0
 initial_density1 = 0.5
@@ -99,3 +98,4 @@ DeltaB=0.5
 #option="constant-helical"
 pitch=1.0
 rhostar = 1.0
+Er_constant = 0.0
diff --git a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_8_vpa_8_vperp_4_diss.toml b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_8_vpa_8_vperp_4_diss.toml
index f4d4956ab..0f3a24c74 100644
--- a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_8_vpa_8_vperp_4_diss.toml
+++ b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_8_vpa_8_vperp_4_diss.toml
@@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false
 evolve_moments_parallel_pressure = false
 evolve_moments_conservation = false
 force_Er_zero_at_wall = false #true
-Er_constant = 0.0
 T_e = 1.0
 T_wall = 1.0
 initial_density1 = 0.5
@@ -99,3 +98,4 @@ DeltaB=0.5
 #option="constant-helical"
 pitch=1.0
 rhostar = 1.0
+Er_constant = 0.0
diff --git a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_diss.toml b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_diss.toml
index 8bf508ee1..96611183c 100644
--- a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_diss.toml
+++ b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_diss.toml
@@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false
 evolve_moments_parallel_pressure = false
 evolve_moments_conservation = false
 force_Er_zero_at_wall = false #true
-Er_constant = 0.0
 T_e = 1.0
 T_wall = 1.0
 initial_density1 = 0.5
@@ -89,3 +88,5 @@ split_operators = false
 vpa_dissipation_coefficient = 0.01
 #z_dissipation_coefficient = 0.1
 #r_dissipation_coefficient = 0.0
+[geometry]
+Er_constant = 0.0
diff --git a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml
index 1abe96001..346996f69 100644
--- a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml
+++ b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml
@@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false
 evolve_moments_parallel_pressure = false
 evolve_moments_conservation = false
 force_Er_zero_at_wall = false #true
-Er_constant = 0.0
 T_e = 1.0
 T_wall = 1.0
 initial_density1 = 0.5
@@ -93,3 +92,5 @@ split_operators = false
 vpa_dissipation_coefficient = -1.0
 z_dissipation_coefficient = -1.0
 r_dissipation_coefficient = -1.0
+[geometry]
+Er_constant = 0.0
diff --git a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml
index d9ddc1aeb..3378a93b1 100644
--- a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml
+++ b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml
@@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false
 evolve_moments_parallel_pressure = false
 evolve_moments_conservation = false
 force_Er_zero_at_wall = false #true
-Er_constant = 0.0
 T_e = 1.0
 T_wall = 1.0
 initial_density1 = 0.5
@@ -96,3 +95,5 @@ split_operators = false
 vpa_dissipation_coefficient = -1.0
 z_dissipation_coefficient = -1.0
 r_dissipation_coefficient = -1.0
+[geometry]
+Er_constant = 0.0
diff --git a/runs/2D-mirror_MMS_ngrid_5_nel_r_16_z_16_vpa_16_vperp_8_diss.toml b/runs/2D-mirror_MMS_ngrid_5_nel_r_16_z_16_vpa_16_vperp_8_diss.toml
index a4a31ccbb..29542c43f 100644
--- a/runs/2D-mirror_MMS_ngrid_5_nel_r_16_z_16_vpa_16_vperp_8_diss.toml
+++ b/runs/2D-mirror_MMS_ngrid_5_nel_r_16_z_16_vpa_16_vperp_8_diss.toml
@@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false
 evolve_moments_parallel_pressure = false
 evolve_moments_conservation = false
 force_Er_zero_at_wall = false #true
-Er_constant = 0.0
 T_e = 1.0
 T_wall = 1.0
 initial_density1 = 0.5
@@ -99,3 +98,4 @@ DeltaB=0.5
 #option="constant-helical"
 pitch=0.5
 rhostar = 1.0
+Er_constant = 0.0
diff --git a/runs/2D-mirror_MMS_ngrid_5_nel_r_32_z_32_vpa_32_vperp_16_diss.toml b/runs/2D-mirror_MMS_ngrid_5_nel_r_32_z_32_vpa_32_vperp_16_diss.toml
index 01dd1b487..3561f1d11 100644
--- a/runs/2D-mirror_MMS_ngrid_5_nel_r_32_z_32_vpa_32_vperp_16_diss.toml
+++ b/runs/2D-mirror_MMS_ngrid_5_nel_r_32_z_32_vpa_32_vperp_16_diss.toml
@@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false
 evolve_moments_parallel_pressure = false
 evolve_moments_conservation = false
 force_Er_zero_at_wall = false #true
-Er_constant = 0.0
 T_e = 1.0
 T_wall = 1.0
 initial_density1 = 0.5
@@ -99,3 +98,4 @@ DeltaB=0.5
 #option="constant-helical"
 pitch=0.5
 rhostar = 1.0
+Er_constant = 0.0
diff --git a/runs/2D-mirror_MMS_ngrid_5_nel_r_4_z_4_vpa_4_vperp_2_diss.toml b/runs/2D-mirror_MMS_ngrid_5_nel_r_4_z_4_vpa_4_vperp_2_diss.toml
index b13205d70..a8add1f3e 100644
--- a/runs/2D-mirror_MMS_ngrid_5_nel_r_4_z_4_vpa_4_vperp_2_diss.toml
+++ b/runs/2D-mirror_MMS_ngrid_5_nel_r_4_z_4_vpa_4_vperp_2_diss.toml
@@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false
 evolve_moments_parallel_pressure = false
 evolve_moments_conservation = false
 force_Er_zero_at_wall = false #true
-Er_constant = 0.0
 T_e = 1.0
 T_wall = 1.0
 initial_density1 = 0.5
@@ -99,3 +98,4 @@ DeltaB=0.5
 #option="constant-helical"
 pitch=0.5
 rhostar = 1.0
+Er_constant = 0.0
diff --git a/runs/2D-mirror_MMS_ngrid_5_nel_r_8_z_8_vpa_8_vperp_4_diss.toml b/runs/2D-mirror_MMS_ngrid_5_nel_r_8_z_8_vpa_8_vperp_4_diss.toml
index f19b91e39..2eb57f6ee 100644
--- a/runs/2D-mirror_MMS_ngrid_5_nel_r_8_z_8_vpa_8_vperp_4_diss.toml
+++ b/runs/2D-mirror_MMS_ngrid_5_nel_r_8_z_8_vpa_8_vperp_4_diss.toml
@@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false
 evolve_moments_parallel_pressure = false
 evolve_moments_conservation = false
 force_Er_zero_at_wall = false #true
-Er_constant = 0.0
 T_e = 1.0
 T_wall = 1.0
 initial_density1 = 0.5
@@ -99,3 +98,4 @@ DeltaB=0.5
 #option="constant-helical"
 pitch=0.5
 rhostar = 1.0
+Er_constant = 0.0
diff --git a/runs/2D-wall_MMS_nel_r_32_z_32_vpa_16_vperp_1_diss.toml b/runs/2D-wall_MMS_nel_r_32_z_32_vpa_16_vperp_1_diss.toml
index 274fc59e6..fa1f6c220 100644
--- a/runs/2D-wall_MMS_nel_r_32_z_32_vpa_16_vperp_1_diss.toml
+++ b/runs/2D-wall_MMS_nel_r_32_z_32_vpa_16_vperp_1_diss.toml
@@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false
 evolve_moments_parallel_pressure = false
 evolve_moments_conservation = false
 force_Er_zero_at_wall = false #true
-Er_constant = 0.0
 T_e = 1.0
 T_wall = 1.0
 initial_density1 = 0.5
@@ -95,3 +94,4 @@ r_dissipation_coefficient = 0.1
 option="constant-helical"
 pitch=0.5
 rhostar = 1.0
+Er_constant = 0.0
diff --git a/util/precompile_run.jl b/util/precompile_run.jl
index f5a7fb512..43f293ee5 100644
--- a/util/precompile_run.jl
+++ b/util/precompile_run.jl
@@ -66,17 +66,24 @@ for input ∈ [base_input, cheb_input, wall_bc_input, wall_bc_cheb_input]
     push!(inputs_list, x)
 end
 
-collisions_input = merge(wall_bc_cheb_input, Dict("n_neutral_species" => 0,
+collisions_input1 = merge(wall_bc_cheb_input, Dict("n_neutral_species" => 0,
                                                   "krook_collisions" => Dict{String,Any}("use_krook" => true),
                                                   "fokker_planck_collisions" => Dict{String,Any}("use_fokker_planck" => true, "self_collisions" => true, "slowing_down_test" => true),
                                                   "vperp_discretization" => "gausslegendre_pseudospectral",
                                                   "vpa_discretization" => "gausslegendre_pseudospectral",
                                                  ))
+collisions_input2 = merge(wall_bc_cheb_input, Dict("n_neutral_species" => 0,
+                                                  "krook_collisions" => Dict{String,Any}("use_krook" => true),
+                                                  "fokker_planck_collisions" => Dict{String,Any}("use_fokker_planck" => true, "self_collisions" => true, "slowing_down_test" => true),
+                                                  "vperp_discretization" => "gausslegendre_pseudospectral",
+                                                  "vpa_discretization" => "gausslegendre_pseudospectral",
+                                                  "vperp_bc" => "zero-impose-regularity",
+                                                 ))
 # add an additional input for every geometry option available in addition to the default
 geo_input1 = merge(wall_bc_cheb_input, Dict("n_neutral_species" => 0,
                                             "geometry" => Dict{String,Any}("option" => "1D-mirror", "DeltaB" => 0.5, "pitch" => 0.5, "rhostar" => 1.0))) 
 
-push!(inputs_list, collisions_input, geo_input1)
+push!(inputs_list, collisions_input1, collisions_input2, geo_input1)
 
 for input in inputs_list
     run_moment_kinetics(input)