diff --git a/src/charge_exchange.jl b/src/charge_exchange.jl
index f9b27e807..6fa15a22e 100644
--- a/src/charge_exchange.jl
+++ b/src/charge_exchange.jl
@@ -90,7 +90,7 @@ function charge_exchange_collisions_3V!(f_out, f_neutral_out, f_neutral_gav_in,
                     # apply CX collisions to all ion species 
                     # for each ion species, obtain affect of charge exchange collisions
                     # with all of the neutral species
-                    f_out[ivpa,1,iz,ir,is] +=
+                    f_out[ivpa,ivperp,iz,ir,is] +=
                         dt*charge_exchange_frequency*(
                             f_neutral_gav_in[ivpa,ivperp,iz,ir,isn]*fvec_in.density[iz,ir,is]
                             - fvec_in.pdf[ivpa,ivperp,iz,ir,is]*fvec_in.density_neutral[iz,ir,isn])
diff --git a/src/chebyshev.jl b/src/chebyshev.jl
index fc36f5bea..9fef27b4a 100644
--- a/src/chebyshev.jl
+++ b/src/chebyshev.jl
@@ -322,7 +322,7 @@ with all grid points for Clenshaw-Curtis quadrature
 """
 function clenshaw_curtis_weights(ngrid, nelement, n, imin, imax, scale_factor)
     # create array containing the integration weights
-    wgts = zeros(n)
+    wgts = zeros(mk_float, n)
     # calculate the modified Chebshev moments of the first kind
     μ = chebyshevmoments(ngrid)
     @inbounds begin
diff --git a/src/coordinates.jl b/src/coordinates.jl
index 027e2e064..60b70e4ce 100644
--- a/src/coordinates.jl
+++ b/src/coordinates.jl
@@ -7,15 +7,15 @@ export equally_spaced_grid
 
 using ..type_definitions: mk_float, mk_int
 using ..array_allocation: allocate_float, allocate_int
-using ..file_io: open_output_file
-using ..chebyshev: scaled_chebyshev_grid
+using ..calculus: derivative!
+using ..chebyshev: scaled_chebyshev_grid, setup_chebyshev_pseudospectral
 using ..quadrature: composite_simpson_weights
 using ..input_structs: advection_input
 
 """
 structure containing basic information related to coordinates
 """
-struct coordinate
+struct coordinate{Tadvection<:Union{advection_input,Nothing}}
     # name is the name of the variable associated with this coordiante
     name::String
     # n is the total number of grid points associated with this coordinate
@@ -60,7 +60,7 @@ struct coordinate
     # nelement entries
     scratch_2d::Array{mk_float,2}
     # struct containing advection speed options/inputs
-    advection::advection_input
+    advection::Tadvection
 end
 
 """
@@ -95,10 +95,25 @@ function define_coordinate(input, composition=nothing)
     # struct containing the advection speed options/inputs for this coordinate
     advection = input.advection
 
-    return coordinate(input.name, n, input.ngrid, input.nelement, input.L, grid,
+    coord = coordinate(input.name, n, input.ngrid, input.nelement, input.L, grid,
         cell_width, igrid, ielement, imin, imax, input.discretization, input.fd_option,
         input.bc, wgts, uniform_grid, duniform_dgrid, scratch, copy(scratch),
         scratch_2d, advection)
+
+    if input.discretization == "chebyshev_pseudospectral" && coord.ngrid > 1
+        # create arrays needed for explicit Chebyshev pseudospectral treatment in this
+        # coordinate and create the plans for the forward and backward fast Chebyshev
+        # transforms
+        spectral = setup_chebyshev_pseudospectral(coord)
+        # obtain the local derivatives of the uniform grid with respect to the used grid
+        derivative!(coord.duniform_dgrid, coord.uniform_grid, coord, spectral)
+    else
+        # create dummy Bool variable to return in place of the above struct
+        spectral = false
+        coord.duniform_dgrid .= 1.0
+    end
+
+    return coord, spectral
 end
 
 """
diff --git a/src/file_io.jl b/src/file_io.jl
index ca0c08f78..1e88883e4 100644
--- a/src/file_io.jl
+++ b/src/file_io.jl
@@ -10,6 +10,7 @@ export write_data_to_binary
 
 using NCDatasets
 using ..communication: _block_synchronize
+using ..coordinates: coordinate
 using ..looping
 using ..moment_kinetics_structs: scratch_pdf, em_fields_struct
 using ..type_definitions: mk_float, mk_int
@@ -29,56 +30,37 @@ struct ios
     fields::IOStream
 end
 
-# Use this long-winded type (found by using `typeof(v)` where `v` is a variable
-# returned from `NCDatasets.defVar()`) because compiler does not seem to be
-# able to pick up the return types of `defVar()` at compile time, so without
-# using it the result returned from `setup_file_io()` is not a concrete type.
-nc_var_type{N} = Union{
-   NCDatasets.CFVariable{mk_float, N,
-                         NCDatasets.Variable{mk_float, N, NCDatasets.NCDataset},
-                         NCDatasets.Attributes{NCDatasets.NCDataset{Nothing}},
-                         NamedTuple{(:fillvalue, :scale_factor, :add_offset,
-                                     :calendar, :time_origin, :time_factor),
-                                    NTuple{6, Nothing}}},
-   NCDatasets.CFVariable{mk_float, N,
-                         NCDatasets.Variable{mk_float, N,
-                                             NCDatasets.NCDataset{Nothing}},
-                         NCDatasets.Attributes{NCDatasets.NCDataset{Nothing}},
-                         NamedTuple{(:fillvalue, :scale_factor, :add_offset,
-                                     :calendar, :time_origin, :time_factor),
-                                    NTuple{6, Nothing}}}}
-
 """
 structure containing the data/metadata needed for netcdf file i/o
 """
-struct netcdf_info
+struct netcdf_info{Ttime, Tfi, Tphi, Tmom, Tfn}
     # file identifier for the netcdf file to which data is written
     fid::NCDataset
     # handle for the time variable
-    time::nc_var_type{1}
+    time::Ttime
     # handle for the charged species distribution function variable
-    f::nc_var_type{6}
+    f::Tfi
     # handle for the electrostatic potential variable
-    phi::nc_var_type{3}
+    phi::Tphi
     # handle for the charged species density
-    density::nc_var_type{4}
+    density::Tmom
     # handle for the charged species parallel flow
-    parallel_flow::nc_var_type{4}
+    parallel_flow::Tmom
     # handle for the charged species parallel pressure
-    parallel_pressure::nc_var_type{4}
+    parallel_pressure::Tmom
     # handle for the charged species parallel heat flux
-    parallel_heat_flux::nc_var_type{4}
+    parallel_heat_flux::Tmom
     # handle for the charged species thermal speed
-    thermal_speed::nc_var_type{4}
+    thermal_speed::Tmom
     
     # handle for the neutral species distribution function variable
-    f_neutral::nc_var_type{7}
+    f_neutral::Tfn
     # handle for the neutral species density
-    density_neutral::nc_var_type{4}
-    uz_neutral::nc_var_type{4}
-    pz_neutral::nc_var_type{4}
-    qz_neutral::nc_var_type{4}
-    thermal_speed_neutral::nc_var_type{4}
+    density_neutral::Tmom
+    uz_neutral::Tmom
+    pz_neutral::Tmom
+    qz_neutral::Tmom
+    thermal_speed_neutral::Tmom
 
 end
 
@@ -142,6 +124,10 @@ function define_dimensions!(fid, nvz, nvr, nvzeta, nvpa, nvperp, nz, nr, n_speci
     end
     # define the time dimension, with an expandable size (denoted by Inf)
     defDim(fid, "ntime", Inf)
+    # define a length-1 dimension for storing strings. Don't know why they cannot be
+    # stored as scalars, maybe did not find the right function/method, maybe a missing
+    # feature or bug in NCDatasets.jl?
+    defDim(fid, "str_dim", 1)
 
     return nothing
 end
@@ -152,97 +138,71 @@ end
 Define static (i.e. time-independent) variables for an output file.
 """
 function define_static_variables!(fid,vz,vr,vzeta,vpa,vperp,z,r,composition,collisions)
-    # create and write the "r" variable to file
-    varname = "r"
-    attributes = Dict("description" => "radial coordinate")
-    dims = ("nr",)
-    vartype = mk_float
-    var = defVar(fid, varname, vartype, dims, attrib=attributes)
-    var[:] = r.grid
-    # create and write the "r_wgts" variable to file
-    varname = "r_wgts"
-    attributes = Dict("description" => "integration weights for radial coordinate")
-    vartype = mk_float
-    var = defVar(fid, varname, vartype, dims, attrib=attributes)
-    var[:] = r.wgts
-    # create and write the "z" variable to file
-    varname = "z"
-    attributes = Dict("description" => "parallel coordinate")
-    dims = ("nz",)
-    vartype = mk_float
-    var = defVar(fid, varname, vartype, dims, attrib=attributes)
-    var[:] = z.grid
-    # create and write the "z_wgts" variable to file
-    varname = "z_wgts"
-    attributes = Dict("description" => "integration weights for parallel coordinate")
-    vartype = mk_float
-    var = defVar(fid, varname, vartype, dims, attrib=attributes)
-    var[:] = z.wgts
-    # create and write the "vperp" variable to file
-    varname = "vperp"
-    attributes = Dict("description" => "parallel velocity")
-    dims = ("nvperp",)
-    vartype = mk_float
-    var = defVar(fid, varname, vartype, dims, attrib=attributes)
-    var[:] = vperp.grid
-    # create and write the "vperp_wgts" variable to file
-    varname = "vperp_wgts"
-    attributes = Dict("description" => "integration weights for parallel velocity coordinate")
-    vartype = mk_float
-    var = defVar(fid, varname, vartype, dims, attrib=attributes)
-    var[:] = vperp.wgts
-    # create and write the "vpa" variable to file
-    varname = "vpa"
-    attributes = Dict("description" => "parallel velocity")
-    dims = ("nvpa",)
-    vartype = mk_float
-    var = defVar(fid, varname, vartype, dims, attrib=attributes)
-    var[:] = vpa.grid
-    # create and write the "vpa_wgts" variable to file
-    varname = "vpa_wgts"
-    attributes = Dict("description" => "integration weights for parallel velocity coordinate")
-    vartype = mk_float
-    var = defVar(fid, varname, vartype, dims, attrib=attributes)
-    var[:] = vpa.wgts
-    # create and write the "vzeta" variable to file
-    varname = "vzeta"
-    attributes = Dict("description" => "parallel velocity")
-    dims = ("nvzeta",)
-    vartype = mk_float
-    var = defVar(fid, varname, vartype, dims, attrib=attributes)
-    var[:] = vzeta.grid
-    # create and write the "vzeta_wgts" variable to file
-    varname = "vzeta_wgts"
-    attributes = Dict("description" => "integration weights for parallel velocity coordinate")
-    vartype = mk_float
-    var = defVar(fid, varname, vartype, dims, attrib=attributes)
-    var[:] = vzeta.wgts
-    # create and write the "vr" variable to file
-    varname = "vr"
-    attributes = Dict("description" => "parallel velocity")
-    dims = ("nvr",)
-    vartype = mk_float
-    var = defVar(fid, varname, vartype, dims, attrib=attributes)
-    var[:] = vr.grid
-    # create and write the "vr_wgts" variable to file
-    varname = "vr_wgts"
-    attributes = Dict("description" => "integration weights for parallel velocity coordinate")
-    vartype = mk_float
-    var = defVar(fid, varname, vartype, dims, attrib=attributes)
-    var[:] = vr.wgts
-    # create and write the "vz" variable to file
-    varname = "vz"
-    attributes = Dict("description" => "parallel velocity")
-    dims = ("nvz",)
-    vartype = mk_float
-    var = defVar(fid, varname, vartype, dims, attrib=attributes)
-    var[:] = vz.grid
-    # create and write the "vz_wgts" variable to file
-    varname = "vz_wgts"
-    attributes = Dict("description" => "integration weights for parallel velocity coordinate")
-    vartype = mk_float
-    var = defVar(fid, varname, vartype, dims, attrib=attributes)
-    var[:] = vz.wgts
+    function save_coordinate(coord::coordinate, description::String)
+        # Create and write the grid for coord
+        varname = coord.name
+        attributes = Dict("description" => description)
+        dims = ("n$(coord.name)",)
+        vartype = mk_float
+        var = defVar(fid, varname, vartype, dims, attrib=attributes)
+        var[:] = coord.grid
+
+        # create and write the weights for coord
+        varname = "$(coord.name)_wgts"
+        attributes = Dict("description" => "integration weights for $(coord.name) coordinate")
+        vartype = mk_float
+        var = defVar(fid, varname, vartype, dims, attrib=attributes)
+        var[:] = coord.wgts
+
+        # create and write ngrid for coord
+        varname = "$(coord.name)_ngrid"
+        attributes = Dict("description" => "ngrid for $(coord.name) coordinate")
+        dims = ()
+        vartype = mk_int
+        var = defVar(fid, varname, vartype, dims, attrib=attributes)
+        var[:] = coord.ngrid
+
+        # create and write nelement for coord
+        varname = "$(coord.name)_nelement"
+        attributes = Dict("description" => "nelement for $(coord.name) coordinate")
+        dims = ()
+        vartype = mk_int
+        var = defVar(fid, varname, vartype, dims, attrib=attributes)
+        var[:] = coord.nelement
+
+        # create and write discretization for coord
+        varname = "$(coord.name)_discretization"
+        attributes = Dict("description" => "discretization for $(coord.name) coordinate")
+        dims = ("str_dim",)
+        vartype = String
+        var = defVar(fid, varname, vartype, dims, attrib=attributes)
+        var[:] = coord.discretization
+
+        # create and write fd_option for coord
+        varname = "$(coord.name)_fd_option"
+        attributes = Dict("description" => "fd_option for $(coord.name) coordinate")
+        dims = ("str_dim",)
+        vartype = String
+        var = defVar(fid, varname, vartype, dims, attrib=attributes)
+        var[:] = coord.fd_option
+
+        # create and write bc for coord
+        varname = "$(coord.name)_bc"
+        attributes = Dict("description" => "bc for $(coord.name) coordinate")
+        dims = ("str_dim",)
+        vartype = String
+        var = defVar(fid, varname, vartype, dims, attrib=attributes)
+        var[:] = coord.bc
+    end
+
+    # create and write the coordinate variables
+    save_coordinate(r, "radial coordinate")
+    save_coordinate(z, "parallel coordinate")
+    save_coordinate(vperp, "perpendicular velocity")
+    save_coordinate(vpa, "parallel velocity")
+    save_coordinate(vzeta, "neutral binormal velocity")
+    save_coordinate(vr, "neutral radial velocity")
+    save_coordinate(vz, "neutral parallel velocity")
     # create and write the "T_e" variable to file
     varname = "T_e"
     attributes = Dict("description" => "electron temperature")
diff --git a/src/initial_conditions.jl b/src/initial_conditions.jl
index 94a53eb4c..2e7927c06 100644
--- a/src/initial_conditions.jl
+++ b/src/initial_conditions.jl
@@ -184,24 +184,18 @@ function init_pdf!(pdf, moments, vz, vr, vzeta, vpa, vperp, z, r, n_ion_species,
 end
 
 function init_pdf_moments_manufactured_solns!(pdf, moments, vz, vr, vzeta, vpa, vperp, z, r, n_ion_species, n_neutral_species, geometry,composition)
-    manufactured_solns_list = manufactured_solutions(r.L,z.L,r.bc,z.bc,geometry,composition,r.n) 
+    manufactured_solns_list = manufactured_solutions(r.L,z.L,vperp.L,vpa.L,vzeta.L,vr.L,vz.L,r.bc,z.bc,geometry,composition,r.n) 
     dfni_func = manufactured_solns_list.dfni_func
     densi_func = manufactured_solns_list.densi_func
     dfnn_func = manufactured_solns_list.dfnn_func
     densn_func = manufactured_solns_list.densn_func
     #nb manufactured functions not functions of species
-    for is in 1:n_ion_species
-        for ir in 1:r.n
-            for iz in 1:z.n
-                moments.charged.dens[iz,ir,is] = densi_func(z.grid[iz],r.grid[ir],0.0)
-                for ivperp in 1:vperp.n
-                    for ivpa in 1:vpa.n
-                        pdf.charged.unnorm[ivpa,ivperp,iz,ir,is] = dfni_func(vpa.grid[ivpa],vperp.grid[ivperp],z.grid[iz],r.grid[ir],0.0)
-                        pdf.charged.norm[ivpa,ivperp,iz,ir,is] = pdf.charged.unnorm[ivpa,ivperp,iz,ir,is]
-                        
-                    end
-                end
-            end
+    begin_s_r_z_region()
+    @loop_s_r_z is ir iz begin
+        moments.charged.dens[iz,ir,is] = densi_func(z.grid[iz],r.grid[ir],0.0)
+        @loop_vperp_vpa ivperp ivpa begin
+            pdf.charged.unnorm[ivpa,ivperp,iz,ir,is] = dfni_func(vpa.grid[ivpa],vperp.grid[ivperp],z.grid[iz],r.grid[ir],0.0)
+            pdf.charged.norm[ivpa,ivperp,iz,ir,is] = pdf.charged.unnorm[ivpa,ivperp,iz,ir,is]
         end
     end
     # update upar, ppar, qpar, vth consistent with manufactured solns
@@ -210,6 +204,7 @@ function init_pdf_moments_manufactured_solns!(pdf, moments, vz, vr, vzeta, vpa,
     # get particle flux
     update_upar!(moments.charged.upar, pdf.charged.unnorm, vpa, vperp, z, r, composition)
     # convert from particle particle flux to parallel flow
+    begin_s_r_z_region()
     @loop_s_r_z is ir iz begin
         moments.charged.upar[iz,ir,is] /= moments.charged.dens[iz,ir,is]
     # update the thermal speed
@@ -217,20 +212,12 @@ function init_pdf_moments_manufactured_solns!(pdf, moments, vz, vr, vzeta, vpa,
     end
     
     if n_neutral_species > 0
-        for isn in 1:n_neutral_species
-            for ir in 1:r.n
-                for iz in 1:z.n
-                    moments.neutral.dens[iz,ir,isn] = densn_func(z.grid[iz],r.grid[ir],0.0)
-                    for ivzeta in 1:vzeta.n
-                        for ivr in 1:vr.n
-                            for ivz in 1:vz.n
-                                pdf.neutral.unnorm[ivz,ivr,ivzeta,iz,ir,isn] = dfnn_func(vz.grid[ivz],vr.grid[ivr],vzeta.grid[ivzeta],z.grid[iz],r.grid[ir],0.0)
-                                pdf.neutral.norm[ivz,ivr,ivzeta,iz,ir,isn] = pdf.neutral.unnorm[ivz,ivr,ivzeta,iz,ir,isn]
-                            end
-                        end
-                    end
-                    
-                end
+        begin_sn_r_z_region()
+        @loop_sn_r_z isn ir iz begin
+            moments.neutral.dens[iz,ir,isn] = densn_func(z.grid[iz],r.grid[ir],0.0)
+            @loop_vzeta_vr_vz ivzeta ivr ivz begin
+                pdf.neutral.unnorm[ivz,ivr,ivzeta,iz,ir,isn] = dfnn_func(vz.grid[ivz],vr.grid[ivr],vzeta.grid[ivzeta],z.grid[iz],r.grid[ir],0.0)
+                pdf.neutral.norm[ivz,ivr,ivzeta,iz,ir,isn] = pdf.neutral.unnorm[ivz,ivr,ivzeta,iz,ir,isn]
             end
         end
         # get consistent moments with manufactured solutions 
@@ -241,8 +228,9 @@ function init_pdf_moments_manufactured_solns!(pdf, moments, vz, vr, vzeta, vpa,
         update_neutral_pzeta!(moments.neutral.pzeta, pdf.neutral.unnorm, vz, vr, vzeta, z, r, composition)
         #update ptot (isotropic pressure)
         if r.n > 1 #if 2D geometry
+            begin_sn_r_z_region()
             @loop_sn_r_z isn ir iz begin            
-            moments.neutral.ptot[iz,ir,isn] = (moments.neutral.pz[iz,ir,isn] + moments.neutral.pr[iz,ir,isn] + moments.neutral.pzeta[iz,ir,isn])/3.0
+                moments.neutral.ptot[iz,ir,isn] = (moments.neutral.pz[iz,ir,isn] + moments.neutral.pr[iz,ir,isn] + moments.neutral.pzeta[iz,ir,isn])/3.0
             end
         else #1D model
             moments.neutral.ptot .= moments.neutral.pz
@@ -252,6 +240,7 @@ function init_pdf_moments_manufactured_solns!(pdf, moments, vz, vr, vzeta, vpa,
         update_neutral_ur!(moments.neutral.ur, pdf.neutral.unnorm, vz, vr, vzeta, z, r, composition)
         update_neutral_uzeta!(moments.neutral.uzeta, pdf.neutral.unnorm, vz, vr, vzeta, z, r, composition)
         # now convert from particle particle flux to parallel flow
+        begin_sn_r_z_region()
         @loop_sn_r_z isn ir iz begin
             moments.neutral.uz[iz,ir,isn] /= moments.neutral.dens[iz,ir,isn]
             moments.neutral.ur[iz,ir,isn] /= moments.neutral.dens[iz,ir,isn]
diff --git a/src/input_structs.jl b/src/input_structs.jl
index f17d5f34b..eb8cdff70 100644
--- a/src/input_structs.jl
+++ b/src/input_structs.jl
@@ -38,6 +38,26 @@ struct time_input
     use_manufactured_solns::Bool
 end
 
+"""
+"""
+mutable struct advance_info
+    vpa_advection::Bool
+    z_advection::Bool
+    r_advection::Bool
+    neutral_z_advection::Bool
+    neutral_r_advection::Bool
+    cx_collisions::Bool
+    cx_collisions_1V::Bool
+    ionization_collisions::Bool
+    ionization_collisions_1V::Bool
+    source_terms::Bool
+    continuity::Bool
+    force_balance::Bool
+    energy::Bool
+    rk_coefs::Array{mk_float,2}
+    manufactured_solns_test::Bool
+end
+
 """
 """
 mutable struct advection_input_mutable
@@ -92,7 +112,7 @@ end
 
 """
 """
-struct grid_input
+struct grid_input{Tadvection<:Union{advection_input,Nothing}}
     # name of the variable associated with this coordinate
     name::String
     # number of grid points per element
@@ -108,7 +128,7 @@ struct grid_input
     # boundary option
     bc::String
     # struct containing advection speed options
-    advection::advection_input
+    advection::Tadvection
 end
 
 """
diff --git a/src/load_data.jl b/src/load_data.jl
index 5f38b9000..7bb71e9d1 100644
--- a/src/load_data.jl
+++ b/src/load_data.jl
@@ -11,6 +11,9 @@ export load_pdf_data
 export load_neutral_pdf_data
 export load_neutral_coordinate_data
 
+using ..coordinates: define_coordinate
+using ..input_structs: grid_input
+
 using NCDatasets
 
 """
@@ -35,49 +38,50 @@ end
 """
 function load_coordinate_data(fid)
     print("Loading coordinate data...")
-    # define a handle for the r coordinate
-    cdfvar = fid["r"]
-    # get the number of r grid points
-    nr = length(cdfvar)
-    # load the data for r
-    r = cdfvar.var[:]
-    # get the weights associated with the r coordinate
-    cdfvar = fid["r_wgts"]
-    r_wgts = cdfvar.var[:]
-    # Lr = r box length
-    Lr = r[end]-r[1]
-    
-    # define a handle for the z coordinate
-    cdfvar = fid["z"]
-    # get the number of z grid points
-    nz = length(cdfvar)
-    # load the data for z
-    z = cdfvar.var[:]
-    # get the weights associated with the z coordinate
-    cdfvar = fid["z_wgts"]
-    z_wgts = cdfvar.var[:]
-    # Lz = z box length
-    Lz = z[end]-z[1]
-
-    # define a handle for the vperp coordinate
-    cdfvar = fid["vperp"]
-    # get the number of vperp grid points
-    nvperp = length(cdfvar)
-    # load the data for vperp
-    vperp = cdfvar.var[:]
-    # get the weights associated with the vperp coordinate
-    cdfvar = fid["vperp_wgts"]
-    vperp_wgts = cdfvar.var[:]
-
-    # define a handle for the vpa coordinate
-    cdfvar = fid["vpa"]
-    # get the number of vpa grid points
-    nvpa = length(cdfvar)
-    # load the data for vpa
-    vpa = cdfvar.var[:]
-    # get the weights associated with the vpa coordinate
-    cdfvar = fid["vpa_wgts"]
-    vpa_wgts = cdfvar.var[:]
+
+    function load_coordinate(name::String)
+        # define a handle for the coordinate
+        cdfvar = fid[name]
+        # get the number of grid points
+        n = length(cdfvar)
+        # load the data for grid point positions
+        grid = cdfvar.var[:]
+        # get the weights associated with the coordinate
+        cdfvar = fid["$(name)_wgts"]
+        wgts = cdfvar.var[:]
+        # Lr = box length for coordinate
+        L = grid[end]-grid[1]
+
+        cdfvar = fid["$(name)_ngrid"]
+        ngrid = cdfvar.var[:]
+
+        cdfvar = fid["$(name)_nelement"]
+        nelement = cdfvar.var[:]
+
+        cdfvar = fid["$(name)_discretization"]
+        discretization = cdfvar.var[1]
+
+        cdfvar = fid["$(name)_fd_option"]
+        fd_option = cdfvar.var[1]
+
+        cdfvar = fid["$(name)_bc"]
+        bc = cdfvar.var[1]
+
+        input = grid_input(name, ngrid, nelement, L, discretization, fd_option, bc, nothing)
+
+        coord, spectral = define_coordinate(input)
+
+        # grid is recreated in define_coordinate, so check it is consistent with the
+        # saved grid positions
+        @assert isapprox(grid, coord.grid, rtol=1.e-14)
+
+        return coord
+    end
+
+    r = load_coordinate("r")
+    z = load_coordinate("z")
+    vperp = load_coordinate("vperp")
+    vpa = load_coordinate("vpa")
 
     # define a handle for the time coordinate
     cdfvar = fid["time"]
@@ -92,7 +96,7 @@ function load_coordinate_data(fid)
     n_neutral_species = fid.dim["n_neutral_species"]
     println("done.")
 
-    return nvpa, vpa, vpa_wgts, nvperp, vperp, vperp_wgts, nz, z, z_wgts, Lz, nr, r, r_wgts, Lr, ntime, time, n_ion_species, n_neutral_species
+    return vpa, vperp, z, r, ntime, time
 end
 
 function load_neutral_coordinate_data(fid)
diff --git a/src/manufactured_solns.jl b/src/manufactured_solns.jl
index 59fc7e820..3ccb7d351 100644
--- a/src/manufactured_solns.jl
+++ b/src/manufactured_solns.jl
@@ -4,8 +4,16 @@ module manufactured_solns
 
 export manufactured_solutions
 export manufactured_sources
+export manufactured_solutions_as_arrays
+export manufactured_rhs_as_array
 
+using SpecialFunctions
 using Symbolics
+using ..array_allocation: allocate_shared_float
+using ..coordinates: coordinate
+using ..input_structs: geometry_input, advance_info, species_composition, collisions_input
+using ..looping
+using ..type_definitions
 
     @variables r z vpa vperp t vz vr vzeta
 
@@ -44,7 +52,7 @@ using Symbolics
     # neutral density symbolic function
     function densn_sym(Lr,Lz,r_bc,z_bc,geometry,composition)
         if r_bc == "periodic" && z_bc == "periodic" 
-            densn = 1.5 +  0.1*(cos(2.0*pi*r/Lr) + cos(2.0*pi*z/Lz))*sin(2.0*pi*t)  
+            densn = 1.5 +  0.1*(cos(2.0*pi*r/Lr) + cos(2.0*pi*z/Lz))*cos(2.0*pi*t)
         elseif (r_bc == "periodic" && z_bc == "wall")
             T_wall = composition.T_wall
             Bzed = geometry.Bzed
@@ -60,40 +68,73 @@ using Symbolics
     end
 
     # neutral distribution symbolic function
-    function dfnn_sym(Lr,Lz,r_bc,z_bc,geometry,composition)
+    function dfnn_sym(Lr,Lz,Lvzeta,Lvr,Lvz,r_bc,z_bc,geometry,composition)
         densn = densn_sym(Lr,Lz,r_bc,z_bc,geometry,composition)
         if (r_bc == "periodic" && z_bc == "periodic")
-            dfnn = densn * exp( - vz^2 - vr^2 - vzeta^2)
+            # Ad-hoc odd-in-v* component of dfnn which gives non-zero u* and dfnn
+            # positive everywhere, while the odd components integrate to 0 so do not
+            # need to be accounted for in normalisation.
+            dfnn = densn * (exp( - vz^2 - vr^2 - vzeta^2)
+                            + (sin(two*pi*r/Lr) + sin(2*pi*z/Lz))/ten
+                              * (vz + vr + vzeta) * exp(-two*(vz^2+vr^2+vzeta^2)))
+            # Force the symbolic function dfnn to vanish when v*=±Lv*/2, so that it is
+            # consistent with "zero" bc in v*.
+            # Normalisation factors ensure that when this f is
+            # integrated on -Lv*/2≤v*≤Lv*/2 the result (taking into account the
+            # normalisations used by moment_kinetics) is exactly n.
+            # Note that:
+            #  ∫_-Lv/2^Lv/2  [1-(2*v/Lv)^2]exp(-v^2) dv
+            #  = [sqrt(π)(Lv^2-2)*erf(Lv/2) + 2*exp(-Lv^2/4)Lv] / Lv^2
+            dfnn = (one - (two*vzeta/Lvzeta)^2) * dfnn *
+                      Lvzeta / ((Lvzeta^2-2)*erf(Lvzeta/two) + two/sqrt(π)*exp(-Lvzeta^2/four)*Lvzeta) # vzeta normalization
+            dfnn = (one - (two*vr/Lvr)^2) * dfnn *
+                      Lvr / ((Lvr^2-2)*erf(Lvr/two) + two/sqrt(π)*exp(-Lvr^2/four)*Lvr) # vr normalization
+            dfnn = (one - (two*vz/Lvz)^2) * dfnn *
+                      Lvz / ((Lvz^2-2)*erf(Lvz/two) + two/sqrt(π)*exp(-Lvz^2/four)*Lvz) # vz normalization
         elseif (r_bc == "periodic" && z_bc == "wall")
-            Hplus = 0.5*(sign(vz) + 1.0)
-            Hminus = 0.5*(sign(-vz) + 1.0)
+            Hplus = half*(sign(vz) + one)
+            Hminus = half*(sign(-vz) + one)
             FKw = knudsen_cosine(composition)
             Bzed = geometry.Bzed
             Bmag = geometry.Bmag
-            Gamma_minus = 0.5*(Bzed/Bmag)*nminus_sym(Lr,r_bc)/sqrt(pi)
-            Gamma_plus = 0.5*(Bzed/Bmag)*nplus_sym(Lr,r_bc)/sqrt(pi)
-            dfnn = Hplus *( ( 0.5 - z/Lz)*Gamma_minus + 1.0 )*FKw + Hminus*( ( 0.5 + z/Lz)*Gamma_plus + 1.0 )*FKw 
+            Gamma_minus = half*(Bzed/Bmag)*nminus_sym(Lr,r_bc)/sqrt(pi)
+            Gamma_plus = half*(Bzed/Bmag)*nplus_sym(Lr,r_bc)/sqrt(pi)
+            dfnn = Hplus *( ( half - z/Lz)*Gamma_minus + one )*FKw + Hminus*( ( half + z/Lz)*Gamma_plus + one )*FKw 
         end
         return dfnn
     end
-    function gyroaveraged_dfnn_sym(Lr,Lz,r_bc,z_bc,geometry,composition)
-        densn = densn_sym(Lr,Lz,r_bc,z_bc,geometry,composition)
-        #if (r_bc == "periodic" && z_bc == "periodic")
-            dfnn = densn * exp( - vpa^2 - vperp^2 )
-        #end
-        return dfnn
+    function gyroaveraged_dfnn_sym(Lr,Lz,Lvzeta,Lvr,Lvz,r_bc,z_bc,geometry,composition)
+        dfnn = dfnn_sym(Lr,Lz,Lvzeta,Lvr,Lvz,r_bc,z_bc,geometry,composition)
+        # Because of particular form of odd component in dfnn_sym, if we set vzeta=-vr
+        # then we eliminate the odd component in the perpendicular velocities,
+        # effectively gyro-averaging
+        gav_dfnn = substitute(dfnn, Dict(vzeta=>vperp/sqrt(two), vr=>-vperp/sqrt(two), vz=>vpa))
+        return gav_dfnn
     end
     
-    # ion density symbolic function
+    # Define some constants in mk_float. Avoids loss of precision due to implicit
+    # conversions from FLoat64->Float128 if we use Float128.
+    const half = 1/mk_float(2)
+    const one = mk_float(1)
+    const three_over_two = mk_float(3)/mk_float(2)
+    const two = mk_float(2)
+    const four = mk_float(4)
+    const five = mk_float(5)
+    const ten = mk_float(10)
+
     function densi_sym(Lr,Lz,r_bc,z_bc)
+        # Note: explicitly convert numerical factors to mk_float so the output gets full
+        # precision if we use quad-precision (Float128)
         if r_bc == "periodic" && z_bc == "periodic"
-            densi = 1.5 +  0.1*(sin(2.0*pi*r/Lr) + sin(2.0*pi*z/Lz))*sin(2.0*pi*t)  
+            densi = three_over_two +  (sin(two*pi*r/Lr) + sin(two*pi*z/Lz))/ten*cos(two*pi*t)
         elseif r_bc == "Dirichlet" && z_bc == "periodic"
-            #densi = 1.0 +  0.5*sin(2.0*pi*z/Lz)*(r/Lr + 0.5) + 0.2*sin(2.0*pi*r/Lr)*sin(2.0*pi*t)
-            #densi = 1.0 +  0.5*sin(2.0*pi*z/Lz)*(r/Lr + 0.5) + sin(2.0*pi*r/Lr)*sin(2.0*pi*t)
-            densi = 1.0 +  0.5*(r/Lr + 0.5) 
+            #densi = 1 +  1//2*sin(2*pi*z/Lz)*(r/Lr + 1//2) + 0.2*sin(2*pi*r/Lr)*sin(2*pi*t)
+            #densi = 1 +  1//2*sin(2*pi*z/Lz)*(r/Lr + 1//2) + sin(2*pi*r/Lr)*sin(2*pi*t)
+            densi = one +  half*(r/Lr + half)/two
         elseif r_bc == "periodic" && z_bc == "wall"
-            densi = 0.25*(0.5 - z/Lz)*nminus_sym(Lr,r_bc) + 0.25*(z/Lz + 0.5)*nplus_sym(Lr,r_bc) + (z/Lz + 0.5)*(0.5 - z/Lz)*nzero_sym(Lr,r_bc)  #+  0.5*(r/Lr + 0.5) + 0.5*(z/Lz + 0.5)
+            densi = (half - z/Lz)/four + (z/Lz + half)/four + (z/Lz + half)*(half - z/Lz)/five  #+  1//2*(r/Lr + 1//2) + 1//2*(z/Lz + 1//2)
+        else
+            error("Unsupported options r_bc=$r_bc, z_bc=$z_bc")
         end
         return densi
     end
@@ -101,7 +142,9 @@ using Symbolics
     
 
     # ion distribution symbolic function
-    function dfni_sym(Lr,Lz,r_bc,z_bc,geometry,nr)
+    function dfni_sym(Lr,Lz,Lvperp,Lvpa,r_bc,z_bc,geometry,nr)
+        # Note: explicitly convert numerical factors to mk_float so the output gets full
+        # precision if we use quad-precision (Float128)
         densi = densi_sym(Lr,Lz,r_bc,z_bc)
         # define derivative operators
         Dr = Differential(r) 
@@ -122,31 +165,68 @@ using Symbolics
         end
         
         if (r_bc == "periodic" && z_bc == "periodic") || (r_bc == "Dirichlet" && z_bc == "periodic")
-            dfni = densi * exp( - vpa^2 - vperp^2) 
+            ## This version with upar is very expensive to evaluate, probably due to
+            ## spatially-varying erf()?
+            #upar = (sin(two*pi*r/Lr) + sin(2*pi*z/Lz))/ten
+            #dfni = densi * exp(- (vpa-upar)^2 - vperp^2) #/ pi^1.5
+            ## Force the symbolic function dfni to vanish when vpa=±Lvpa/2, so that it is
+            ## consistent with "zero" bc in vpa.
+            ## Normalisation factors on 2nd and 3rd lines ensure that when this f is
+            ## integrated on -Lvpa/2≤vpa≤Lvpa/2 and 0≤vperp≤Lvperp the result (taking
+            ## into account the normalisations used by moment_kinetics) is exactly n.
+            ## Note that:
+            ##  ∫_-Lvpa/2^Lvpa/2  [1-(2*(vpa-upar)/Lvpa)^2]exp(-(vpa-upar)^2) dvpa
+            ##  = [sqrt(π)(Lvpa^2-4*upar^2-2)*(erf(Lvpa/2-upar)+erf(Lvpa/2+upar))
+            ##     +2*exp(-(Lvpa+2upar)^2/4)*(exp(2*Lvpa*upar)*(Lvpa+2upar)+Lvpa-2*upar)] / 2 / Lvpa^2
+            ##
+            ##  ∫_0^Lvperp vperp*exp(-vperp^2) = (1 - exp(-Lvperp^2))/2
+            #dfni = (one - (two*vpa/Lvpa)^2) * dfni *
+            #          two*Lvpa / ((Lvpa^2-four*upar^2-2)*(erf(Lvpa/two-upar)+erf(Lvpa/two+upar)) + two/sqrt(π)*exp(-(Lvpa+two*upar)^2/four)*(exp(two*Lvpa*upar)+Lvpa-two*upar)) / # vpa normalization
+            #         (one - exp(-Lvperp^2)) # vperp normalisation
+
+            # Ad-hoc odd-in-vpa component of dfni which gives non-zero upar and dfni
+            # positive everywhere, while the odd component integrates to 0 so does not
+            # need to be accounted for in normalisation.
+            dfni = densi * (exp(- vpa^2 - vperp^2)
+                            + (sin(two*pi*r/Lr) + sin(2*pi*z/Lz))/ten
+                              * vpa * exp(-two*(vpa^2+vperp^2))) #/ pi^1.5
+            # Force the symbolic function dfni to vanish when vpa=±Lvpa/2, so that it is
+            # consistent with "zero" bc in vpa.
+            # Normalisation factors on 2nd and 3rd lines ensure that when this f is
+            # integrated on -Lvpa/2≤vpa≤Lvpa/2 and 0≤vperp≤Lvperp the result (taking
+            # into account the normalisations used by moment_kinetics) is exactly n.
+            # Note that:
+            #  ∫_-Lvpa/2^Lvpa/2  [1-(2*vpa/Lvpa)^2]exp(-vpa^2) dvpa
+            #  = [sqrt(π)(Lvpa^2-2)*erf(Lvpa/2) + 2*exp(-Lvpa^2/4)Lvpa] / Lvpa^2
+            #
+            #  ∫_0^Lvperp vperp*exp(-vperp^2) = (1 - exp(-Lvperp^2))/2
+            dfni = (one - (two*vpa/Lvpa)^2) * dfni *
+                      Lvpa / ((Lvpa^2-2)*erf(Lvpa/two) + two/sqrt(π)*exp(-Lvpa^2/four)*Lvpa) / # vpa normalization
+                     (one - exp(-Lvperp^2)) # vperp normalisation
         elseif r_bc == "periodic" && z_bc == "wall"
-            vpabar = vpa - (rhostar/2.0)*(Bmag/Bzed)*expand_derivatives(Er)*rfac # effective velocity in z direction * (Bmag/Bzed)
-            Hplus = 0.5*(sign(vpabar) + 1.0)
-            Hminus = 0.5*(sign(-vpabar) + 1.0)
+            vpabar = vpa - (rhostar/two)*(Bmag/Bzed)*expand_derivatives(Er)*rfac # effective velocity in z direction * (Bmag/Bzed)
+            Hplus = half*(sign(vpabar) + one)
+            Hminus = half*(sign(-vpabar) + one)
             ffa =  exp(- vperp^2)
-            dfni = ffa * ( nminus_sym(Lr,r_bc)* (0.5 - z/Lz) * Hminus * vpabar^2 + nminus_sym(Lr,r_bc)*(z/Lz + 0.5) * Hplus * vpabar^2 + nzero_sym(Lr,r_bc)*(z/Lz + 0.5)*(0.5 - z/Lz) ) * exp( - vpabar^2 )
+            dfni = ffa * ( nminus_sym(Lr,r_bc)* (half - z/Lz) * Hminus * vpabar^2 + nminus_sym(Lr,r_bc)*(z/Lz + half) * Hplus * vpabar^2 + nzero_sym(Lr,r_bc)*(z/Lz + half)*(half - z/Lz) ) * exp( - vpabar^2 )
+        else
+            error("Unsupported options r_bc=$r_bc, z_bc=$z_bc")
         end
         return dfni
     end
-    function cartesian_dfni_sym(Lr,Lz,r_bc,z_bc)
-        densi = densi_sym(Lr,Lz,r_bc,z_bc)
-        #if (r_bc == "periodic" && z_bc == "periodic") || (r_bc == "Dirichlet" && z_bc == "periodic")
-            dfni = densi * exp( - vz^2 - vr^2 - vzeta^2) 
-        #end
-        return dfni
+    function cartesian_dfni_sym(Lr,Lz,Lvperp,Lvpa,r_bc,z_bc,geometry,nr)
+        dfni = dfni_sym(Lr,Lz,Lvperp,Lvpa,r_bc,z_bc,geometry,nr)
+        vrvzvzeta_dfni = substitute(dfni, Dict(vperp=>sqrt(vzeta^2+vr^2), vpa=>vz))
+        return vrvzvzeta_dfni
     end
 
-    function manufactured_solutions(Lr,Lz,r_bc,z_bc,geometry,composition,nr)
+    function manufactured_solutions(Lr,Lz,Lvperp,Lvpa,Lvzeta,Lvr,Lvz,r_bc,z_bc,geometry,composition,nr)
 
         densi = densi_sym(Lr,Lz,r_bc,z_bc)
-        dfni = dfni_sym(Lr,Lz,r_bc,z_bc,geometry,nr)
+        dfni = dfni_sym(Lr,Lz,Lvperp,Lvpa,r_bc,z_bc,geometry,nr)
         
         densn = densn_sym(Lr,Lz,r_bc,z_bc,geometry,composition)
-        dfnn = dfnn_sym(Lr,Lz,r_bc,z_bc,geometry,composition)
+        dfnn = dfnn_sym(Lr,Lz,Lvzeta,Lvr,Lvz,r_bc,z_bc,geometry,composition)
         
         #build julia functions from these symbolic expressions
         # cf. https://docs.juliahub.com/Symbolics/eABRO/3.4.0/tutorials/symbolic_functions/
@@ -166,24 +246,31 @@ using Symbolics
         return manufactured_solns_list
     end 
 
-    function manufactured_sources(Lr,Lz,r_bc,z_bc,composition,geometry,collisions,nr)
-        
+    function manufactured_rhs_sym(Lr::mk_float, Lz::mk_float, Lvperp::mk_float,
+                                  Lvpa::mk_float,Lvzeta::mk_float, Lvr::mk_float,
+                                  Lvz::mk_float, r_bc::String, z_bc::String,
+                                  composition::species_composition,
+                                  geometry::geometry_input,
+                                  collisions::collisions_input, nr::mk_int,
+                                  advance::Union{advance_info,Nothing}=nothing)
+        # Note: explicitly convert numerical factors to mk_float so the output gets full
+        # precision if we use quad-precision (Float128)
+
         # ion manufactured solutions
         densi = densi_sym(Lr,Lz,r_bc,z_bc)
-        dfni = dfni_sym(Lr,Lz,r_bc,z_bc,geometry,nr)
-        vrvzvzeta_dfni = cartesian_dfni_sym(Lr,Lz,r_bc,z_bc) #dfni in vr vz vzeta coordinates
+        dfni = dfni_sym(Lr,Lz,Lvperp,Lvpa,r_bc,z_bc,geometry,nr)
+        vrvzvzeta_dfni = cartesian_dfni_sym(Lr,Lz,Lvperp,Lvpa,r_bc,z_bc,geometry,nr) #dfni in vr vz vzeta coordinates
         
         # neutral manufactured solutions
         densn = densn_sym(Lr,Lz,r_bc,z_bc,geometry,composition)
-        dfnn = dfnn_sym(Lr,Lz,r_bc,z_bc,geometry,composition)
-        gav_dfnn = gyroaveraged_dfnn_sym(Lr,Lz,r_bc,z_bc,geometry,composition) # gyroaverage < dfnn > in vpa vperp coordinates
+        dfnn = dfnn_sym(Lr,Lz,Lvzeta,Lvr,Lvz,r_bc,z_bc,geometry,composition)
+        gav_dfnn = gyroaveraged_dfnn_sym(Lr,Lz,Lvzeta,Lvr,Lvz,r_bc,z_bc,geometry,composition) # gyroaverage < dfnn > in vpa vperp coordinates
         
         # define derivative operators
         Dr = Differential(r) 
         Dz = Differential(z) 
         Dvpa = Differential(vpa) 
         Dvperp = Differential(vperp) 
-        Dt = Differential(t) 
     
         # get geometric/composition data
         Bzed = geometry.Bzed
@@ -209,13 +296,81 @@ using Symbolics
         Er = -Dr(phi)
         Ez = -Dz(phi)
     
-        # the ion source to maintain the manufactured solution
-        Si = ( Dt(dfni) + ( vpa * (Bzed/Bmag) - 0.5*rhostar*Er*rfac ) * Dz(dfni) + ( 0.5*rhostar*Ez*rfac ) * Dr(dfni) + ( 0.5*Ez*Bzed/Bmag ) * Dvpa(dfni)
-               + cx_frequency*( densn*dfni - densi*gav_dfnn ) ) - ionization_frequency*dense*gav_dfnn 
+        rhs_ion = 0 * z
+        if advance === nothing || advance.vpa_advection
+            rhs_ion += - ( half*Ez*Bzed/Bmag ) * Dvpa(dfni)
+        end
+        if advance === nothing || advance.z_advection
+            rhs_ion += -( vpa * (Bzed/Bmag) - half*rhostar*Er ) * Dz(dfni)
+        end
+        if advance === nothing || advance.r_advection
+            rhs_ion += - ( half*rhostar*Ez ) * Dr(dfni)
+        end
+        if advance === nothing || advance.cx_collisions
+            rhs_ion += cx_frequency * (gav_dfnn*densi - dfni*densn)
+        end
+        #if advance === nothing || advance.ionization_collsions
+        #    # placeholder
+        #end
+        #if advance === nothing || advance.source_terms
+        #    # placeholder
+        #end
+        #if advance === nothing || advance.continuity
+        #    # placeholder
+        #end
+        #if advance === nothing || advance.force_balance
+        #    # placeholder
+        #end
+        #if advance === nothing || advance.energy
+        #    # placeholder
+        #end
+        rhs_neutral = 0 * z
+        if advance == nothing || advance.neutral_z_advection
+            rhs_neutral += -vz * Dz(dfnn)
+        end
+        if advance == nothing || advance.neutral_r_advection
+            rhs_neutral += - rfac*vr * Dr(dfnn)
+        end
+        if advance == nothing || advance.cx_collisions
+            rhs_neutral += - cx_frequency* (densi*dfnn - densn*vrvzvzeta_dfni)
+        end
+        if advance == nothing || advance.ionization_collisions
+            rhs_neutral += - ionization_frequency*dense*dfnn
+        end
+
+        return expand_derivatives(rhs_ion), expand_derivatives(rhs_neutral)
+    end
+
+    function manufactured_rhs(Lr::mk_float, Lz::mk_float, Lvperp::mk_float,
+                              Lvpa::mk_float, Lvzeta::mk_float, Lvr::mk_float,
+                              Lvz::mk_float, r_bc::String, z_bc::String,
+                              composition::species_composition,
+                              geometry::geometry_input, collisions::collisions_input,
+                              nr::mk_int, advance::Union{advance_info,Nothing}=nothing)
+        rhs_ion_sym, rhs_neutral_sym = manufactured_rhs_sym(Lr, Lz, Lvperp, Lvpa,
+            Lvzeta, Lvr, Lvz, r_bc, z_bc, composition, geometry, collisions, nr,
+            advance)
+        return build_function(rhs_ion_sym, vpa, vperp, z, r, t, expression=Val{false}),
+               build_function(rhs_neutral_sym, vz, vr, vzeta, z, r, t, expression=Val{false})
+    end
+
+    function manufactured_sources(Lr::mk_float, Lz::mk_float, Lvperp::mk_float,
+                                  Lvpa::mk_float, Lvzeta::mk_float, Lvr::mk_float,
+                                  Lvz::mk_float, r_bc::String, z_bc::String,
+                                  composition::species_composition,
+                                  geometry::geometry_input,
+                                  collisions::collisions_input, nr::mk_int)
+
+        dfni = dfni_sym(Lr,Lz,Lvperp,Lvpa,r_bc,z_bc,geometry,nr)
+        dfnn = dfnn_sym(Lr,Lz,Lvzeta,Lvr,Lvz,r_bc,z_bc,geometry,composition)
+        rhs_ion, rhs_neutral = manufactured_rhs_sym(Lr, Lz, Lvperp, Lvpa, Lvzeta, Lvr,
+            Lvz, r_bc, z_bc, composition, geometry, collisions, nr)
+
+        Dt = Differential(t)
+
+        Si = Dt(dfni) - rhs_ion
         Source_i = expand_derivatives(Si)
-        
-        # the neutral source to maintain the manufactured solution
-        Sn = Dt(dfnn) + vz * Dz(dfnn) + rfac*vr * Dr(dfnn) + cx_frequency* (densi*dfnn - densn*vrvzvzeta_dfni) + ionization_frequency*dense*dfnn
+        Sn = Dt(dfnn) - rhs_neutral
         Source_n = expand_derivatives(Sn)
         
         Source_i_func = build_function(Source_i, vpa, vperp, z, r, t, expression=Val{false})
@@ -225,5 +380,83 @@ using Symbolics
         
         return manufactured_sources_list
     end 
-    
+
+    """
+        manufactured_solutions_as_arrays(
+            t::mk_float, r::coordinate, z::coordinate, vperp::coordinate,
+            vpa::coordinate)
+
+    Create array filled with manufactured solutions.
+
+    Returns
+    -------
+    (densi, phi, dfni)
+    """
+    function manufactured_solutions_as_arrays(
+        t::mk_float, r::coordinate, z::coordinate, vperp::coordinate,
+        vpa::coordinate, vzeta::coordinate, vr::coordinate, vz::coordinate)
+
+        dfni_func, densi_func = manufactured_solutions(r.L, z.L, vperp.L, vpa.L, vzeta.L,
+                                                       vr.L, vz.L, r.bc, z.bc)
+
+        densi = allocate_shared_float(z.n, r.n)
+        dfni = allocate_shared_float(vpa.n, vperp.n, z.n, r.n)
+
+        begin_r_z_region()
+        @loop_r_z ir iz begin
+            densi[iz,ir] = densi_func(z.grid[iz], r.grid[ir], t)
+            @loop_vperp_vpa ivperp ivpa begin
+                dfni[ivpa,ivperp,iz,ir] = dfni_func(vpa.grid[ivpa], vperp.grid[ivperp], z.grid[iz],
+                                        r.grid[ir], t)
+            end
+        end
+
+        phi = log.(densi)
+
+        return densi, phi, dfni
+    end
+
+    """
+        manufactured_rhs_as_array(
+            t::mk_float, r::coordinate, z::coordinate, vperp::coordinate, vpa::coordinate,
+            composition::species_composition, geometry::geometry_input,
+            collisions::collisions_input, advance::Union{advance_info,Nothing})
+
+    Create arrays filled with manufactured rhs.
+
+    Returns
+    -------
+    rhs_ion, rhs_neutral
+    """
+    function manufactured_rhs_as_array(
+        t::mk_float, r::coordinate, z::coordinate, vperp::coordinate, vpa::coordinate,
+        vzeta::coordinate, vr::coordinate, vz::coordinate,
+        composition::species_composition, geometry::geometry_input,
+        collisions::collisions_input, advance::Union{advance_info,Nothing})
+
+        rhs_ion_func, rhs_neutral_func = manufactured_rhs(r.L, z.L, vperp.L, vpa.L,
+            vzeta.L, vr.L, vz.L, r.bc, z.bc, composition, geometry, collisions, r.n,
+            advance)
+
+        rhs_ion = allocate_shared_float(vpa.n, vperp.n, z.n, r.n)
+
+        begin_r_z_vperp_vpa_region()
+        @loop_r_z_vperp_vpa ir iz ivperp ivpa begin
+            rhs_ion[ivpa,ivperp,iz,ir] =
+            rhs_ion_func(vpa.grid[ivpa], vperp.grid[ivperp], z.grid[iz],
+                         r.grid[ir], t)
+        end
+
+        rhs_neutral = allocate_shared_float(vz.n, vr.n, vzeta.n, z.n, r.n)
+
+        begin_r_z_vzeta_vr_vz_region()
+        @loop_r_z_vzeta_vr_vz ir iz ivzeta ivr ivz begin
+            rhs_neutral[ivz,ivr,ivzeta,iz,ir] =
+            rhs_neutral_func(vz.grid[ivz], vr.grid[ivr], vzeta.grid[ivzeta],
+                             z.grid[iz], r.grid[ir], t)
+        end
+
+        return rhs_ion, rhs_neutral
+    end
+
 end
diff --git a/src/moment_kinetics.jl b/src/moment_kinetics.jl
index f48c336e3..d844769e8 100644
--- a/src/moment_kinetics.jl
+++ b/src/moment_kinetics.jl
@@ -22,9 +22,9 @@ include("chebyshev.jl")
 include("finite_differences.jl")
 include("quadrature.jl")
 include("calculus.jl")
-include("file_io.jl")
 include("input_structs.jl")
 include("coordinates.jl")
+include("file_io.jl")
 include("velocity_moments.jl")
 include("velocity_grid_transforms.jl")
 include("em_fields.jl")
@@ -129,7 +129,7 @@ Perform all the initialization steps for a run.
 function setup_moment_kinetics(input_dict::Dict)
     # Set up MPI
     initialize_comms!()
-    
+
     input = mk_input(input_dict)
     # obtain input options from moment_kinetics_input.jl
     # and check input to catch errors
@@ -139,21 +139,29 @@ function setup_moment_kinetics(input_dict::Dict)
         vz_input, vr_input, vzeta_input, 
         composition, species, collisions, geometry, drive_input = input
     # initialize z grid and write grid point locations to file
-    z = define_coordinate(z_input, composition)
+    z, z_spectral = define_coordinate(z_input, composition)
     # initialize r grid and write grid point locations to file
-    r = define_coordinate(r_input, composition)
+    r, r_spectral = define_coordinate(r_input, composition)
     # initialize vpa grid and write grid point locations to file
-    vpa = define_coordinate(vpa_input, composition)
+    vpa, vpa_spectral = define_coordinate(vpa_input, composition)
     # initialize vperp grid and write grid point locations to file
-    vperp = define_coordinate(vperp_input, composition)
+    vperp, vperp_spectral = define_coordinate(vperp_input, composition)
     # initialize gyrophase grid and write grid point locations to file
-    gyrophase = define_coordinate(gyrophase_input, composition)
+    gyrophase, gyrophase_spectral = define_coordinate(gyrophase_input, composition)
     # initialize vz grid and write grid point locations to file
-    vz = define_coordinate(vz_input, composition)
+    vz, vz_spectral = define_coordinate(vz_input, composition)
     # initialize vr grid and write grid point locations to file
-    vr = define_coordinate(vr_input, composition)
+    vr, vr_spectral = define_coordinate(vr_input, composition)
     # initialize vr grid and write grid point locations to file
-    vzeta = define_coordinate(vzeta_input, composition)
+    vzeta, vzeta_spectral = define_coordinate(vzeta_input, composition)
+
+    ##
+    # construct named list of spectral objects to compactify arguments
+    ##
+    spectral_objects = (vz_spectral = vz_spectral, vr_spectral = vr_spectral,
+                        vzeta_spectral = vzeta_spectral, vpa_spectral = vpa_spectral,
+                        vperp_spectral = vperp_spectral, z_spectral = z_spectral,
+                        r_spectral = r_spectral)
     # Create loop range variables for shared-memory-parallel loops
     if composition.n_neutral_species == 0
         n_neutral_loop_size = 1 # Need this to have looping setup. Avoid neutral loops with if statements.
@@ -169,9 +177,10 @@ function setup_moment_kinetics(input_dict::Dict)
     code_time = 0.
     # create arrays and do other work needed to setup
     # the main time advance loop -- including normalisation of f by density if requested
-    moments, fields, spectral_objects, advect_objects, 
-    scratch, advance, scratch_dummy, manufactured_source_list = setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, composition,
-        drive_input, moments, t_input, collisions, species, geometry, boundary_distributions)
+    moments, fields, advect_objects, scratch,
+        advance, scratch_dummy, manufactured_source_list = setup_time_advance!(pdf, vz,
+            vr, vzeta, vpa, vperp, z, r, spectral_objects, composition, drive_input,
+            moments, t_input, collisions, species, geometry, boundary_distributions)
     # setup i/o
     io, cdf = setup_file_io(output_dir, run_name, vz, vr, vzeta, vpa, vperp, z, r, composition, collisions)
     # write initial data to ascii files
diff --git a/src/moment_kinetics_input.jl b/src/moment_kinetics_input.jl
index 7ffad5798..ebd099ac8 100644
--- a/src/moment_kinetics_input.jl
+++ b/src/moment_kinetics_input.jl
@@ -362,8 +362,9 @@ function mk_input(scan_input=Dict())
     end
 
     # check input to catch errors/unsupported options
-    check_input(io, output_dir, nstep, dt, use_semi_lagrange,
-        z_immutable, vpa_immutable, composition, species_immutable, evolve_moments)
+    check_input(io, output_dir, nstep, dt, use_semi_lagrange, r_immutable, z_immutable,
+        vpa_immutable, vperp_immutable, gyrophase_immutable, vz_immutable, vr_immutable,
+        vzeta_immutable, composition, species_immutable, evolve_moments)
 
     # return immutable structs for z, vpa, species and composition
     all_inputs = (run_name, output_dir, evolve_moments, t, 
@@ -777,17 +778,22 @@ end
 """
 check various input options to ensure they are all valid/consistent
 """
-function check_input(io, output_dir, nstep, dt, use_semi_lagrange, z, vpa,
-    composition, species, evolve_moments)
+function check_input(io, output_dir, nstep, dt, use_semi_lagrange, r, z, vpa, vperp, gyrophase,
+    vz, vr, vzeta, composition, species, evolve_moments)
     # copy the input file to the output directory to be saved
     if block_rank[] == 0
         cp(joinpath(@__DIR__, "moment_kinetics_input.jl"), joinpath(output_dir, "moment_kinetics_input.jl"), force=true)
     end
     # open ascii file in which informtaion about input choices will be written
     check_input_time_advance(nstep, dt, use_semi_lagrange, io)
-    check_input_z(z, io)
-    check_input_vpa(vpa, io)
-    #check_input_initialization(composition, species, io) MRH Need to update
+    check_coordinate_input(r, "r", io)
+    check_coordinate_input(z, "z", io)
+    check_coordinate_input(vpa, "vpa", io)
+    check_coordinate_input(vperp, "vperp", io)
+    check_coordinate_input(gyrophase, "gyrophase", io)
+    check_coordinate_input(vz, "vz", io)
+    check_coordinate_input(vr, "vr", io)
+    check_coordinate_input(vzeta, "vzeta", io)
     # if the parallel flow is evolved separately, then the density must also be evolved separately
     if evolve_moments.parallel_flow && !evolve_moments.density
         print(io,">evolve_moments.parallel_flow = true, but evolve_moments.density = false.")
@@ -811,69 +817,46 @@ function check_input_time_advance(nstep, dt, use_semi_lagrange, io)
 end
 
 """
+Check input for a coordinate
 """
-function check_input_z(z, io)
-    println(io)
-    println(io,"######## z-grid ########")
-    println(io)
-    # discretization_option determines discretization in z
-    # supported options are chebyshev_pseudospectral and finite_difference
-    if z.discretization == "chebyshev_pseudospectral"
-        print(io,">z.discretization = 'chebyshev_pseudospectral'.  ")
-        println(io,"using a Chebyshev pseudospectral method in z.")
-    elseif z.discretization == "finite_difference"
-        println(io,">z.discretization = 'finite_difference', ",
-            "and z.fd_option = ", z.fd_option,
-            "  using finite differences on an equally spaced grid in z.")
-        fd_check_option(z.fd_option, z.ngrid)
-    else
-        input_option_error("z.discretization", z.discretization)
+function check_coordinate_input(coord, coord_name, io)
+    if coord.ngrid * coord.nelement == 1
+        # Coordinate is not being used for this run
+        return
     end
-    # boundary_option determines z boundary condition
-    # supported options are "constant" and "periodic"
-    if z.bc == "constant"
-        println(io,">z.bc = 'constant'.  enforcing constant incoming BC in z.")
-    elseif z.bc == "periodic"
-        println(io,">z.bc = 'periodic'.  enforcing periodicity in z.")
-    elseif z.bc == "wall"
-        println(io,">z.bc = 'wall'.  enforcing wall BC in z.")
-    else
-        input_option_error("z.bc", z.bc)
-    end
-    println(io,">using ", z.ngrid, " grid points per z element on ", z.nelement,
-        " elements across the z domain [", -0.5*z.L, ",", 0.5*z.L, "].")
-end
 
-"""
-"""
-function check_input_vpa(vpa, io)
     println(io)
-    println(io,"######## vpa-grid ########")
+    println(io,"######## $coord_name-grid ########")
     println(io)
-    # discretization_option determines discretization in vpa
+    # discretization_option determines discretization in coord
     # supported options are chebyshev_pseudospectral and finite_difference
-    if vpa.discretization == "chebyshev_pseudospectral"
-        print(io,">vpa.discretization = 'chebyshev_pseudospectral'.  ")
-        println(io,"using a Chebyshev pseudospectral method in vpa.")
-    elseif vpa.discretization == "finite_difference"
-        println(io,">vpa.discretization = 'finite_difference', and ",
-            "vpa.fd_option = ", vpa.fd_option,
-            "  using finite differences on an equally spaced grid in vpa.")
-        fd_check_option(vpa.fd_option, vpa.ngrid)
+    if coord.discretization == "chebyshev_pseudospectral"
+        print(io,">$coord_name.discretization = 'chebyshev_pseudospectral'.  ")
+        println(io,"using a Chebyshev pseudospectral method in $coord_name.")
+    elseif coord.discretization == "finite_difference"
+        println(io,">$coord_name.discretization = 'finite_difference', ",
+            "and $coord_name.fd_option = ", coord.fd_option,
+            "  using finite differences on an equally spaced grid in $coord_name.")
+        fd_check_option(coord.fd_option, coord.ngrid)
     else
-        input_option_error("vpa.discretization", vpa.discretization)
+        input_option_error("$coord_name.discretization", coord.discretization)
     end
-    # boundary_option determines vpa boundary condition
-    # supported options are "zero" and "periodic"
-    if vpa.bc == "zero"
-        println(io,">vpa.bc = 'zero'.  enforcing zero incoming BC in vpa.")
-    elseif vpa.bc == "periodic"
-        println(io,">vpa.bc = 'periodic'.  enforcing periodicity in vpa.")
+    # boundary_option determines coord boundary condition
+    # supported options are "constant" and "periodic" for spatial dimensions
+    # and "zero" and "periodic" for velocity dimensions
+    if !occursin("v", coord_name) && coord.bc == "constant"
+        println(io,">$coord_name.bc = 'constant'.  enforcing constant incoming BC in $coord_name.")
+    elseif occursin("v", coord_name) && coord.bc == "zero"
+        println(io,">$coord_name.bc = 'zero'.  enforcing constant incoming BC in $coord_name.")
+    elseif coord.bc == "periodic"
+        println(io,">$coord_name.bc = 'periodic'.  enforcing periodicity in $coord_name.")
+    elseif coord_name == "z" && coord.bc == "wall"
+        println(io,">$coord_name.bc = 'wall'.  enforcing wall BC in $coord_name.")
     else
-        input_option_error("vpa.bc", vpa.bc)
+        input_option_error("$coord_name.bc", coord.bc)
     end
-    println(io,">using ", vpa.ngrid, " grid points per vpa element on ", vpa.nelement,
-        " elements across the vpa domain [", -0.5*vpa.L, ",", 0.5*vpa.L, "].")
+    println(io,">using ", coord.ngrid, " grid points per $coord_name element on ", coord.nelement,
+        " elements across the $coord_name domain [", -0.5*coord.L, ",", 0.5*coord.L, "].")
 end
 
 """
diff --git a/src/neutral_advection.jl b/src/neutral_advection.jl
index 97d02b5a4..f61a6f700 100644
--- a/src/neutral_advection.jl
+++ b/src/neutral_advection.jl
@@ -28,7 +28,7 @@ function neutral_advection_r!(f_out, fvec_in, advect, r, z, vzeta, vr, vz, dt, r
         @loop_z_vzeta_vr_vz iz ivzeta ivr ivz begin
             # take the normalized pdf contained in fvec_in.pdf and remove the normalization,
             # returning the true (un-normalized) particle distribution function in r.scratch
-            @. r.scratch = fvec_in.pdf_neutral[ivz,ivr,ivzeta,iz,:,isn]
+            @. r.scratch = @views fvec_in.pdf_neutral[ivz,ivr,ivzeta,iz,:,isn]
 
             @views advance_f_local!(f_out[ivz,ivr,ivzeta,iz,:,isn], r.scratch,
                                     advect[isn], ivz, ivr, ivzeta, iz,
@@ -86,7 +86,7 @@ function neutral_advection_z!(f_out, fvec_in, advect, r, z, vzeta, vr, vz, dt, z
         @loop_r_vzeta_vr_vz ir ivzeta ivr ivz begin
             # take the normalized pdf contained in fvec_in.pdf and remove the normalization,
             # returning the true (un-normalized) particle distribution function in r.scratch
-            @. z.scratch = fvec_in.pdf_neutral[ivz,ivr,ivzeta,:,ir,isn]
+            @. z.scratch = @views fvec_in.pdf_neutral[ivz,ivr,ivzeta,:,ir,isn]
 
             @views advance_f_local!(f_out[ivz,ivr,ivzeta,:,ir,isn], z.scratch,
                                     advect[isn], ivz, ivr, ivzeta, ir,
diff --git a/src/post_processing.jl b/src/post_processing.jl
index f6ac5bc24..d239e1fd3 100644
--- a/src/post_processing.jl
+++ b/src/post_processing.jl
@@ -55,6 +55,28 @@ function moving_average(v::AbstractVector, n::mk_int)
     return result
 end
 
+"""
+    L2_error_norm(a, b)
+
+Calculate the L2 norm of the error between a and b: sqrt(mean((a-b)^2))
+"""
+function L2_error_norm(a, b)
+    @assert size(a) == size(b)
+    error = @. (a-b)
+    return sqrt(mean(error.^2))
+end
+
+"""
+    L_infinity_error_norm(a, b)
+
+Calculate the L_infinity norm of the error between a and b: maximum(|a-b|)
+"""
+function L_infinity_error_norm(a, b)
+    @assert size(a) == size(b)
+    error = @. (a-b)
+    return maximum(abs.(error))
+end
+
 """
 """
 function analyze_and_plot_data(path)
diff --git a/src/r_advection.jl b/src/r_advection.jl
index 8e2d8a1d4..2574399d0 100644
--- a/src/r_advection.jl
+++ b/src/r_advection.jl
@@ -27,7 +27,7 @@ function r_advection!(f_out, fvec_in, fields, advect, r, z, vperp, vpa,
         @loop_z_vperp_vpa iz ivperp ivpa begin
             # take the normalized pdf contained in fvec_in.pdf and remove the normalization,
             # returning the true (un-normalized) particle distribution function in r.scratch
-            @. r.scratch = fvec_in.pdf[ivpa,ivperp,iz,:,is]
+            @. r.scratch = @views fvec_in.pdf[ivpa,ivperp,iz,:,is]
 
             @views advance_f_local!(f_out[ivpa,ivperp,iz,:,is], r.scratch,
                                     advect[is], ivpa, ivperp, iz,
@@ -49,7 +49,7 @@ function update_speed_r!(advect, fields, vpa, vperp, z, r, geometry)
         ExBfac = 0.5*geometry.rstar
         @inbounds begin
             @loop_z_vperp_vpa iz ivperp ivpa begin
-                @views advect.speed[:,ivpa,ivperp,iz] .= ExBfac*fields.Ez[iz,:]
+                @views @. advect.speed[:,ivpa,ivperp,iz] = ExBfac*fields.Ez[iz,:]
             end
         end
     elseif r.advection.option == "default" && r.n == 1 
diff --git a/src/time_advance.jl b/src/time_advance.jl
index 9935546e8..e10cb28c3 100644
--- a/src/time_advance.jl
+++ b/src/time_advance.jl
@@ -12,7 +12,6 @@ using ..debugging
 using ..file_io: write_data_to_ascii, write_data_to_binary, debug_dump
 using ..looping
 using ..moment_kinetics_structs: scratch_pdf
-using ..chebyshev: setup_chebyshev_pseudospectral
 using ..chebyshev: chebyshev_derivative!
 using ..velocity_moments: update_moments!, reset_moments_status!
 using ..velocity_moments: enforce_moment_constraints!
@@ -25,6 +24,7 @@ using ..initial_conditions: enforce_z_boundary_condition!, enforce_boundary_cond
 using ..initial_conditions: enforce_vpa_boundary_condition!, enforce_r_boundary_condition!
 using ..initial_conditions: enforce_neutral_boundary_conditions!
 using ..initial_conditions: enforce_neutral_z_boundary_condition!, enforce_neutral_r_boundary_condition!
+using ..input_structs: advance_info, time_input
 using ..advection: setup_advection, update_boundary_indices!
 using ..z_advection: update_speed_z!, z_advection!
 using ..r_advection: update_speed_r!, r_advection!
@@ -42,29 +42,8 @@ using ..em_fields: setup_em_fields, update_phi!
 
 using ..manufactured_solns: manufactured_sources
 using ..advection: advection_info
-
 @debug_detect_redundant_block_synchronize using ..communication: debug_detect_redundant_is_active
 
-"""
-"""
-mutable struct advance_info
-    vpa_advection::Bool
-    z_advection::Bool
-    r_advection::Bool
-    neutral_z_advection::Bool
-    neutral_r_advection::Bool
-    cx_collisions::Bool
-    cx_collisions_1V::Bool
-    ionization_collisions::Bool
-    ionization_collisions_1V::Bool
-    source_terms::Bool
-    continuity::Bool
-    force_balance::Bool
-    energy::Bool
-    rk_coefs::Array{mk_float,2}
-    manufactured_solns_test::Bool
-end
-
 mutable struct scratch_dummy_arrays
     dummy_sr::Array{mk_float,2}
     dummy_vpavperp::Array{mk_float,2}
@@ -99,8 +78,9 @@ this includes creating and populating structs
 for Chebyshev transforms, velocity space moments,
 EM fields, semi-Lagrange treatment, and advection terms
 """
-function setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, composition, drive_input, moments,
-                             t_input, collisions, species, geometry, boundary_distributions)
+function setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, spectral_objects,
+                             composition, drive_input, moments, t_input, collisions,
+                             species, geometry, boundary_distributions)
     # define some local variables for convenience/tidiness
     n_species = composition.n_species
     n_ion_species = composition.n_ion_species
@@ -158,91 +138,6 @@ function setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, composition,
                            manufactured_solns_test)
 
     
-    if z.discretization == "chebyshev_pseudospectral"
-        # create arrays needed for explicit Chebyshev pseudospectral treatment in vpa
-        # and create the plans for the forward and backward fast Chebyshev transforms
-        z_spectral = setup_chebyshev_pseudospectral(z)
-        # obtain the local derivatives of the uniform z-grid with respect to the used z-grid
-        chebyshev_derivative!(z.duniform_dgrid, z.uniform_grid, z_spectral, z)
-    else
-        # create dummy Bool variable to return in place of the above struct
-        z_spectral = false
-        z.duniform_dgrid .= 1.0
-    end
-    
-    if r.discretization == "chebyshev_pseudospectral"
-        # create arrays needed for explicit Chebyshev pseudospectral treatment in vpa
-        # and create the plans for the forward and backward fast Chebyshev transforms
-        r_spectral = setup_chebyshev_pseudospectral(r)
-        # obtain the local derivatives of the uniform r-grid with respect to the used r-grid
-        chebyshev_derivative!(r.duniform_dgrid, r.uniform_grid, r_spectral, r)
-    else
-        # create dummy Bool variable to return in place of the above struct
-        r_spectral = false
-        r.duniform_dgrid .= 1.0
-    end
-    
-    if vpa.discretization == "chebyshev_pseudospectral"
-        # create arrays needed for explicit Chebyshev pseudospectral treatment in vpa
-        # and create the plans for the forward and backward fast Chebyshev transforms
-        vpa_spectral = setup_chebyshev_pseudospectral(vpa)
-        # obtain the local derivatives of the uniform vpa-grid with respect to the used vpa-grid
-        chebyshev_derivative!(vpa.duniform_dgrid, vpa.uniform_grid, vpa_spectral, vpa)
-    else
-        # create dummy Bool variable to return in place of the above struct
-        vpa_spectral = false
-        vpa.duniform_dgrid .= 1.0
-    end
-    
-    # MRH CONSIDER REMOVING -> vperp_spectral never used
-    if vperp.discretization == "chebyshev_pseudospectral" && vperp.n > 1
-        # create arrays needed for explicit Chebyshev pseudospectral treatment in vperp
-        # and create the plans for the forward and backward fast Chebyshev transforms
-        vperp_spectral = setup_chebyshev_pseudospectral(vperp)
-        # obtain the local derivatives of the uniform vperp-grid with respect to the used vperp-grid
-        chebyshev_derivative!(vperp.duniform_dgrid, vperp.uniform_grid, vperp_spectral, vperp)
-    else
-        # create dummy Bool variable to return in place of the above struct
-        vperp_spectral = false
-        vperp.duniform_dgrid .= 1.0
-    end
-    
-    if vz.discretization == "chebyshev_pseudospectral" 
-        # create arrays needed for explicit Chebyshev pseudospectral treatment in vz
-        # and create the plans for the forward and backward fast Chebyshev transforms
-        vz_spectral = setup_chebyshev_pseudospectral(vz)
-        # obtain the local derivatives of the uniform vz-grid with respect to the used vz-grid
-        chebyshev_derivative!(vz.duniform_dgrid, vz.uniform_grid, vz_spectral, vz)
-    else
-        # create dummy Bool variable to return in place of the above struct
-        vz_spectral = false
-        vz.duniform_dgrid .= 1.0
-    end
-    
-    if vr.discretization == "chebyshev_pseudospectral" && vr.n > 1
-        # create arrays needed for explicit Chebyshev pseudospectral treatment in vr
-        # and create the plans for the forward and backward fast Chebyshev transforms
-        vr_spectral = setup_chebyshev_pseudospectral(vr)
-        # obtain the local derivatives of the uniform vr-grid with respect to the used vr-grid
-        chebyshev_derivative!(vr.duniform_dgrid, vr.uniform_grid, vr_spectral, vr)
-    else
-        # create dummy Bool variable to return in place of the above struct
-        vr_spectral = false
-        vr.duniform_dgrid .= 1.0
-    end
-    
-    if vzeta.discretization == "chebyshev_pseudospectral" && vzeta.n > 1
-        # create arrays needed for explicit Chebyshev pseudospectral treatment in vzeta
-        # and create the plans for the forward and backward fast Chebyshev transforms
-        vzeta_spectral = setup_chebyshev_pseudospectral(vzeta)
-        # obtain the local derivatives of the uniform vzeta-grid with respect to the used vzeta-grid
-        chebyshev_derivative!(vzeta.duniform_dgrid, vzeta.uniform_grid, vzeta_spectral, vzeta)
-    else
-        # create dummy Bool variable to return in place of the above struct
-        vzeta_spectral = false
-        vzeta.duniform_dgrid .= 1.0
-    end
-    
     # create an array of structs containing scratch arrays for the pdf and low-order moments
     # that may be evolved separately via fluid equations
     scratch = setup_scratch_arrays(moments, pdf.charged.norm, pdf.neutral.norm, t_input.n_rk_stages)
@@ -256,7 +151,8 @@ function setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, composition,
     fields = setup_em_fields(z.n, r.n, drive_input.force_phi, drive_input.amplitude, drive_input.frequency)
     # initialize the electrostatic potential
     begin_serial_region()
-    update_phi!(fields, scratch[1], z, r, composition, z_spectral, r_spectral)
+    update_phi!(fields, scratch[1], z, r, composition, spectral_objects.z_spectral,
+                spectral_objects.r_spectral)
     @serial_region begin
         # save the initial phi(z) for possible use later (e.g., if forcing phi)
         fields.phi0 .= fields.phi
@@ -371,17 +267,15 @@ function setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, composition,
     end
     
     ##
-    # construct named list of advect & spectral objects to compactify arguments
+    # construct named list of advect objects to compactify arguments
     ##
     
     #advect_objects = (vpa_advect = vpa_advect, vperp_advect = vperp_advect, z_advect = z_advect, 
     # r_advect = r_advect, neutral_z_advect = neutral_z_advect, neutral_r_advect = neutral_r_advect)
     advect_objects = advect_object_struct(vpa_advect, vperp_advect, z_advect, r_advect, neutral_z_advect, neutral_r_advect)
-    #spectral_objects = (vz_spectral = vz_spectral, vr_spectral = vr_spectral, vzeta_spectral = vzeta_spectral,
-    # vpa_spectral = vpa_spectral, vperp_spectral = vperp_spectral, z_spectral = z_spectral, r_spectral = r_spectral)
-    spectral_objects = spectral_object_struct(vz_spectral, vr_spectral, vzeta_spectral, vpa_spectral, vperp_spectral, z_spectral, r_spectral)
     if(t_input.use_manufactured_solns)
-        manufactured_source_list = manufactured_sources(r.L,z.L,r.bc,z.bc,composition,geometry,collisions,r.n)
+        manufactured_source_list = manufactured_sources(r.L, z.L, vperp.L, vpa.L,
+            vzeta.L, vr.L, vz.L, r.bc, z.bc, composition, geometry, collisions, r.n)
     else
         manufactured_source_list = false # dummy Bool to be passed as argument instead of list
     end
@@ -389,8 +283,8 @@ function setup_time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, composition,
     # Ensure all processes are synchronized at the end of the setup
     _block_synchronize()
 
-    return moments, fields, spectral_objects, advect_objects, 
-    scratch, advance, scratch_dummy, manufactured_source_list
+    return moments, fields, advect_objects, scratch, advance, scratch_dummy,
+           manufactured_source_list
 end
 
 """
@@ -636,7 +530,7 @@ function ssp_rk!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyrophase,
            composition, collisions, geometry, boundary_distributions, 
            advance, scratch_dummy, manufactured_source_list,  istep)#pdf_in,
     
-    begin_s_r_z_vperp_region()
+    begin_s_r_z_region()
     
     n_rk_stages = t_input.n_rk_stages
 
@@ -871,4 +765,35 @@ function update_pdf_unnorm!(pdf, moments, scratch, composition, vpa, vperp)
     end
 end
 
+"""
+Hacky way of extracting df/dt (and dn/dt, etc. if calculated) by wrapping
+euler_time_advance!()
+"""
+function evaluate_ddt!(fvec_out, fvec_in, pdf, fields, moments, advect_objects, vz, vr,
+        vzeta, vpa, vperp, gyrophase, z, r, t, t_input, spectral_objects, composition,
+        collisions, geometry, boundary_distributions, scratch_dummy,
+        manufactured_source_list, advance)
+
+    # Set initial state of fvec_out.pdf and fvec_out.pdf_neutral to 0.0, so that they
+    # will only contain the change in f
+    begin_s_r_z_vperp_vpa_region()
+    @loop_s_r_z_vperp_vpa is ir iz ivperp ivpa begin
+        fvec_out.pdf[ivpa,ivperp,iz,ir,is] = 0.0
+    end
+    begin_sn_r_z_vzeta_vr_vz_region()
+    @loop_sn_r_z_vzeta_vr_vz isn ir iz ivzeta ivr ivz begin
+        fvec_out.pdf_neutral[ivz,ivr,ivzeta,iz,ir,isn] = 0.0
+    end
+
+    # modify t_input to set dt=1
+    modified_t_input = time_input(t_input.nstep, 1.0, t_input.nwrite,
+        t_input.use_semi_lagrange, t_input.n_rk_stages, t_input.split_operators,
+        t_input.use_manufactured_solns)
+
+    euler_time_advance!(fvec_out, fvec_in, pdf, fields, moments, advect_objects, vz, vr,
+        vzeta, vpa, vperp, gyrophase, z, r, t, modified_t_input, spectral_objects,
+        composition, collisions, geometry, boundary_distributions, scratch_dummy,
+        manufactured_source_list, advance, 1)
+end
+
 end
diff --git a/src/z_advection.jl b/src/z_advection.jl
index 2d088e4b6..d7ca2e238 100644
--- a/src/z_advection.jl
+++ b/src/z_advection.jl
@@ -24,7 +24,7 @@ function z_advection!(f_out, fvec_in, fields, advect, z, vpa, vperp, r, dt, t, z
 
         # advance z-advection equation
         @loop_r_vperp_vpa ir ivperp ivpa begin
-            @. z.scratch = fvec_in.pdf[ivpa,ivperp,:,ir,is]
+            @. z.scratch = @views fvec_in.pdf[ivpa,ivperp,:,ir,is]
             @views advance_f_local!(f_out[ivpa,ivperp,:,ir,is], z.scratch, advect[is], ivpa, ivperp, ir,
                                     z, dt, z_spectral)
         end
@@ -47,7 +47,7 @@ function update_speed_z!(advect, fields, vpa, vperp, z, r, t, geometry)
         @inbounds begin
             
             @loop_r_vperp_vpa ir ivperp ivpa begin
-                    @views advect.speed[:,ivpa,ivperp,ir] .= vpa.grid[ivpa]*kpar .+ ExBfac*fields.Er[:,ir]
+                @views @. advect.speed[:,ivpa,ivperp,ir] = vpa.grid[ivpa]*kpar + ExBfac*fields.Er[:,ir]
             end
         
         end
diff --git a/test/calculus_tests.jl b/test/calculus_tests.jl
index faf35b836..8bd1545e6 100644
--- a/test/calculus_tests.jl
+++ b/test/calculus_tests.jl
@@ -4,20 +4,16 @@ include("setup.jl")
 
 using moment_kinetics.input_structs: grid_input, advection_input
 using moment_kinetics.coordinates: define_coordinate
-using moment_kinetics.chebyshev: setup_chebyshev_pseudospectral
 using moment_kinetics.calculus: derivative!, integral
 
 using Random
 
-fd_fake_setup(x) = return false
-
 function runtests()
     @testset "calculus" verbose=use_verbose begin
         println("calculus tests")
         @testset "fundamental theorem of calculus" begin
-            @testset "$discretization $ngrid $nelement" for (discretization, setup_func) ∈
-                    (("finite_difference", fd_fake_setup),
-                     ("chebyshev_pseudospectral", setup_chebyshev_pseudospectral)),
+            @testset "$discretization $ngrid $nelement" for
+                    discretization ∈ ("finite_difference", "chebyshev_pseudospectral"),
                     ngrid ∈ (5,6,7,8,9,10), nelement ∈ (1, 2, 3, 4, 5)
 
                 if discretization == "finite_difference" && (ngrid - 1) * nelement % 2 == 1
@@ -41,11 +37,7 @@ function runtests()
                 input = grid_input("coord", ngrid, nelement, L,
                     discretization, fd_option, bc, adv_input)
                 # create the coordinate struct 'x'
-                x = define_coordinate(input)
-                # create arrays needed for Chebyshev pseudospectral treatment in x
-                # and create the plans for the forward and backward fast Chebyshev
-                # transforms
-                spectral = setup_func(x)
+                x, spectral = define_coordinate(input)
                 # create array for the function f(x) to be differentiated/integrated
                 f = Array{Float64,1}(undef, x.n)
                 # create array for the derivative df/dx
@@ -87,7 +79,7 @@ function runtests()
                 input = grid_input("coord", ngrid, nelement, L,
                     "finite_difference", fd_option, bc, adv_input)
                 # create the coordinate struct 'x'
-                x = define_coordinate(input)
+                x, spectral = define_coordinate(input)
 
                 # create array for the derivative df/dx and the expected result
                 df = Array{Float64,1}(undef, x.n)
@@ -129,7 +121,7 @@ function runtests()
                 input = grid_input("coord", ngrid, nelement, L,
                     "finite_difference", fd_option, bc, adv_input)
                 # create the coordinate struct 'x'
-                x = define_coordinate(input)
+                x, spectral = define_coordinate(input)
 
                 # create array for the derivative df/dx and the expected result
                 df = Array{Float64,1}(undef, x.n)
@@ -167,7 +159,7 @@ function runtests()
                 input = grid_input("coord", ngrid, nelement, L,
                     "finite_difference", fd_option, bc, adv_input)
                 # create the coordinate struct 'x'
-                x = define_coordinate(input)
+                x, spectral = define_coordinate(input)
 
                 # create array for the derivative df/dx and the expected result
                 df = Array{Float64,1}(undef, x.n)
@@ -213,7 +205,7 @@ function runtests()
                 input = grid_input("coord", ngrid, nelement, L,
                     "finite_difference", fd_option, bc, adv_input)
                 # create the coordinate struct 'x'
-                x = define_coordinate(input)
+                x, spectral = define_coordinate(input)
 
                 # create array for the derivative df/dx and the expected result
                 df = Array{Float64,1}(undef, x.n)
@@ -416,11 +408,7 @@ function runtests()
                 input = grid_input("coord", ngrid, nelement, L,
                     "chebyshev_pseudospectral", fd_option, bc, adv_input)
                 # create the coordinate struct 'x'
-                x = define_coordinate(input)
-                # create arrays needed for Chebyshev pseudospectral treatment in x
-                # and create the plans for the forward and backward fast Chebyshev
-                # transforms
-                spectral = setup_chebyshev_pseudospectral(x)
+                x, spectral = define_coordinate(input)
 
                 offset = randn(rng)
                 f = @. sinpi(2.0 * x.grid / L) + offset
@@ -609,11 +597,7 @@ function runtests()
                 input = grid_input("coord", ngrid, nelement, L,
                     "chebyshev_pseudospectral", fd_option, bc, adv_input)
                 # create the coordinate struct 'x'
-                x = define_coordinate(input)
-                # create arrays needed for Chebyshev pseudospectral treatment in x
-                # and create the plans for the forward and backward fast Chebyshev
-                # transforms
-                spectral = setup_chebyshev_pseudospectral(x)
+                x, spectral = define_coordinate(input)
 
                 offset = randn(rng)
                 f = @. sinpi(2.0 * x.grid / L) + offset
@@ -650,11 +634,7 @@ function runtests()
                 input = grid_input("coord", ngrid, nelement, L,
                     "chebyshev_pseudospectral", fd_option, bc, adv_input)
                 # create the coordinate struct 'x'
-                x = define_coordinate(input)
-                # create arrays needed for Chebyshev pseudospectral treatment in x
-                # and create the plans for the forward and backward fast Chebyshev
-                # transforms
-                spectral = setup_chebyshev_pseudospectral(x)
+                x, spectral = define_coordinate(input)
 
                 # test polynomials up to order ngrid-1
                 for n ∈ 0:ngrid-1
@@ -699,11 +679,7 @@ function runtests()
                 input = grid_input("coord", ngrid, nelement, L,
                     "chebyshev_pseudospectral", fd_option, bc, adv_input)
                 # create the coordinate struct 'x'
-                x = define_coordinate(input)
-                # create arrays needed for Chebyshev pseudospectral treatment in x
-                # and create the plans for the forward and backward fast Chebyshev
-                # transforms
-                spectral = setup_chebyshev_pseudospectral(x)
+                x, spectral = define_coordinate(input)
 
                 # test polynomials up to order ngrid-1
                 for n ∈ 0:ngrid-1
diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl
index ded734b7c..ec70cf647 100644
--- a/test/interpolation_tests.jl
+++ b/test/interpolation_tests.jl
@@ -4,12 +4,9 @@ include("setup.jl")
 
 using moment_kinetics.input_structs: grid_input, advection_input
 using moment_kinetics.coordinates: define_coordinate
-using moment_kinetics.chebyshev: setup_chebyshev_pseudospectral
 using moment_kinetics.interpolation:
     interpolate_to_grid_1d, interpolate_to_grid_z, interpolate_to_grid_vpa
 
-fd_fake_setup(z) = return false
-
 # periodic test function
 # returns an array whose shape is the outer product of the 2nd, 3rd, ... arguments
 test_function(L, coords...) =
@@ -28,20 +25,15 @@ adv_input = advection_input("default", 1.0, 0.0, 0.0)
 function runtests()
     @testset "interpolation" verbose=use_verbose begin
         @testset "$discretization, $ntest, $nelement, $zlim" for
-                (discretization, setup_func, rtol) ∈
-                    (("finite_difference", fd_fake_setup, 1.e-5),
-                     ("chebyshev_pseudospectral", setup_chebyshev_pseudospectral, 1.e-8)),
+                (discretization, rtol) ∈
+                    (("finite_difference", 1.e-5), ("chebyshev_pseudospectral", 1.e-8)),
                     ntest ∈ (3, 14), nelement ∈ (2, 8), zlim ∈ (L/2.0, L/5.0)
 
             # create the 'input' struct containing input info needed to create a coordinate
             input = grid_input("coord", ngrid, nelement, L,
                 discretization, fd_option, bc, adv_input)
             # create the coordinate struct 'z'
-            z = define_coordinate(input)
-            # For Chebyshev method, create arrays needed for Chebyshev pseudospectral
-            # treatment in z and create the plans for the forward and backward fast
-            # Chebyshev transforms. Just get `false` for finite difference.
-            spectral = setup_func(z)
+            z, spectral = define_coordinate(input)
 
             test_grid = [z for z in range(-zlim, zlim, length=ntest)]
 
diff --git a/test/manufactured_ddt_tests.jl b/test/manufactured_ddt_tests.jl
new file mode 100644
index 000000000..01bc35688
--- /dev/null
+++ b/test/manufactured_ddt_tests.jl
@@ -0,0 +1,499 @@
+"""
+Test cases using the method of manufactured solutions (MMS)
+"""
+module ManufacturedDDTTests
+
+include("setup.jl")
+include("mms_utils.jl")
+
+using moment_kinetics: setup_moment_kinetics, cleanup_moment_kinetics!
+using moment_kinetics.looping
+using moment_kinetics.manufactured_solns
+using moment_kinetics.post_processing: L2_error_norm, L_infinity_error_norm
+using moment_kinetics.time_advance: evaluate_ddt!, advance_info
+using moment_kinetics.type_definitions
+using Statistics: mean
+
+# Create a temporary directory for test output
+test_output_directory = tempname()
+mkpath(test_output_directory)
+
+ngrid = 7
+const input_sound_wave_periodic = Dict(
+    "use_manufactured_solns" => true,
+    "n_ion_species" => 1,
+    "n_neutral_species" => 1,
+    "boltzmann_electron_response" => true,
+    "run_name" => "MMS-rperiodic",
+    "base_directory" => test_output_directory,
+    "evolve_moments_density" => false,
+    "evolve_moments_parallel_flow" => false,
+    "evolve_moments_parallel_pressure" => false,
+    "evolve_moments_conservation" => false,
+    "T_e" => 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.62831853071,
+    "ionization_frequency" => 0.0,
+    #"nstep" => 10, #1700,
+    #"dt" => 0.002,
+    #"nwrite" => 10, #1700,
+    "nstep" => 1700, #1700,
+    "dt" => 0.0002, #0.002,
+    "nwrite" => 1700, #1700,
+    "use_semi_lagrange" => false,
+    "n_rk_stages" => 1,
+    "split_operators" => false,
+    "z_ngrid" => ngrid,
+    "z_nelement" => 2,
+    "z_bc" => "periodic",
+    "z_discretization" => "chebyshev_pseudospectral",
+    "r_ngrid" => ngrid,
+    "r_nelement" => 2,
+    "r_bc" => "periodic",
+    "r_discretization" => "chebyshev_pseudospectral",
+    "vpa_ngrid" => ngrid,
+    "vpa_nelement" => 4,
+    "vpa_L" => 8.0,
+    "vpa_bc" => "zero",
+    "vpa_discretization" => "chebyshev_pseudospectral",
+    "vperp_ngrid" => ngrid,
+    "vperp_nelement" => 4,
+    "vperp_L" => 8.0,
+    "vperp_bc" => "zero",
+    "vperp_discretization" => "chebyshev_pseudospectral",
+    "gyrophase_ngrid" => ngrid,
+    "gyrophase_nelement" => 4,
+    "vperp_discretization" => "chebyshev_pseudospectral",
+    "vz_ngrid" => ngrid,
+    "vz_nelement" => 4,
+    "vz_L" => 8.0,
+    "vz_bc" => "zero",
+    "vz_discretization" => "chebyshev_pseudospectral",
+    "vr_ngrid" => ngrid,
+    "vr_nelement" => 4,
+    "vr_L" => 8.0,
+    "vr_bc" => "zero",
+    "vr_discretization" => "chebyshev_pseudospectral",
+    "vzeta_ngrid" => ngrid,
+    "vzeta_nelement" => 4,
+    "vzeta_L" => 8.0,
+    "vzeta_bc" => "zero",
+    "vzeta_discretization" => "chebyshev_pseudospectral",
+)
+
+const input_sound_wave_periodic_no_neutrals =
+    merge(input_sound_wave_periodic, Dict( "n_neutral_species" => 0))
+
+"""
+    evaluate_initial_ddt(input_dict::Dict, advance_input::advance_info)
+
+Evaluate df/dt for the initial state of f.
+
+Very similar to combination of run_moment_kinetics function and various parts of
+time_advance module.
+"""
+function evaluate_initial_ddt(input_dict::Dict, advance_input::advance_info)
+    # set up all the structs, etc. needed for a run
+    mk_state = setup_moment_kinetics(input_dict)
+
+    # Split mk_state tuple into separate variables
+    pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyrophase, z, r, moments,
+    fields, spectral_objects, advect_objects, composition, collisions, geometry,
+    boundary_distributions, advance, scratch_dummy, manufactured_source_list, io,
+    cdf = mk_state
+
+    # Put initial state into scratch[1]
+    begin_s_r_z_region()
+    first_scratch = scratch[1]
+    output = scratch[2]
+    @loop_s_r_z_vperp_vpa is ir iz ivperp ivpa begin
+        first_scratch.pdf[ivpa,ivperp,iz,ir,is] = pdf.charged.norm[ivpa,ivperp,iz,ir,is]
+        output.pdf[ivpa,ivperp,iz,ir,is] = 0.0
+    end
+    @loop_s_r_z is ir iz begin
+        first_scratch.density[iz,ir,is] = moments.charged.dens[iz,ir,is]
+        first_scratch.upar[iz,ir,is] = moments.charged.upar[iz,ir,is]
+        first_scratch.ppar[iz,ir,is] = moments.charged.ppar[iz,ir,is]
+        output.density[iz,ir,is] = 0.0
+        output.upar[iz,ir,is] = 0.0
+        output.ppar[iz,ir,is] = 0.0
+    end
+    begin_sn_r_z_region()
+    @loop_sn_r_z_vzeta_vr_vz isn ir iz ivzeta ivr ivz begin
+        first_scratch.pdf_neutral[ivz,ivr,ivzeta,iz,ir,isn] = pdf.neutral.norm[ivz,ivr,ivzeta,iz,ir,isn]
+        output.pdf_neutral[ivz,ivr,ivzeta,iz,ir,isn] = 0.0
+    end
+    @loop_sn_r_z isn ir iz begin
+        first_scratch.density_neutral[iz,ir,isn] = moments.neutral.dens[iz,ir,isn]
+        output.density_neutral[iz,ir,isn] = 0.0
+    end
+
+    # modify advance_info argument so that no manufactured sources are added when
+    # evaluating the time advance, as in this test we only want manufactured initial
+    # conditions. Also take the terms to evaluate from advance_input so that the calling
+    # function can control which term to use.
+    modified_advance = advance_info(advance_input.vpa_advection,
+        advance_input.z_advection, advance_input.r_advection,
+        advance_input.neutral_z_advection, advance_input.neutral_r_advection,
+        advance_input.cx_collisions, advance_input.cx_collisions_1V,
+        advance_input.ionization_collisions, advance_input.ionization_collisions_1V,
+        advance_input.source_terms, advance_input.continuity,
+        advance_input.force_balance, advance_input.energy, advance.rk_coefs, false)
+
+    evaluate_ddt!(output, first_scratch, pdf, fields, moments, advect_objects, vz, vr,
+        vzeta, vpa, vperp, gyrophase, z, r, t, t_input, spectral_objects, composition,
+        collisions, geometry, boundary_distributions, scratch_dummy,
+        manufactured_source_list, modified_advance)
+
+    # Make a copy of the array with df/dt so it doesn't get messed up by MPI cleanup in
+    # cleanup_moment_kinetics!()
+    dfdt_ion = nothing
+    dfdt_neutral = nothing
+    if global_rank[] == 0
+        dfdt_ion = copy(output.pdf)
+        dfdt_neutral = copy(output.pdf_neutral)
+    end
+
+    rhs_ion_manf, rhs_neutral_manf =
+        manufactured_rhs_as_array(mk_float(0.0), r, z, vperp, vpa,
+           vzeta, vr, vz, composition, geometry, collisions, modified_advance)
+
+    if global_rank[] == 0
+        rhs_ion_manf = copy(rhs_ion_manf)
+        rhs_neutral_manf = copy(rhs_neutral_manf)
+    else
+        rhs_ion_manf = nothing
+        rhs_neutral_manf = nothing
+    end
+
+    # clean up i/o and communications
+    cleanup_moment_kinetics!(io, cdf)
+
+    return dfdt_ion, dfdt_neutral, rhs_ion_manf, rhs_neutral_manf
+end
+
+"""
+    remove_boundary_points(f::AbstractArray{mk_float,4})
+
+Remove points which are modified by boundary conditions* and therefore not
+expected to match the manufactured rhs.
+
+* Boundary conditions are applied at the end of the euler_time_advance!()
+function. In normal operation, these are applied to the distribution function
+(and possibly the moments). Here we hack around to get df/dt stored in the
+fvec_out argument, so the bcs are applied to our df/dt. Physically this makes
+sense for zero-value bcs (as df/dt=0 on the boundary for them), but means that
+the calculated df/dt does not match the one evaluated from the manufactured
+solution, for which evaluating the rhs doesn't necessarily give zero at the
+boundary.
+"""
+function remove_ion_boundary_points(f::AbstractArray{mk_float,4})
+    return f[2:end-1,:,:,:]
+end
+
+"""
+    runcase(input::Dict, advance::advance_info, returnstuff=false)
+
+Run a simulation with parameters set by `input` using manufactured sources and return
+the errors in each variable compared to the manufactured solution.
+"""
+function runcase(input::Dict, advance::advance_info; neutrals, returnstuff=false)
+    dfdt_ion = nothing
+    dfdt_neutral = nothing
+    rhs_ion_manf = nothing
+    rhs_neutral_manf = nothing
+    quietoutput() do
+        dfdt_ion, dfdt_neutral, rhs_ion_manf, rhs_neutral_manf = evaluate_initial_ddt(input, advance)
+    end
+
+    # last 2 elements of mk_state are `io` and `cdf`
+    error_2 = nothing
+    error_inf = nothing
+    if global_rank[] == 0
+        # Only one species, so get rid of species index
+        dfdt_ion = dfdt_ion[:,:,:,:,1]
+        if neutrals
+            dfdt_neutral = dfdt_neutral[:,:,:,:,:,1]
+        end
+
+        # Get rid of vpa boundary points which are (possibly) overwritten by
+        # boundary conditions in the numerically evaluated result.
+        # Not necessary for neutrals as there are no v-space boundary conditions for
+        # neutrals.
+        dfdt_ion = remove_ion_boundary_points(dfdt_ion)
+        rhs_ion_manf = remove_ion_boundary_points(rhs_ion_manf)
+
+        error_2_ion = L2_error_norm(dfdt_ion, rhs_ion_manf)
+        error_inf_ion = L_infinity_error_norm(dfdt_ion, rhs_ion_manf)
+        println("ion error ", error_2_ion, " ", error_inf_ion)
+        if neutrals
+            error_2_neutral = L2_error_norm(dfdt_neutral, rhs_neutral_manf)
+            error_inf_neutral = L_infinity_error_norm(dfdt_neutral, rhs_neutral_manf)
+
+            # Create combined errors
+            error_2 = mean([error_2_ion, error_2_neutral])
+            error_inf = max(error_inf_ion, error_inf_neutral)
+
+            println("neutral error ", error_2_neutral, " ", error_inf_neutral)
+            println("combined error ", error_2, " ", error_inf)
+        else
+            error_2 = error_2_ion
+            error_inf = error_inf_ion
+        end
+    end
+
+    if returnstuff
+        return error_2, error_inf, dfdt_ion, rhs_ion_manf, dfdt_neutral, rhs_neutral_manf
+    else
+        return error_2, error_inf
+    end
+end
+
+"""
+    calculate_convergence_orders(nelements, errors)
+
+Calculate estimated convergence order between consecutive steps of the nelement
+scan.
+
+Don't assume the final error is the 'best' value in case it is affected by
+rounding errors, etc.
+"""
+function calculate_convergence_orders(nelements, errors)
+    @assert size(nelements) == size(errors)
+
+    error_ratios = errors[1:end-1] ./ errors[2:end]
+    nelements_ratios = nelements[2:end] ./ nelements[1:end-1]
+    orders = @. log(error_ratios) / log(nelements_ratios)
+    return orders
+end
+
+"""
+    testconvergence(input::Dict, advance::advance_info; returnstuff::Bool)
+
+Test convergence with spatial resolution
+
+The parameters for the run are given in `input::Dict`. 
+`which_term` controls which term to include in the evolution equation, or include all terms.
+`returnstuff` can be set to true to return the calculated and manfactured RHS of df/dt=(...).
+"""
+function testconvergence(input::Dict, advance::advance_info; ngrid=nothing, returnstuff=false)
+    errors_2 = Vector{mk_float}(undef, 0)
+    errors_inf = Vector{mk_float}(undef, 0)
+
+    if ngrid === nothing
+        ngrid = get_and_check_ngrid(input)
+    else
+        set_ngrid(input, ngrid)
+    end
+
+    # Are there any neutrals?
+    neutrals = input["n_neutral_species"] > 0
+
+    global_rank[] == 0 && println("ngrid=$ngrid")
+
+    n_max = neutrals ? 35 : 50
+    nelement_values = collect(2:2:n_max(ngrid-1))
+    if returnstuff
+        nelement_values = [nelement_values[end]]
+    end
+    lastrhs_ion, lastrhs_manf_ion, lastrhs_neutral, lastrhs_manf_neutral = nothing, nothing, nothing, nothing
+    for nelement ∈ nelement_values
+        global_rank[] == 0 && println("testing nelement=$nelement")
+        case_input = increase_resolution(input, nelement)
+
+        if returnstuff
+            error_2, error_inf, lastrhs_ion, lastrhs_manf_ion, lastrhs_neutral,
+            lastrhs_manf_neutral = runcase(case_input, advance,
+                                           neutrals=neutrals,
+                                           returnstuff=returnstuff)
+        else
+            error_2, error_inf = runcase(case_input, advance,
+                                         neutrals=neutrals)
+        end
+
+        if global_rank[] == 0
+            push!(errors_2, error_2)
+            push!(errors_inf, error_inf)
+        end
+    end
+
+    if global_rank[] == 0
+        convergence_2 = errors_2[1:end] ./ errors_2[end]
+        convergence_inf = errors_inf[1:end] ./ errors_inf[end]
+        expected_convergence = @. (nelement_values[end] / nelement_values[1:end])^(ngrid - 1)
+        println("errors")
+        println(errors_2)
+        println(errors_inf)
+        println("convergence")
+        println(convergence_2)
+        println(convergence_inf)
+        println("expected convergence")
+        println(expected_convergence)
+        println("convergence orders")
+        println(calculate_convergence_orders(nelement_values, errors_2))
+        println(calculate_convergence_orders(nelement_values, errors_inf))
+    end
+
+    if returnstuff
+        return lastrhs_ion, lastrhs_manf_ion, lastrhs_neutral, lastrhs_manf_neutral
+    else
+        return nothing
+    end
+end
+
+function setup_advance(which_term, neutrals=true)
+    advance = advance_info(false, false, false, false, false, false, false, false,
+        false, false, false, false, false, zeros(1,1), false)
+
+    if which_term == :all
+        advance.vpa_advection = true
+        advance.z_advection = true
+        advance.r_advection = true
+        advance.neutral_z_advection = neutrals
+        advance.neutral_r_advection = neutrals
+        advance.cx_collisions = neutrals
+        #advance.ionization_collisions = neutrals
+        advance.source_terms = true
+        advance.continuity = true
+        advance.force_balance = true
+        advance.energy = true
+    elseif which_term == :vpa_advection
+        advance.vpa_advection = true
+    elseif which_term == :z_advection
+        advance.z_advection = true
+    elseif which_term == :r_advection
+        advance.r_advection = true
+    elseif which_term == :neutral_z_advection
+        advance.neutral_z_advection = true
+    elseif which_term == :neutral_r_advection
+        advance.neutral_r_advection = true
+    elseif which_term == :cx_collisions
+        advance.cx_collisions = true
+    elseif which_term == :ionization_collisions
+        advance.ionization_collisions = true
+    elseif which_term == :source_terms
+        advance.source_terms = true
+    elseif which_term == :continuity
+        advance.continuity = true
+    elseif which_term == :force_balance
+        advance.force_balance = true
+    elseif which_term == :energy
+        advance.energy = true
+    end
+
+    return advance
+end
+
+function runtests(ngrid=nothing)
+    @testset "MMS" verbose=use_verbose begin
+        global_rank[] == 0 && println("MMS tests")
+
+        for (input, label, neutrals) ∈ (
+                                        (input_sound_wave_periodic_no_neutrals, "no neutrals", false),
+                                        (input_sound_wave_periodic, "default", true),
+                                       )
+            global_rank[] == 0 && println("$label case")
+            @testset "r-periodic, z-periodic $label" begin
+                @testset "vpa_advection" begin
+                    global_rank[] == 0 && println("\nvpa_advection")
+                    testconvergence(input,
+                                    setup_advance(:vpa_advection, neutrals), ngrid=ngrid)
+                end
+                @testset "z_advection" begin
+                    global_rank[] == 0 && println("\nz_advection")
+                    testconvergence(input, setup_advance(:z_advection, neutrals),
+                                    ngrid=ngrid)
+                end
+                @testset "r_advection" begin
+                    global_rank[] == 0 && println("\nr_advection")
+                    testconvergence(input, setup_advance(:r_advection, neutrals),
+                                    ngrid=ngrid)
+                end
+                if neutrals
+                    @testset "neutral_z_advection" begin
+                        global_rank[] == 0 && println("\nneutral_z_advection")
+                        testconvergence(input,
+                                        setup_advance(:neutral_z_advection, neutrals), ngrid=ngrid)
+                    end
+                    @testset "neutral_r_advection" begin
+                        global_rank[] == 0 && println("\nneutral_r_advection")
+                        testconvergence(input,
+                                        setup_advance(:neutral_r_advection, neutrals), ngrid=ngrid)
+                    end
+                    @testset "cx_collisions" begin
+                        global_rank[] == 0 && println("\ncx_collisions")
+                        testconvergence(input,
+                                        setup_advance(:cx_collisions, neutrals), ngrid=ngrid)
+                    end
+                    @testset "ionization_collisions" begin
+                        global_rank[] == 0 && println("\nionization_collisions")
+                        testconvergence(input,
+                                        setup_advance(:ionization_collisions, neutrals), ngrid=ngrid)
+                    end
+                end
+                #@testset "continuity" begin
+                #    global_rank[] == 0 && println("\ncontinuity")
+                #    testconvergence(input, setup_advance(:continuity, neutrals),
+                #                    ngrid=ngrid)
+                #end
+                #@testset "force_balance" begin
+                #    global_rank[] == 0 && println("\nforce_balance")
+                #    testconvergence(input,
+                #                    setup_advance(:force_balance, neutrals), ngrid=ngrid)
+                #end
+                #@testset "energy" begin
+                #    global_rank[] == 0 && println("\nenergy")
+                #    testconvergence(input, setup_advance(:energy, neutrals),
+                #                    ngrid=ngrid)
+                #end
+                @testset "all" begin
+                    global_rank[] == 0 && println("\nall terms")
+                    testconvergence(input, setup_advance(:all, neutrals),
+                                    ngrid=ngrid)
+                end
+            end
+        end
+    end
+
+    return nothing
+end
+
+function runtests_return(which_term::Symbol=:all; ngrid=nothing, neutrals=true)
+    global_rank[] == 0 && println("MMS tests with return")
+
+    if neutrals
+        input = input_sound_wave_periodic
+    else
+        input = input_sound_wave_periodic_no_neutrals
+    end
+
+    results = testconvergence(input, setup_advance(which_term, neutrals);
+                              ngrid=ngrid, returnstuff=true)
+
+    return results
+end
+
+end # ManufacturedSolutionsTests
+
+
+using .ManufacturedDDTTests
+
+ManufacturedDDTTests.runtests()
diff --git a/test/manufactured_solutions_tests.jl b/test/manufactured_solutions_tests.jl
new file mode 100644
index 000000000..d4f88cac7
--- /dev/null
+++ b/test/manufactured_solutions_tests.jl
@@ -0,0 +1,224 @@
+"""
+Test cases using the method of manufactured solutions (MMS)
+"""
+module ManufacturedSolutionsTests
+
+include("setup.jl")
+include("mms_utils.jl")
+
+using moment_kinetics.post_processing: L2_error_norm, L_infinity_error_norm
+using moment_kinetics.manufactured_solns
+using moment_kinetics.type_definitions
+
+# Create a temporary directory for test output
+test_output_directory = tempname()
+mkpath(test_output_directory)
+
+const input_sound_wave_periodic = Dict(
+    "use_manufactured_solns" => true,
+    "n_ion_species" => 1,
+    "n_neutral_species" => 0,
+    "boltzmann_electron_response" => true,
+    "run_name" => "MMS-rperiodic",
+    "base_directory" => test_output_directory,
+    "evolve_moments_density" => false,
+    "evolve_moments_parallel_flow" => false,
+    "evolve_moments_parallel_pressure" => false,
+    "evolve_moments_conservation" => false,
+    "T_e" => 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.62831853071,
+    "ionization_frequency" => 0.0,
+    #"nstep" => 10, #1700,
+    #"dt" => 0.002,
+    #"nwrite" => 10, #1700,
+    "nstep" => 1700, #1700,
+    "dt" => 0.0002, #0.002,
+    "nwrite" => 1700, #1700,
+    "use_semi_lagrange" => false,
+    "n_rk_stages" => 4,
+    "split_operators" => false,
+    "z_ngrid" => 4,
+    "z_nelement" => 2,
+    "z_bc" => "periodic",
+    "z_discretization" => "chebyshev_pseudospectral",
+    "r_ngrid" => 4,
+    "r_nelement" => 2,
+    "r_bc" => "periodic",
+    "r_discretization" => "chebyshev_pseudospectral",
+    "vpa_ngrid" => 4,
+    "vpa_nelement" => 4,
+    "vpa_L" => 8.0,
+    "vpa_bc" => "periodic",
+    "vpa_discretization" => "chebyshev_pseudospectral",
+    "vperp_ngrid" => 4,
+    "vperp_nelement" => 4,
+    "vperp_L" => 8.0,
+    "vperp_bc" => "periodic",
+    "vperp_discretization" => "chebyshev_pseudospectral",
+)
+
+"""
+    runcase(input::Dict)
+
+Run a simulation with parameters set by `input` using manufactured sources and return
+the errors in each variable compared to the manufactured solution.
+"""
+function runcase(input::Dict)
+    in_manf, phi_manf, f_manf = nothing, nothing, nothing
+    # call setup_moment_kinetics(), time_advance!(), cleanup_moment_kinetics!()
+    # separately so we can run manufactured_solutions_as_arrays() in parallel
+    quietoutput() do
+        # run simulation
+        pdf, vz, vr, vzeta, vpa, vperp, z, r, spectral_objects, composition,
+        drive_input, moments, t_input, collisions, species, geometry,
+        boundary_distributions = setup_moment_kinetics(input_dict)
+        time_advance!(pdf, vz, vr, vzeta, vpa, vperp, z, r, spectral_objects,
+                      composition, drive_input, moments, t_input, collisions,
+                      species, geometry, boundary_distributions)
+
+        n_manf, phi_manf, f_manf =
+            manufactured_solutions_as_arrays(t, r, z, vperp, vpa)
+
+        if global_rank[] == 0
+            # Need to copy as cleanup_moment_kinetics!() will invalidate
+            # shared-memory arrays
+            n_manf = copy(n_manf)
+            phi_manf = copy(phi_manf)
+            f_manf = copy(f_manf)
+        end
+
+        cleanup_moment_kinetics!(io, cdf)
+    end
+
+    n_error_2 = nothing
+    n_error_inf = nothing
+    phi_error_2 = nothing
+    phi_error_inf = nothing
+    f_error_2 = nothing
+    f_error_inf = nothing
+    if global_rank[] == 0
+        output = load_test_output(input, (:phi, :moments, :f))
+
+        t = output["time"][end]
+        n = output["density"][:,:,1,end]
+        phi = output["phi"][:,:,end]
+        f = output["f"][:,:,:,:,1,end]
+        f0 = f[size(f,1)÷2, 1, :, :]
+
+        f0_manf = f_manf[size(f,1)÷2, 1, :, :]
+
+        n_error_2 = L2_error_norm(n, n_manf)
+        n_error_inf = L_infinity_error_norm(n, n_manf)
+
+        phi_error_2 = L2_error_norm(phi, phi_manf)
+        phi_error_inf = L_infinity_error_norm(phi, phi_manf)
+
+        f_error_2 = L2_error_norm(f, f_manf)
+        f_error_inf = L_infinity_error_norm(f, f_manf)
+
+        f0_error_2 = L2_error_norm(f0, f0_manf)
+        f0_error_inf = L_infinity_error_norm(f0, f0_manf)
+
+        println("n ", n_error_2, " ", n_error_inf)
+        println("phi ", phi_error_2, " ", phi_error_inf)
+        println("f ", f_error_2, " ", f_error_inf)
+        println("f0 ", f0_error_2, " ", f0_error_inf)
+    end
+
+    return n_error_2, n_error_inf, phi_error_2, phi_error_inf, f_error_2, f_error_inf
+end
+
+"""
+    testconvergence(input::Dict)
+
+Test convergence with spatial resolution
+
+The parameters for the run are given in `input::Dict`.
+"""
+function testconvergence(input::Dict)
+    n_errors_2 = Vector{mk_float}(undef, 0)
+    n_errors_inf = Vector{mk_float}(undef, 0)
+    phi_errors_2 = Vector{mk_float}(undef, 0)
+    phi_errors_inf = Vector{mk_float}(undef, 0)
+    f_errors_2 = Vector{mk_float}(undef, 0)
+    f_errors_inf = Vector{mk_float}(undef, 0)
+
+    ngrid = get_and_check_ngrid(input)
+
+    nelement_values = [2, 3, 4]
+    for nelement ∈ nelement_values
+        global_rank[] == 0 && println("testing nelement=$nelement")
+        case_input = increase_resolution(input, nelement)
+
+        n_error_2, n_error_inf, phi_error_2, phi_error_inf, f_error_2,
+        f_error_inf = runcase(case_input)
+
+        if global_rank[] == 0
+            push!(n_errors_2, n_error_2)
+            push!(n_errors_inf, n_error_inf)
+            push!(phi_errors_2, phi_error_2)
+            push!(phi_errors_inf, phi_error_inf)
+            push!(f_errors_2, f_error_2)
+            push!(f_errors_inf, f_error_inf)
+        end
+    end
+
+    if global_rank[] == 0
+        n_convergence_2 = n_errors_2[1] ./ n_errors_2[2:end]
+        n_convergence_inf = n_errors_inf[1] ./ n_errors_inf[2:end]
+        phi_convergence_2 = phi_errors_2[1] ./ phi_errors_2[2:end]
+        phi_convergence_inf = phi_errors_inf[1] ./ phi_errors_inf[2:end]
+        f_convergence_2 = f_errors_2[1] ./ f_errors_2[2:end]
+        f_convergence_inf = f_errors_inf[1] ./ f_errors_inf[2:end]
+        expected_convergence = @. (nelement_values[2:end] / nelement_values[1])^(ngrid - 1)
+        println("n convergence")
+        println(n_convergence_2)
+        println(n_convergence_inf)
+        println("phi convergence")
+        println(phi_convergence_2)
+        println(phi_convergence_inf)
+        println("f convergence")
+        println(f_convergence_2)
+        println(f_convergence_inf)
+        println("expected convergence")
+        println(expected_convergence)
+    end
+end
+
+function runtests()
+    @testset "MMS" verbose=use_verbose begin
+        global_rank[] == 0 && println("MMS tests")
+
+        @testset "r-periodic, z-periodic" begin
+            testconvergence(input_sound_wave_periodic)
+        end
+    end
+
+    return nothing
+end
+
+end # ManufacturedSolutionsTests
+
+
+using .ManufacturedSolutionsTests
+
+ManufacturedSolutionsTests.runtests()
diff --git a/test/mms_utils.jl b/test/mms_utils.jl
new file mode 100644
index 000000000..d574b60f4
--- /dev/null
+++ b/test/mms_utils.jl
@@ -0,0 +1,93 @@
+"""
+Some shared functions used by MMS tests
+"""
+module MMSTestUtils
+
+export increase_resolution, get_and_check_ngrid, set_ngrid, test_error_series
+
+using moment_kinetics.type_definitions
+
+"""
+    increase_resolution(input::Dict, factor)
+
+Increase resolution of simulation by multiplying the numbers of elements `*_nelement` in
+the `input` settings by `factor`.
+"""
+function increase_resolution(input::Dict, nelement)
+    result = copy(input)
+    result["run_name"] = input["run_name"] * "_$nelement"
+    for key ∈ keys(result)
+        if occursin("_nelement", key)
+            if occursin("v", key) || occursin("gyrophase", key)
+                result[key] = 2 * nelement
+            else
+                result[key] = nelement
+            end
+        end
+    end
+
+    return result
+end
+
+"""
+    get_and_check_ngrid(input::Dict)
+
+Get value of `ngrid` and check that it is the same for all dimensions. `ngrid` needs to
+be the same as it sets the convergence order, and we want this to be the same for all
+operators.
+"""
+function get_and_check_ngrid(input::Dict)::mk_int
+    ngrid = nothing
+
+    for key ∈ keys(input)
+        if occursin("_ngrid", key)
+            if ngrid === nothing
+                ngrid = input[key]
+            else
+                if ngrid != input[key]
+                    error("*_ngrid should all be the same, but $key=$(input[key]) when "
+                          * "we already found ngrid=$ngrid")
+                end
+            end
+        end
+    end
+
+    return ngrid
+end
+
+"""
+    set_ngrid(input::Dict, ngrid::mk_int)
+
+Set value of `ngrid`, the same for all dimensions.
+"""
+function set_ngrid(input::Dict, ngrid::mk_int)
+    for key ∈ keys(input)
+        if occursin("_ngrid", key)
+            input[key] = ngrid
+        end
+    end
+
+    return nothing
+end
+
+"""
+    test_error_series(errors::Vector{mk_float}, resolution_factors::Vector,
+                      expected_order, expected_lowest)
+
+Test whether the error norms in `errors` converge as expected with increases in
+resolution by `resolution_factors`. `expected_order` is the order p such that the error
+is expected to be proportional to h^p. `expected_lowest` is the expected value of the
+error at the lowest resolution (used as a regression test).
+
+Note the entries in `errors` and `resolution_factors` should be sorted in increasing
+order of `resolution_factors`.
+"""
+function test_error_series(errors::Vector{mk_float}, resolution_factors::Vector,
+                           expected_order, expected_lowest)
+    error_factors = errors[1:end-1] ./ errors[2:end]
+    expected_factors = resolution_factors[2:end].^expected_order
+end
+
+end # MMSTestUtils
+
+using .MMSTestUtils
diff --git a/test/nonlinear_sound_wave_tests.jl b/test/nonlinear_sound_wave_tests.jl
index 7240d469f..ef7e935a0 100644
--- a/test/nonlinear_sound_wave_tests.jl
+++ b/test/nonlinear_sound_wave_tests.jl
@@ -5,7 +5,6 @@ include("setup.jl")
 using Base.Filesystem: tempname
 using TimerOutputs
 
-using moment_kinetics.chebyshev: setup_chebyshev_pseudospectral
 using moment_kinetics.coordinates: define_coordinate
 using moment_kinetics.input_structs: grid_input, advection_input
 using moment_kinetics.load_data: open_netcdf_file, load_coordinate_data,
@@ -257,9 +256,6 @@ function run_test(test_input, rtol; args...)
             # open the netcdf file and give it the handle 'fid'
             fid = open_netcdf_file(path)
 
-            # load space-time coordinate data
-            nvpa, vpa, vpa_wgts, nz, z, z_wgts, Lz, nr, r, r_wgts, Lr, ntime, time, n_ion_species, n_neutral_species = load_coordinate_data(fid)
-
             # load fields data
             phi_zrt = load_fields_data(fid)
 
@@ -297,21 +293,11 @@ function run_test(test_input, rtol; args...)
                            z_L, test_input["z_discretization"], "",
                            "periodic", #test_input["z_bc"],
                            adv_input)
-        z = define_coordinate(input)
-        if test_input["z_discretization"] == "chebyshev_pseudospectral"
-            z_spectral = setup_chebyshev_pseudospectral(z)
-        else
-            z_spectral = false
-        end
+        z, z_spectral = define_coordinate(input)
         input = grid_input("coord", test_input["vpa_ngrid"], test_input["vpa_nelement"],
                            vpa_L, test_input["vpa_discretization"], "",
                            test_input["vpa_bc"], adv_input)
-        vpa = define_coordinate(input)
-        if test_input["vpa_discretization"] == "chebyshev_pseudospectral"
-            vpa_spectral = setup_chebyshev_pseudospectral(vpa)
-        else
-            vpa_spectral = false
-        end
+        vpa, vpa_spectral = define_coordinate(input)
 
         # Test against values interpolated onto 'expected' grid which is fairly coarse no we
         # do not have to save too much data in this file
diff --git a/test/setup.jl b/test/setup.jl
index 413978622..3ac1272f4 100644
--- a/test/setup.jl
+++ b/test/setup.jl
@@ -12,10 +12,14 @@ using moment_kinetics
 
 module MKTestUtilities
 
-export use_verbose, @long, quietoutput, global_rank, maxabs_norm, @testset_skip
+export use_verbose, @long, quietoutput, global_rank, load_test_output, maxabs_norm,
+       @testset_skip
 
 using moment_kinetics.communication: global_rank
 using moment_kinetics.command_line_options: get_options
+using moment_kinetics.load_data: open_netcdf_file
+using moment_kinetics.load_data: load_coordinate_data, load_fields_data,
+                                 load_charged_particle_moments_data, load_pdf_data
 
 const use_verbose = get_options()["verbose"]
 
@@ -119,6 +123,49 @@ macro testset_skip(args...)
     return ex
 end
 
+"""
+    load_test_output(input::Dict, to_load::Tuple{Symbol})
+
+Load the output of a test that was run with settings in `input`. `to_load` specifies
+which variables to load - it can include potential `:phi`, moments `:moments`,
+distribution function `:f`. Coordinate data is always loaded.
+
+Returns a `Dict` whose keys are `String` containing all loaded data.
+"""
+function load_test_output(input::Dict, to_load::Tuple{Vararg{Symbol}})
+    output = Dict{String, Any}()
+
+    path = joinpath(realpath(input["base_directory"]), input["run_name"], input["run_name"])
+
+    # open the netcdf file and give it the handle 'fid'
+    fid = open_netcdf_file(path)
+
+    # load space-time coordinate data
+    output["vpa"], output["vperp"], output["z"], output["r"], output["ntime"],
+    output["time"] = load_coordinate_data(fid)
+
+    if :phi ∈ to_load
+        # Load EM fields
+        output["phi"] = load_fields_data(fid)
+    end
+
+    if :moments ∈ to_load
+        # Load moments
+        output["density"], output["parallel_flow"], output["parallel_pressure"],
+        output["parallel_heat_flux"], output["thermal_speed"], output["n_species"],
+        output["evolve_ppar"] = load_charged_particle_moments_data(fid)
+    end
+
+    if :f ∈ to_load
+        # load full (vpa,vperp,z,r,species,t) particle distribution function (pdf) data
+        output["f"] = load_pdf_data(fid)
+    end
+
+    close(fid)
+
+    return output
+end
+
 end
 
 using .MKTestUtilities
diff --git a/test/sound_wave_tests.jl b/test/sound_wave_tests.jl
index 5bb16b4f1..a2f07e0c2 100644
--- a/test/sound_wave_tests.jl
+++ b/test/sound_wave_tests.jl
@@ -7,8 +7,6 @@ using TimerOutputs
 #using Plots: plot, plot!, gui
 
 using moment_kinetics.array_allocation: allocate_float
-using moment_kinetics.load_data: open_netcdf_file
-using moment_kinetics.load_data: load_coordinate_data, load_fields_data
 using moment_kinetics.analysis: analyze_fields_data
 using moment_kinetics.post_processing: fit_delta_phi_mode
 
@@ -162,44 +160,38 @@ function run_test(test_input, analytic_frequency, analytic_growth_rate,
             # Load and analyse output
             #########################
 
-            path = joinpath(realpath(input["base_directory"]), name, name)
+            output = load_test_output(input, (:phi,))
 
-            # open the netcdf file and give it the handle 'fid'
-            fid = open_netcdf_file(path)
+            z = output["z"]
+            ntime = output["ntime"]
+            time = output["time"]
 
-            # load space-time coordinate data
-            nvpa, vpa, vpa_wgts, nvperp, vperp, vperp_wgts,
-            nz, z, z_wgts, Lz, nr, r, r_wgts, Lr, ntime, time, n_ion_species, n_neutral_species = load_coordinate_data(fid)
+            phi_zrt = output["phi"]
 
-            # load fields data
-            phi_zrt = load_fields_data(fid)
-
-            close(fid)
-            
             ir0 = 1 
             
             phi = phi_zrt[:,ir0,:]
 
             # analyze the fields data
-            phi_fldline_avg, delta_phi = analyze_fields_data(phi, ntime, nz, z_wgts, Lz)
+            phi_fldline_avg, delta_phi = analyze_fields_data(phi, ntime, z.n, z.wgts, z.L)
 
             # use a fit to calculate the damping rate and growth rate of the perturbed
             # electrostatic potential
             itime_max = ntime
-            iz0 = cld(nz, 3)
+            iz0 = cld(z.n, 3)
             shifted_time = allocate_float(ntime)
             @. shifted_time = time - time[itime_min]
-            @views phi_fit = fit_delta_phi_mode(shifted_time[itime_min:itime_max], z,
+            @views phi_fit = fit_delta_phi_mode(shifted_time[itime_min:itime_max], z.grid,
                                                 delta_phi[:, itime_min:itime_max])
             ## The following plot code (copied from post_processing.jl) may be helpful for
             ## debugging tests. Uncomment to use, and also uncomment
             ## `using Plots: plot, plot!, gui at the top of the file.
-            #L = z[end] - z[begin]
+            #L = z.grid[end] - z.grid[begin]
             #fitted_delta_phi =
-            #    @. (phi_fit.amplitude0 * cos(2.0 * π * (z[iz0] + phi_fit.offset0) / L)
+            #    @. (phi_fit.amplitude0 * cos(2.0 * π * (z.grid[iz0] + phi_fit.offset0) / L)
             #        * exp(phi_fit.growth_rate * shifted_time)
             #        * cos(phi_fit.frequency * shifted_time + phi_fit.phase))
-            #@views plot(time, abs.(delta_phi[iz0,:]), xlabel="t*Lz/vti", ylabel="δϕ", yaxis=:log)
+            #@views plot(time, abs.(delta_phi[iz0,:]), xlabel="t*z.L/vti", ylabel="δϕ", yaxis=:log)
             #plot!(time, abs.(fitted_delta_phi))
             #gui()
         end
diff --git a/test/wall_bc_tests.jl b/test/wall_bc_tests.jl
index 5dfd7415a..b52214e6d 100644
--- a/test/wall_bc_tests.jl
+++ b/test/wall_bc_tests.jl
@@ -8,7 +8,6 @@ include("setup.jl")
 using Base.Filesystem: tempname
 using TimerOutputs
 
-using moment_kinetics.chebyshev: setup_chebyshev_pseudospectral
 using moment_kinetics.coordinates: define_coordinate
 using moment_kinetics.input_structs: grid_input, advection_input
 using moment_kinetics.interpolation: interpolate_to_grid_z
@@ -153,9 +152,6 @@ function run_test(test_input, expected_phi, tolerance; args...)
             # open the netcdf file and give it the handle 'fid'
             fid = open_netcdf_file(path)
 
-            # load space-time coordinate data
-            nvpa, vpa, vpa_wgts, nz, z, z_wgts, Lz, nr, r, r_wgts, Lr, ntime, time, n_ion_species, n_neutral_species = load_coordinate_data(fid)
-
             # load fields data
             phi_zrt = load_fields_data(fid)
 
@@ -182,12 +178,7 @@ function run_test(test_input, expected_phi, tolerance; args...)
         input = grid_input("coord", test_input["z_ngrid"], test_input["z_nelement"], 1.0,
                            test_input["z_discretization"], "", test_input["z_bc"],
                            adv_input)
-        z = define_coordinate(input)
-        if test_input["z_discretization"] == "chebyshev_pseudospectral"
-            z_spectral = setup_chebyshev_pseudospectral(z)
-        else
-            z_spectral = false
-        end
+        z, z_spectral = define_coordinate(input)
 
         # Cross comparison of all discretizations to same benchmark
         phi_interp = interpolate_to_grid_z(cross_compare_points, phi[:, end], z, z_spectral)