From 8f120810bc63c8e530afb1b6cc07792f6d877fff Mon Sep 17 00:00:00 2001
From: Michael Hardman <michael.hardman@physics.ox.ac.uk>
Date: Fri, 8 Sep 2023 17:03:44 +0100
Subject: [PATCH 01/12] Attempt to start making a nonuniformly spaced elemental
 grid

---
 src/coordinates.jl | 52 ++++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 50 insertions(+), 2 deletions(-)

diff --git a/src/coordinates.jl b/src/coordinates.jl
index 704f67023..7e4cb04f5 100644
--- a/src/coordinates.jl
+++ b/src/coordinates.jl
@@ -4,6 +4,7 @@ module coordinates
 
 export define_coordinate, write_coordinate
 export equally_spaced_grid
+export set_element_boundaries
 
 using ..type_definitions: mk_float, mk_int
 using ..array_allocation: allocate_float, allocate_int
@@ -84,6 +85,12 @@ struct coordinate
     local_io_range::UnitRange{Int64}
     # global range to write into in output file
     global_io_range::UnitRange{Int64}
+    # scale for each element
+    element_scale::Array{mk_float,1}
+    # shift for each element
+    element_shift::Array{mk_float,1}
+    # boundaries for each element
+    element_boundaries::Array{mk_float,1}
 end
 
 """
@@ -105,6 +112,14 @@ function define_coordinate(input, parallel_io::Bool=false)
     # obtain (local) index mapping from the grid within each element
     # to the full grid
     imin, imax = elemental_to_full_grid_map(input.ngrid, input.nelement_local)
+    # initialise the data used to construct the grid
+    # boundaries for each element
+    element_boundaries = allocate_float(nelement_global)
+    element_boundaries_local = allocate_float(nelement_local)
+    
+    element_scale = allocate_float(nelement_local)
+    # shift for each element
+    element_shift = allocate_float(nelement_local)
     # initialize the grid and the integration weights associated with the grid
     # also obtain the Chebyshev theta grid and spacing if chosen as discretization option
     grid, wgts, uniform_grid = init_grid(input.ngrid, input.nelement_global,
@@ -126,7 +141,6 @@ function define_coordinate(input, parallel_io::Bool=false)
     # endpoints, so only two pieces of information must be shared
     send_buffer = allocate_float(1)
     receive_buffer = allocate_float(1)
-
     # Add some ranges to support parallel file io
     if !parallel_io
         # No parallel io, just write everything
@@ -149,7 +163,7 @@ function define_coordinate(input, parallel_io::Bool=false)
         cell_width, igrid, ielement, imin, imax, input.discretization, input.fd_option,
         input.bc, wgts, uniform_grid, duniform_dgrid, scratch, copy(scratch), copy(scratch),
         scratch_2d, copy(scratch_2d), advection, send_buffer, receive_buffer, input.comm,
-        local_io_range, global_io_range)
+        local_io_range, global_io_range, element_scale, element_shift, element_boundaries)
 
     if input.discretization == "chebyshev_pseudospectral" && coord.n > 1
         # create arrays needed for explicit Chebyshev pseudospectral treatment in this
@@ -167,6 +181,40 @@ function define_coordinate(input, parallel_io::Bool=false)
     return coord, spectral
 end
 
+function set_element_boundaries(nelement_global, L, element_spacing_option)
+    # set global element boundaries
+    element_boundaries = allocate_float(nelement_global+1)
+    if element_spacing_option == "sqrt"
+        # number of boundaries of sqrt grid
+        nsqrt = floor(mk_int,(nelement_global+1)/2)
+        println("nsqrt",nsqrt)
+        # number of boundaries of uniform grid
+        nuniform = nelement_global + 3 - 2*nsqrt
+        println("nuniform",nuniform)
+        DL = L/6.0 # 1/3 of the domain is uniformly spaced
+        delta = 2.0*DL/(nuniform-1) # length of each element in the uniform section
+        for j in 1:nsqrt
+            element_boundaries[j] = -(L/2.0) + ((L/2.0) - DL)*((j-1)/(nsqrt-1))^2
+        end
+        println(element_boundaries)
+        for j in 2:nuniform-1 #nsqrt+1:nelement_global + 1 - nsqrt
+            element_boundaries[nsqrt-1+j] = -DL + delta*(j-1) 
+        end
+        println(element_boundaries)
+        for j in 1:nsqrt
+            element_boundaries[(nelement_global+1)+ 1 - j] = (L/2.0) - ((L/2.0) - DL)*((j-1)/(nsqrt-1))^2
+        end
+        println(element_boundaries)
+    elseif element_spacing_option == "uniform" # uniform spacing 
+        for j in 1:nelement_global+1
+            element_boundaries[j] = L*((j-1)/(nelement_global) - 0.5)
+        end
+    else 
+        println("ERROR: element_spacing_option: ",element_spacing_option, " not supported")
+    end
+    return element_boundaries
+end
+
 """
 setup a grid with n_global grid points on the interval [-L/2,L/2]
 """

From 5475949b8c4602e2d0ed9c70d0267e2699fac32d Mon Sep 17 00:00:00 2001
From: Michael Hardman <michael.hardman@physics.ox.ac.uk>
Date: Fri, 8 Sep 2023 20:17:35 +0100
Subject: [PATCH 02/12] Version of element boundary routine working for
 nelement >= 2

---
 src/coordinates.jl | 44 ++++++++++++++++++++++++++++++--------------
 1 file changed, 30 insertions(+), 14 deletions(-)

diff --git a/src/coordinates.jl b/src/coordinates.jl
index 7e4cb04f5..d38a8e8af 100644
--- a/src/coordinates.jl
+++ b/src/coordinates.jl
@@ -114,12 +114,12 @@ function define_coordinate(input, parallel_io::Bool=false)
     imin, imax = elemental_to_full_grid_map(input.ngrid, input.nelement_local)
     # initialise the data used to construct the grid
     # boundaries for each element
-    element_boundaries = allocate_float(nelement_global)
-    element_boundaries_local = allocate_float(nelement_local)
+    element_boundaries = allocate_float(input.nelement_global)
+    element_boundaries_local = allocate_float(input.nelement_local)
     
-    element_scale = allocate_float(nelement_local)
+    element_scale = allocate_float(input.nelement_local)
     # shift for each element
-    element_shift = allocate_float(nelement_local)
+    element_shift = allocate_float(input.nelement_local)
     # initialize the grid and the integration weights associated with the grid
     # also obtain the Chebyshev theta grid and spacing if chosen as discretization option
     grid, wgts, uniform_grid = init_grid(input.ngrid, input.nelement_global,
@@ -186,23 +186,39 @@ function set_element_boundaries(nelement_global, L, element_spacing_option)
     element_boundaries = allocate_float(nelement_global+1)
     if element_spacing_option == "sqrt"
         # number of boundaries of sqrt grid
-        nsqrt = floor(mk_int,(nelement_global+1)/2)
+        nsqrt = floor(mk_int,(nelement_global)/2) + 1
+        if nelement_global%2 > 0 # odd
+            #fac = 0.0
+            #for j in 2:nsqrt
+            #    fac += 2.0*(((j-2)/(nsqrt-1))^2 - ((j-1)/(nsqrt-1))^2 )
+            #end
+            x = ((nsqrt-2)/(nsqrt-1))^2
+            if nsqrt < 3
+                fac = 2.0/3.0
+            else
+                #fac = 0.5*(sqrt( 1.0 + 4.0*x) - 1.0)/x
+                fac = 1.0/( 3.0/2.0 - 0.5*((nsqrt-2)/(nsqrt-1))^2)
+            end
+        else
+            fac = 1.0
+        end
+        
         println("nsqrt",nsqrt)
         # number of boundaries of uniform grid
-        nuniform = nelement_global + 3 - 2*nsqrt
-        println("nuniform",nuniform)
-        DL = L/6.0 # 1/3 of the domain is uniformly spaced
-        delta = 2.0*DL/(nuniform-1) # length of each element in the uniform section
+        #nuniform = nelement_global + 3 - 2*nsqrt
+        #println("nuniform",nuniform)
+        #DL = L/6.0 # 1/3 of the domain is uniformly spaced
+        #delta = 2.0*DL/(nuniform-1) # length of each element in the uniform section
         for j in 1:nsqrt
-            element_boundaries[j] = -(L/2.0) + ((L/2.0) - DL)*((j-1)/(nsqrt-1))^2
+            element_boundaries[j] = -(L/2.0) + fac*(L/2.0)*((j-1)/(nsqrt-1))^2
         end
         println(element_boundaries)
-        for j in 2:nuniform-1 #nsqrt+1:nelement_global + 1 - nsqrt
-            element_boundaries[nsqrt-1+j] = -DL + delta*(j-1) 
-        end
+        #for j in 2:nuniform-1 #nsqrt+1:nelement_global + 1 - nsqrt
+        #    element_boundaries[nsqrt-1+j] = -DL + delta*(j-1) 
+        #end
         println(element_boundaries)
         for j in 1:nsqrt
-            element_boundaries[(nelement_global+1)+ 1 - j] = (L/2.0) - ((L/2.0) - DL)*((j-1)/(nsqrt-1))^2
+            element_boundaries[(nelement_global+1)+ 1 - j] = (L/2.0) - fac*(L/2.0)*((j-1)/(nsqrt-1))^2
         end
         println(element_boundaries)
     elseif element_spacing_option == "uniform" # uniform spacing 

From e7788087bc52a346916e2386577c47c0c80ed53a Mon Sep 17 00:00:00 2001
From: Michael Hardman <michael.hardman@physics.ox.ac.uk>
Date: Mon, 11 Sep 2023 09:22:00 +0100
Subject: [PATCH 03/12] Merge the new method of element boundary specification
 into the latest master branch. Use the feature by specifying the option
 z_element_spacing_option = "sqrt". The default value is "uniform", which
 reproduces previous behaviour. The new method of specifying element
 boundaries gives the potential to optimize the element spacing for new
 applications in the future.

---
 src/chebyshev.jl                   | 30 ++++++++-------
 src/coordinates.jl                 | 56 ++++++++++++++--------------
 src/file_io.jl                     |  4 ++
 src/input_structs.jl               |  4 ++
 src/load_data.jl                   |  4 +-
 src/moment_kinetics_input.jl       | 59 +++++++++++++++++++-----------
 src/post_processing.jl             |  5 ++-
 test/calculus_tests.jl             | 56 +++++++++++++++++++---------
 test/interpolation_tests.jl        |  6 ++-
 test/nonlinear_sound_wave_tests.jl |  7 ++--
 test/wall_bc_tests.jl              |  5 ++-
 11 files changed, 145 insertions(+), 91 deletions(-)

diff --git a/src/chebyshev.jl b/src/chebyshev.jl
index 0908ac98e..854cdf7a5 100644
--- a/src/chebyshev.jl
+++ b/src/chebyshev.jl
@@ -42,8 +42,8 @@ end
 """
 initialize chebyshev grid scaled to interval [-box_length/2, box_length/2]
 """
-function scaled_chebyshev_grid(ngrid, nelement_global, nelement_local, n,
-			irank, box_length, imin, imax)
+function scaled_chebyshev_grid(ngrid, nelement_local, n,
+			element_scale, element_shift, imin, imax)
     # initialize chebyshev grid defined on [1,-1]
     # with n grid points chosen to facilitate
     # the fast Chebyshev transform (aka the discrete cosine transform)
@@ -55,7 +55,7 @@ function scaled_chebyshev_grid(ngrid, nelement_global, nelement_local, n,
     # setup the scale factor by which the Chebyshev grid on [-1,1]
     # is to be multiplied to account for the full domain [-L/2,L/2]
     # and the splitting into nelement elements with ngrid grid points
-    scale_factor = 0.5*box_length/float(nelement_global)
+    #scale_factor = 0.5*box_length/float(nelement_global)
     # account for the fact that the minimum index needed for the chebyshev_grid
     # within each element changes from 1 to 2 in going from the first element
     # to the remaining elements
@@ -63,8 +63,10 @@ function scaled_chebyshev_grid(ngrid, nelement_global, nelement_local, n,
     @inbounds for j ∈ 1:nelement_local
         #wgts[imin[j]:imax[j]] .= sqrt.(1.0 .- reverse(chebyshev_grid)[k:ngrid].^2) * scale_factor
         # amount by which to shift the centre of this element from zero
-        iel_global = j + irank*nelement_local
-        shift = box_length*((float(iel_global)-0.5)/float(nelement_global) - 0.5)
+        #iel_global = j + irank*nelement_local
+        #shift = box_length*((float(iel_global)-0.5)/float(nelement_global) - 0.5)
+        scale_factor = element_scale[j]
+        shift = element_shift[j]
         # reverse the order of the original chebyshev_grid (ran from [1,-1])
         # and apply the scale factor and shift
         grid[imin[j]:imax[j]] .= (reverse(chebyshev_grid)[k:ngrid] * scale_factor) .+ shift
@@ -72,7 +74,7 @@ function scaled_chebyshev_grid(ngrid, nelement_global, nelement_local, n,
         # to avoid double-counting boundary element
         k = 2
     end
-    wgts = clenshaw_curtis_weights(ngrid, nelement_local, n, imin, imax, scale_factor)
+    wgts = clenshaw_curtis_weights(ngrid, nelement_local, n, imin, imax, element_scale)
     return grid, wgts
 end
 
@@ -90,7 +92,7 @@ function elementwise_derivative!(coord, ff, chebyshev::chebyshev_info)
     @boundscheck nelement == size(df,2) && coord.ngrid == size(df,1) || throw(BoundsError(df))
     # note that one must multiply by 2*nelement/L to get derivative
     # in scaled coordinate
-    scale_factor = 2.0*float(coord.nelement_global)/coord.L
+    #scale_factor = 2.0*float(coord.nelement_global)/coord.L
 	# scale factor is (length of a single element/2)^{-1}
 	
     # variable k will be used to avoid double counting of overlapping point
@@ -112,7 +114,7 @@ function elementwise_derivative!(coord, ff, chebyshev::chebyshev_info)
         # and multiply by scaling factor needed to go
         # from Chebyshev z coordinate to actual z
         for i ∈ 1:coord.ngrid
-            df[i,j] *= scale_factor
+            df[i,j] /= coord.element_scale[j]
         end
         k = 1
     end
@@ -323,23 +325,25 @@ end
 returns wgts array containing the integration weights associated
 with all grid points for Clenshaw-Curtis quadrature
 """
-function clenshaw_curtis_weights(ngrid, nelement_local, n, imin, imax, scale_factor)
+function clenshaw_curtis_weights(ngrid, nelement_local, n, imin, imax, element_scale)
     # create array containing the integration weights
     wgts = zeros(mk_float, n)
     # calculate the modified Chebshev moments of the first kind
     μ = chebyshevmoments(ngrid)
+    # calculate the raw weights for a normalised grid on [-1,1]
+    w = clenshawcurtisweights(μ)
     @inbounds begin
         # calculate the weights within a single element and
         # scale to account for modified domain (not [-1,1])
-        wgts[1:ngrid] = clenshawcurtisweights(μ)*scale_factor
+        wgts[1:ngrid] = w*element_scale[1]
         if nelement_local > 1
             # account for double-counting of points at inner element boundaries
-            wgts[ngrid] *= 2
+            #wgts[ngrid] *= 2
             for j ∈ 2:nelement_local
-                wgts[imin[j]:imax[j]] .= wgts[2:ngrid]
+                wgts[imin[j]-1:imax[j]] .+= w*element_scale[j]
             end
             # remove double-counting of outer element boundary for last element
-            wgts[n] *= 0.5
+            #wgts[n] *= 0.5
         end
     end
     return wgts
diff --git a/src/coordinates.jl b/src/coordinates.jl
index d38a8e8af..b44a7b2ec 100644
--- a/src/coordinates.jl
+++ b/src/coordinates.jl
@@ -4,7 +4,9 @@ module coordinates
 
 export define_coordinate, write_coordinate
 export equally_spaced_grid
+# testing
 export set_element_boundaries
+export init_grid
 
 using ..type_definitions: mk_float, mk_int
 using ..array_allocation: allocate_float, allocate_int
@@ -90,7 +92,9 @@ struct coordinate
     # shift for each element
     element_shift::Array{mk_float,1}
     # boundaries for each element
-    element_boundaries::Array{mk_float,1}
+    #element_boundaries::Array{mk_float,1}
+    # option used to set up element spacing
+    element_spacing_option::String
 end
 
 """
@@ -114,16 +118,12 @@ function define_coordinate(input, parallel_io::Bool=false)
     imin, imax = elemental_to_full_grid_map(input.ngrid, input.nelement_local)
     # initialise the data used to construct the grid
     # boundaries for each element
-    element_boundaries = allocate_float(input.nelement_global)
-    element_boundaries_local = allocate_float(input.nelement_local)
-    
-    element_scale = allocate_float(input.nelement_local)
-    # shift for each element
-    element_shift = allocate_float(input.nelement_local)
+    element_boundaries = set_element_boundaries(input.nelement_global, input.L, input.element_spacing_option)
+    # shift and scale factors for each local element
+    element_scale, element_shift = set_element_scale_and_shift(input.nelement_global, input.nelement_local, input.irank, element_boundaries)
     # initialize the grid and the integration weights associated with the grid
     # also obtain the Chebyshev theta grid and spacing if chosen as discretization option
-    grid, wgts, uniform_grid = init_grid(input.ngrid, input.nelement_global,
-        input.nelement_local, n_global, n_local, input.irank, input.L,
+    grid, wgts, uniform_grid = init_grid(input.ngrid, input.nelement_local, n_global, n_local, input.irank, input.L, element_scale, element_shift,
         imin, imax, igrid, input.discretization, input.name)
     # calculate the widths of the cells between neighboring grid points
     cell_width = grid_spacing(grid, n_local)
@@ -163,7 +163,7 @@ function define_coordinate(input, parallel_io::Bool=false)
         cell_width, igrid, ielement, imin, imax, input.discretization, input.fd_option,
         input.bc, wgts, uniform_grid, duniform_dgrid, scratch, copy(scratch), copy(scratch),
         scratch_2d, copy(scratch_2d), advection, send_buffer, receive_buffer, input.comm,
-        local_io_range, global_io_range, element_scale, element_shift, element_boundaries)
+        local_io_range, global_io_range, element_scale, element_shift, input.element_spacing_option)#, element_boundaries)
 
     if input.discretization == "chebyshev_pseudospectral" && coord.n > 1
         # create arrays needed for explicit Chebyshev pseudospectral treatment in this
@@ -184,7 +184,7 @@ end
 function set_element_boundaries(nelement_global, L, element_spacing_option)
     # set global element boundaries
     element_boundaries = allocate_float(nelement_global+1)
-    if element_spacing_option == "sqrt"
+    if element_spacing_option == "sqrt" && nelement_global > 3
         # number of boundaries of sqrt grid
         nsqrt = floor(mk_int,(nelement_global)/2) + 1
         if nelement_global%2 > 0 # odd
@@ -203,25 +203,14 @@ function set_element_boundaries(nelement_global, L, element_spacing_option)
             fac = 1.0
         end
         
-        println("nsqrt",nsqrt)
-        # number of boundaries of uniform grid
-        #nuniform = nelement_global + 3 - 2*nsqrt
-        #println("nuniform",nuniform)
-        #DL = L/6.0 # 1/3 of the domain is uniformly spaced
-        #delta = 2.0*DL/(nuniform-1) # length of each element in the uniform section
         for j in 1:nsqrt
             element_boundaries[j] = -(L/2.0) + fac*(L/2.0)*((j-1)/(nsqrt-1))^2
         end
-        println(element_boundaries)
-        #for j in 2:nuniform-1 #nsqrt+1:nelement_global + 1 - nsqrt
-        #    element_boundaries[nsqrt-1+j] = -DL + delta*(j-1) 
-        #end
-        println(element_boundaries)
         for j in 1:nsqrt
             element_boundaries[(nelement_global+1)+ 1 - j] = (L/2.0) - fac*(L/2.0)*((j-1)/(nsqrt-1))^2
         end
-        println(element_boundaries)
-    elseif element_spacing_option == "uniform" # uniform spacing 
+        
+    elseif element_spacing_option == "uniform" || nelement_global < 4 # uniform spacing 
         for j in 1:nelement_global+1
             element_boundaries[j] = L*((j-1)/(nelement_global) - 0.5)
         end
@@ -231,10 +220,23 @@ function set_element_boundaries(nelement_global, L, element_spacing_option)
     return element_boundaries
 end
 
+function set_element_scale_and_shift(nelement_global, nelement_local, irank, element_boundaries)
+    element_scale = allocate_float(nelement_local)
+    element_shift = allocate_float(nelement_local)
+    
+    for j in 1:nelement_local
+        iel_global = j + irank*nelement_local
+        upper_boundary = element_boundaries[iel_global+1]
+        lower_boundary = element_boundaries[iel_global]
+        element_scale[j] = 0.5*(upper_boundary-lower_boundary)
+        element_shift[j] = 0.5*(upper_boundary+lower_boundary)
+    end
+    return element_scale, element_shift
+end
 """
 setup a grid with n_global grid points on the interval [-L/2,L/2]
 """
-function init_grid(ngrid, nelement_global, nelement_local, n_global, n_local, irank, L,
+function init_grid(ngrid, nelement_local, n_global, n_local, irank, L, element_scale, element_shift,
                    imin, imax, igrid, discretization, name)
     uniform_grid = equally_spaced_grid(n_global, n_local, irank, L)
     uniform_grid_shifted = equally_spaced_grid_shifted(n_global, n_local, irank, L)
@@ -251,7 +253,7 @@ function init_grid(ngrid, nelement_global, nelement_local, n_global, n_local, ir
     elseif discretization == "chebyshev_pseudospectral"
         if name == "vperp"
             # initialize chebyshev grid defined on [-L/2,L/2]
-            grid, wgts = scaled_chebyshev_grid(ngrid, nelement_global, nelement_local, n_local, irank, L, imin, imax)
+            grid, wgts = scaled_chebyshev_grid(ngrid, nelement_local, n_local, element_scale, element_shift, imin, imax)
             grid .= grid .+ L/2.0 # shift to [0,L] appropriate to vperp variable
             wgts = 2.0 .* wgts .* grid # to include 2 vperp in jacobian of integral
                                         # see note above on normalisation
@@ -262,7 +264,7 @@ function init_grid(ngrid, nelement_global, nelement_local, n_global, n_local, ir
             # needed to obtain Chebyshev spectral coefficients
             # 'wgts' are the integration weights attached to each grid points
             # that are those associated with Clenshaw-Curtis quadrature
-            grid, wgts = scaled_chebyshev_grid(ngrid, nelement_global, nelement_local, n_local, irank, L, imin, imax)
+            grid, wgts = scaled_chebyshev_grid(ngrid, nelement_local, n_local, element_scale, element_shift, imin, imax)
         end
     elseif discretization == "finite_difference"
         if name == "vperp"
diff --git a/src/file_io.jl b/src/file_io.jl
index bb53a90f2..5dd229843 100644
--- a/src/file_io.jl
+++ b/src/file_io.jl
@@ -521,6 +521,10 @@ function define_io_coordinate!(parent, coord, coord_name, description, parallel_
         write_single_value!(group, "bc", coord.bc; parallel_io=parallel_io,
                             description="boundary condition for $coord_name")
 
+        # write the boundary condition for the coordinate
+        write_single_value!(group, "element_spacing_option", coord.element_spacing_option; parallel_io=parallel_io,
+                            description="element_spacing_option for $coord_name")
+
         return group
     end
 
diff --git a/src/input_structs.jl b/src/input_structs.jl
index b013f703b..a08168e58 100644
--- a/src/input_structs.jl
+++ b/src/input_structs.jl
@@ -127,6 +127,8 @@ mutable struct grid_input_mutable
     bc::String
     # mutable struct containing advection speed options
     advection::advection_input_mutable
+    # string option determining boundary spacing
+    element_spacing_option::String
 end
 
 """
@@ -156,6 +158,8 @@ struct grid_input
     advection::advection_input
     # MPI communicator
     comm::MPI.Comm
+    # string option determining boundary spacing
+    element_spacing_option::String
 end
 
 """
diff --git a/src/load_data.jl b/src/load_data.jl
index 206873bbb..5507965b1 100644
--- a/src/load_data.jl
+++ b/src/load_data.jl
@@ -213,6 +213,7 @@ function load_coordinate_data(fid, name; printout=false)
     discretization = load_variable(coord_group, "discretization")
     fd_option = load_variable(coord_group, "fd_option")
     bc = load_variable(coord_group, "bc")
+    element_spacing_option = load_variable(coord_group, "element_spacing_option")
 
     nelement_local = nothing
     if n_local == 1 && ngrid == 1
@@ -225,11 +226,10 @@ function load_coordinate_data(fid, name; printout=false)
     else
         nelement_global = (n_global-1) ÷ (ngrid-1)
     end
-
     # Define input to create coordinate struct
     input = grid_input(name, ngrid, nelement_global, nelement_local, nrank, irank, L,
                        discretization, fd_option, bc, advection_input("", 0.0, 0.0, 0.0),
-                       MPI.COMM_NULL)
+                       MPI.COMM_NULL, element_spacing_option)
 
     coord, spectral = define_coordinate(input)
 
diff --git a/src/moment_kinetics_input.jl b/src/moment_kinetics_input.jl
index 6663d6e11..cf1cf4546 100644
--- a/src/moment_kinetics_input.jl
+++ b/src/moment_kinetics_input.jl
@@ -340,7 +340,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
     # determine the boundary condition to impose in r
     # supported options are "periodic" and "Dirichlet"
     r.bc = get(scan_input, "r_bc", "periodic")
-
+    r.element_spacing_option = get(scan_input, "r_element_spacing_option", "uniform")
+    
     # overwrite some default parameters related to the z grid
     # ngrid is number of grid points per element
     z.ngrid = get(scan_input, "z_ngrid", 9)
@@ -357,7 +358,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
     # determine the boundary condition to impose in z
     # supported options are "constant", "periodic" and "wall"
     z.bc = get(scan_input, "z_bc", "wall")
-
+    z.element_spacing_option = get(scan_input, "z_element_spacing_option", "uniform")
+    
     # overwrite some default parameters related to the vpa grid
     # ngrid is the number of grid points per element
     vpa.ngrid = get(scan_input, "vpa_ngrid", 17)
@@ -374,7 +376,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
     # supported options are "chebyshev_pseudospectral" and "finite_difference"
     vpa.discretization = get(scan_input, "vpa_discretization", "chebyshev_pseudospectral")
     vpa.fd_option = get(scan_input, "vpa_finite_difference_option", "third_order_upwind")
-
+    vpa.element_spacing_option = get(scan_input, "vpa_element_spacing_option", "uniform")
+    
     num_diss_params = setup_numerical_dissipation(
         get(scan_input, "numerical_dissipation", Dict{String,Any}()), true)
 
@@ -395,7 +398,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
     # determine the discretization option for the vperp grid
     # supported options are "finite_difference_vperp" "chebyshev_pseudospectral_vperp"
     vperp.discretization = get(scan_input, "vperp_discretization", "chebyshev_pseudospectral_vperp")
-
+    vperp.element_spacing_option = get(scan_input, "vperp_element_spacing_option", "uniform")
+    
     # overwrite some default parameters related to the gyrophase grid
     # ngrid is the number of grid points per element
     gyrophase.ngrid = get(scan_input, "gyrophase_ngrid", 17)
@@ -421,7 +425,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
         # determine the discretization option for the vz grid
         # supported options are "chebyshev_pseudospectral" and "finite_difference"
         vz.discretization = get(scan_input, "vz_discretization", vpa.discretization)
-
+        vz.element_spacing_option = get(scan_input, "vz_element_spacing_option", "uniform")
+    
         # overwrite some default parameters related to the vr grid
         # ngrid is the number of grid points per element
         vr.ngrid = get(scan_input, "vr_ngrid", 1)
@@ -437,7 +442,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
         # determine the discretization option for the vr grid
         # supported options are "chebyshev_pseudospectral" and "finite_difference"
         vr.discretization = get(scan_input, "vr_discretization", "chebyshev_pseudospectral")
-
+        vr.element_spacing_option = get(scan_input, "vr_element_spacing_option", "uniform")
+    
         # overwrite some default parameters related to the vzeta grid
         # ngrid is the number of grid points per element
         vzeta.ngrid = get(scan_input, "vzeta_ngrid", 1)
@@ -453,6 +459,7 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
         # determine the discretization option for the vzeta grid
         # supported options are "chebyshev_pseudospectral" and "finite_difference"
         vzeta.discretization = get(scan_input, "vzeta_discretization", "chebyshev_pseudospectral")
+        vzeta.element_spacing_option = get(scan_input, "vzeta_element_spacing_option", "uniform")
     end
 
     is_1V = (vperp.ngrid == vperp.nelement_global == 1 && vzeta.ngrid ==
@@ -489,37 +496,37 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
 	z_advection_immutable = advection_input(z.advection.option, z.advection.constant_speed,
         z.advection.frequency, z.advection.oscillation_amplitude)
     z_immutable = grid_input("z", z.ngrid, z.nelement_global, z.nelement_local, nrank_z, irank_z, z.L, 
-        z.discretization, z.fd_option, z.bc, z_advection_immutable, comm_sub_z)
+        z.discretization, z.fd_option, z.bc, z_advection_immutable, comm_sub_z, z.element_spacing_option)
     r_advection_immutable = advection_input(r.advection.option, r.advection.constant_speed,
         r.advection.frequency, r.advection.oscillation_amplitude)
     r_immutable = grid_input("r", r.ngrid, r.nelement_global, r.nelement_local, nrank_r, irank_r, r.L,
-        r.discretization, r.fd_option, r.bc, r_advection_immutable, comm_sub_r)
+        r.discretization, r.fd_option, r.bc, r_advection_immutable, comm_sub_r, r.element_spacing_option)
 	# for dimensions below which do not currently use distributed-memory MPI
 	# assign dummy values to nrank, irank and comm of coord struct
     vpa_advection_immutable = advection_input(vpa.advection.option, vpa.advection.constant_speed,
         vpa.advection.frequency, vpa.advection.oscillation_amplitude)
     vpa_immutable = grid_input("vpa", vpa.ngrid, vpa.nelement_global, vpa.nelement_local, 1, 0, vpa.L,
-        vpa.discretization, vpa.fd_option, vpa.bc, vpa_advection_immutable, MPI.COMM_NULL)
+        vpa.discretization, vpa.fd_option, vpa.bc, vpa_advection_immutable, MPI.COMM_NULL, vpa.element_spacing_option)
     vperp_advection_immutable = advection_input(vperp.advection.option, vperp.advection.constant_speed,
         vperp.advection.frequency, vperp.advection.oscillation_amplitude)
     vperp_immutable = grid_input("vperp", vperp.ngrid, vperp.nelement_global, vperp.nelement_local, 1, 0, vperp.L,
-        vperp.discretization, vperp.fd_option, vperp.bc, vperp_advection_immutable, MPI.COMM_NULL)
+        vperp.discretization, vperp.fd_option, vperp.bc, vperp_advection_immutable, MPI.COMM_NULL, vperp.element_spacing_option)
     gyrophase_advection_immutable = advection_input(gyrophase.advection.option, gyrophase.advection.constant_speed,
         gyrophase.advection.frequency, gyrophase.advection.oscillation_amplitude)
     gyrophase_immutable = grid_input("gyrophase", gyrophase.ngrid, gyrophase.nelement_global, gyrophase.nelement_local, 1, 0, gyrophase.L,
-        gyrophase.discretization, gyrophase.fd_option, gyrophase.bc, gyrophase_advection_immutable, MPI.COMM_NULL)
+        gyrophase.discretization, gyrophase.fd_option, gyrophase.bc, gyrophase_advection_immutable, MPI.COMM_NULL, gyrophase.element_spacing_option)
     vz_advection_immutable = advection_input(vz.advection.option, vz.advection.constant_speed,
         vz.advection.frequency, vz.advection.oscillation_amplitude)
     vz_immutable = grid_input("vz", vz.ngrid, vz.nelement_global, vz.nelement_local, 1, 0, vz.L,
-        vz.discretization, vz.fd_option, vz.bc, vz_advection_immutable, MPI.COMM_NULL)
+        vz.discretization, vz.fd_option, vz.bc, vz_advection_immutable, MPI.COMM_NULL, vz.element_spacing_option)
     vr_advection_immutable = advection_input(vr.advection.option, vr.advection.constant_speed,
         vr.advection.frequency, vr.advection.oscillation_amplitude)
     vr_immutable = grid_input("vr", vr.ngrid, vr.nelement_global, vr.nelement_local, 1, 0, vr.L,
-        vr.discretization, vr.fd_option, vr.bc, vr_advection_immutable, MPI.COMM_NULL)
+        vr.discretization, vr.fd_option, vr.bc, vr_advection_immutable, MPI.COMM_NULL, vr.element_spacing_option)
     vzeta_advection_immutable = advection_input(vzeta.advection.option, vzeta.advection.constant_speed,
         vzeta.advection.frequency, vzeta.advection.oscillation_amplitude)
     vzeta_immutable = grid_input("vzeta", vzeta.ngrid, vzeta.nelement_global, vzeta.nelement_local, 1, 0, vzeta.L,
-        vzeta.discretization, vzeta.fd_option, vzeta.bc, vzeta_advection_immutable, MPI.COMM_NULL)
+        vzeta.discretization, vzeta.fd_option, vzeta.bc, vzeta_advection_immutable, MPI.COMM_NULL, vzeta.element_spacing_option)
     
     species_charged_immutable = Array{species_parameters,1}(undef,n_ion_species)
     species_neutral_immutable = Array{species_parameters,1}(undef,n_neutral_species)
@@ -660,10 +667,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics)
     # mutable struct containing advection speed options/inputs for z
     advection_z = advection_input_mutable(advection_option_z, advection_speed_z,
         frequency_z, oscillation_amplitude_z)
+    element_spacing_option_z = "uniform"
     # create a mutable structure containing the input info related to the z grid
     z = grid_input_mutable("z", ngrid_z, nelement_global_z, nelement_local_z, L_z,
         discretization_option_z, finite_difference_option_z, boundary_option_z,
-        advection_z)
+        advection_z, element_spacing_option_z)
     #################### parameters related to the r grid ######################
     # ngrid_r is number of grid points per element
     ngrid_r = 1
@@ -700,10 +708,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics)
     # mutable struct containing advection speed options/inputs for r
     advection_r = advection_input_mutable(advection_option_r, advection_speed_r,
         frequency_r, oscillation_amplitude_r)
+    element_spacing_option_r = "uniform"
     # create a mutable structure containing the input info related to the r grid
     r = grid_input_mutable("r", ngrid_r, nelement_global_r, nelement_local_r, L_r,
         discretization_option_r, finite_difference_option_r, boundary_option_r,
-        advection_r)
+        advection_r, element_spacing_option_r)
     ############################################################################
     ################### parameters related to the vpa grid #####################
     # ngrid_vpa is the number of grid points per element
@@ -738,10 +747,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics)
     # mutable struct containing advection speed options/inputs for z
     advection_vpa = advection_input_mutable(advection_option_vpa, advection_speed_vpa,
         frequency_vpa, oscillation_amplitude_vpa)
+    element_spacing_option_vpa = "uniform"
     # create a mutable structure containing the input info related to the vpa grid
     vpa = grid_input_mutable("vpa", ngrid_vpa, nelement_vpa, nelement_vpa, L_vpa,
         discretization_option_vpa, finite_difference_option_vpa, boundary_option_vpa,
-        advection_vpa)
+        advection_vpa, element_spacing_option_vpa)
     ############################################################################
     ################### parameters related to the vperp grid #####################
     # ngrid_vperp is the number of grid points per element
@@ -775,10 +785,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics)
     # mutable struct containing advection speed options/inputs for z
     advection_vperp = advection_input_mutable(advection_option_vperp, advection_speed_vperp,
         frequency_vperp, oscillation_amplitude_vperp)
+    element_spacing_option_vperp = "uniform"
     # create a mutable structure containing the input info related to the vperp grid
     vperp = grid_input_mutable("vperp", ngrid_vperp, nelement_vperp, nelement_vperp, L_vperp,
         discretization_option_vperp, finite_difference_option_vperp, boundary_option_vperp,
-        advection_vperp)
+        advection_vperp, element_spacing_option_vperp)
     ############################################################################
     ################### parameters related to the gyrophase grid #####################
     # ngrid_gyrophase is the number of grid points per element
@@ -798,10 +809,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics)
     oscillation_amplitude_gyrophase = 1.0
     advection_gyrophase = advection_input_mutable(advection_option_gyrophase, advection_speed_gyrophase,
         frequency_gyrophase, oscillation_amplitude_gyrophase)
+    element_spacing_option_gyrophase = "uniform"
     # create a mutable structure containing the input info related to the gyrophase grid
     gyrophase = grid_input_mutable("gyrophase", ngrid_gyrophase, nelement_gyrophase, nelement_gyrophase, L_gyrophase,
         discretization_option_gyrophase, finite_difference_option_gyrophase, boundary_option_gyrophase,
-        advection_gyrophase)
+        advection_gyrophase, element_spacing_option_gyrophase)
     ############################################################################
     ################### parameters related to the vr grid #####################
     # ngrid_vr is the number of grid points per element
@@ -833,10 +845,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics)
     # mutable struct containing advection speed options/inputs for z
     advection_vr = advection_input_mutable(advection_option_vr, advection_speed_vr,
         frequency_vr, oscillation_amplitude_vr)
+    element_spacing_option_vr = "uniform"
     # create a mutable structure containing the input info related to the vr grid
     vr = grid_input_mutable("vr", ngrid_vr, nelement_vr, nelement_vr, L_vr,
         discretization_option_vr, finite_difference_option_vr, boundary_option_vr,
-        advection_vr)
+        advection_vr, element_spacing_option_vr)
     ############################################################################
     ################### parameters related to the vz grid #####################
     # ngrid_vz is the number of grid points per element
@@ -868,10 +881,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics)
     # mutable struct containing advection speed options/inputs for z
     advection_vz = advection_input_mutable(advection_option_vz, advection_speed_vz,
         frequency_vz, oscillation_amplitude_vz)
+    element_spacing_option_vz = "uniform"
     # create a mutable structure containing the input info related to the vz grid
     vz = grid_input_mutable("vz", ngrid_vz, nelement_vz, nelement_vz, L_vz,
         discretization_option_vz, finite_difference_option_vz, boundary_option_vz,
-        advection_vz)
+        advection_vz, element_spacing_option_vz)
     ############################################################################
     ################### parameters related to the vzeta grid #####################
     # ngrid_vzeta is the number of grid points per element
@@ -903,10 +917,11 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics)
     # mutable struct containing advection speed options/inputs for z
     advection_vzeta = advection_input_mutable(advection_option_vzeta, advection_speed_vzeta,
         frequency_vzeta, oscillation_amplitude_vzeta)
+    element_spacing_option_vzeta = "uniform"
     # create a mutable structure containing the input info related to the vzeta grid
     vzeta = grid_input_mutable("vzeta", ngrid_vzeta, nelement_vzeta, nelement_vzeta, L_vzeta,
         discretization_option_vzeta, finite_difference_option_vzeta, boundary_option_vzeta,
-        advection_vzeta)
+        advection_vzeta, element_spacing_option_vzeta)
     #############################################################################
     # define default values and create corresponding mutable structs holding
     # information about the composition of the species and their initial conditions
diff --git a/src/post_processing.jl b/src/post_processing.jl
index 69f1021f3..d2d1f7ef3 100644
--- a/src/post_processing.jl
+++ b/src/post_processing.jl
@@ -196,10 +196,11 @@ end
 function construct_global_zr_coords(r_local, z_local)
 
     function make_global_input(coord_local)
+        #element_spacing_option = "" # dummy value
         return grid_input(coord_local.name, coord_local.ngrid,
             coord_local.nelement_global, coord_local.nelement_global, 1, 0, coord_local.L,
             coord_local.discretization, coord_local.fd_option, coord_local.bc,
-            coord_local.advection, MPI.COMM_NULL)
+            coord_local.advection, MPI.COMM_NULL, coord_local.element_spacing_option)
     end
 
     r_global, r_global_spectral = define_coordinate(make_global_input(r_local))
@@ -1054,7 +1055,7 @@ function analyze_and_plot_data(prefix...; run_index=nothing)
         composition, species, collisions,
         geometry, drive_input, num_diss_params, manufactured_solns_input = input
 
-    if !is_1D1V
+    if !is_1D1V || true
         # make plots and animations of the phi, Ez and Er
         plot_charged_moments_2D(density, parallel_flow, parallel_pressure, time,
                                 z_global.grid, r_global.grid, iz0, ir0, n_ion_species,
diff --git a/test/calculus_tests.jl b/test/calculus_tests.jl
index c1bbdd483..eab6e882e 100644
--- a/test/calculus_tests.jl
+++ b/test/calculus_tests.jl
@@ -39,9 +39,11 @@ function runtests()
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
 				comm = MPI.COMM_NULL # dummy value 
+                element_spacing_option = "uniform" # dummy value
 				input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
-                    discretization, fd_option, bc, adv_input, comm)
+                    discretization, fd_option, bc, adv_input, comm,
+                    element_spacing_option)
                 # create the coordinate struct 'x'
                 x, spectral = define_coordinate(input)
                 # create array for the function f(x) to be differentiated/integrated
@@ -86,9 +88,11 @@ function runtests()
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
 				comm = MPI.COMM_NULL # dummy value 
-				input = grid_input("coord", ngrid, nelement,
+				element_spacing_option = "uniform" # dummy value
+                input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
-                    "finite_difference", fd_option, bc, adv_input, comm)
+                    "finite_difference", fd_option, bc, adv_input, comm,
+                    element_spacing_option)
                 # create the coordinate struct 'x'
                 x, spectral = define_coordinate(input)
 
@@ -132,10 +136,12 @@ function runtests()
                 nelement_local = nelement
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
-				comm = MPI.COMM_NULL # dummy value 
+				comm = MPI.COMM_NULL # dummy value
+                element_spacing_option = "uniform" # dummy value
 				input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
-                    "finite_difference", fd_option, bc, adv_input, comm)
+                    "finite_difference", fd_option, bc, adv_input, comm,
+                    element_spacing_option)
                 # create the coordinate struct 'x'
                 x, spectral = define_coordinate(input)
 
@@ -175,10 +181,12 @@ function runtests()
                 nelement_local = nelement
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
-				comm = MPI.COMM_NULL # dummy value 
+				comm = MPI.COMM_NULL # dummy value
+                element_spacing_option = "uniform" # dummy value
 				input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
-                    "finite_difference", fd_option, bc, adv_input, comm)
+                    "finite_difference", fd_option, bc, adv_input, comm,
+                    element_spacing_option)
                 # create the coordinate struct 'x'
                 x, spectral = define_coordinate(input)
 
@@ -226,10 +234,12 @@ function runtests()
                 nelement_local = nelement
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
-				comm = MPI.COMM_NULL # dummy value 
+				comm = MPI.COMM_NULL # dummy value
+                element_spacing_option = "uniform" # dummy value
 				input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
-                    "finite_difference", fd_option, bc, adv_input, comm)
+                    "finite_difference", fd_option, bc, adv_input, comm,
+                    element_spacing_option)
                # create the coordinate struct 'x'
                 x, spectral = define_coordinate(input)
 
@@ -440,10 +450,12 @@ function runtests()
                 nelement_local = nelement
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
-				comm = MPI.COMM_NULL # dummy value 
+				comm = MPI.COMM_NULL # dummy value
+                element_spacing_option = "uniform"
 				input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
-                    "chebyshev_pseudospectral", fd_option, bc, adv_input, comm)
+                    "chebyshev_pseudospectral", fd_option, bc, adv_input, comm,
+                    element_spacing_option)
                 # create the coordinate struct 'x' and info for derivatives, etc.
                 x, spectral = define_coordinate(input)
 
@@ -634,10 +646,12 @@ function runtests()
                 nelement_local = nelement
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
-				comm = MPI.COMM_NULL # dummy value 
+				comm = MPI.COMM_NULL # dummy value
+                element_spacing_option = "uniform"
 				input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
-                    "chebyshev_pseudospectral", fd_option, bc, adv_input, comm)
+                    "chebyshev_pseudospectral", fd_option, bc, adv_input, comm,
+                    element_spacing_option)
                 # create the coordinate struct 'x' and info for derivatives, etc.
                 x, spectral = define_coordinate(input)
 
@@ -676,10 +690,12 @@ function runtests()
                 nelement_local = nelement
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
-				comm = MPI.COMM_NULL # dummy value 
+				comm = MPI.COMM_NULL # dummy value
+                element_spacing_option = "uniform"
 				input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
-                    "chebyshev_pseudospectral", fd_option, bc, adv_input, comm)
+                    "chebyshev_pseudospectral", fd_option, bc, adv_input, comm,
+                    element_spacing_option)
                 # create the coordinate struct 'x' and info for derivatives, etc.
                 x, spectral = define_coordinate(input)
 
@@ -726,10 +742,12 @@ function runtests()
                 nelement_local = nelement
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
-				comm = MPI.COMM_NULL # dummy value 
+				comm = MPI.COMM_NULL # dummy value
+                element_spacing_option = "uniform"
 				input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
-                    "chebyshev_pseudospectral", fd_option, bc, adv_input, comm)
+                    "chebyshev_pseudospectral", fd_option, bc, adv_input, comm,
+                    element_spacing_option)
                 # create the coordinate struct 'x' and info for derivatives, etc.
                 x, spectral = define_coordinate(input)
 
@@ -940,9 +958,11 @@ function runtests()
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
 				comm = MPI.COMM_NULL # dummy value
+                element_spacing_option = "uniform"
 				input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
-                    "chebyshev_pseudospectral", fd_option, bc, adv_input, comm)
+                    "chebyshev_pseudospectral", fd_option, bc, adv_input, comm,
+                    element_spacing_option)
                 # create the coordinate struct 'x' and info for derivatives, etc.
                 x, spectral = define_coordinate(input)
 
diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl
index b40a55ff7..33cd0fbaa 100644
--- a/test/interpolation_tests.jl
+++ b/test/interpolation_tests.jl
@@ -35,10 +35,12 @@ function runtests()
 			nelement_local = nelement
 			nrank_per_block = 0 # dummy value
 			irank = 0 # dummy value
-			comm = MPI.COMM_NULL # dummy value 
+			comm = MPI.COMM_NULL # dummy value
+            element_spacing_option = "uniform"
 			input = grid_input("coord", ngrid, nelement,
 				nelement_local, nrank_per_block, irank, L,
-				discretization, fd_option, bc, adv_input, comm)
+				discretization, fd_option, bc, adv_input, comm,
+                element_spacing_option)
             # create the coordinate struct 'z'
             z, spectral = define_coordinate(input)
 
diff --git a/test/nonlinear_sound_wave_tests.jl b/test/nonlinear_sound_wave_tests.jl
index 2742760cd..703c71b21 100644
--- a/test/nonlinear_sound_wave_tests.jl
+++ b/test/nonlinear_sound_wave_tests.jl
@@ -348,17 +348,18 @@ function run_test(test_input, rtol, atol, upar_rtol=nothing; args...)
         adv_input = advection_input("default", 1.0, 0.0, 0.0)
         nrank_per_block = 0 # dummy value
 		irank = 0 # dummy value
-		comm = MPI.COMM_NULL # dummy value 
+		comm = MPI.COMM_NULL # dummy value
+        element_spacing_option = "uniform"
 		input = grid_input("coord", test_input["z_ngrid"], test_input["z_nelement"], 
                            test_input["z_nelement"], nrank_per_block, irank,
 						   z_L, test_input["z_discretization"], "",
                            "periodic", #test_input["z_bc"],
-                           adv_input,comm)
+                           adv_input,comm, element_spacing_option)
         z, z_spectral = define_coordinate(input)
         input = grid_input("coord", test_input["vpa_ngrid"], test_input["vpa_nelement"],
                            test_input["vpa_nelement"], nrank_per_block, irank,
 						   vpa_L, test_input["vpa_discretization"], "",
-                           test_input["vpa_bc"], adv_input, comm)
+                           test_input["vpa_bc"], adv_input, comm, element_spacing_option)
         vpa, vpa_spectral = define_coordinate(input)
 
         # Test against values interpolated onto 'expected' grid which is fairly coarse no we
diff --git a/test/wall_bc_tests.jl b/test/wall_bc_tests.jl
index 9a5708921..9b51063e2 100644
--- a/test/wall_bc_tests.jl
+++ b/test/wall_bc_tests.jl
@@ -191,11 +191,12 @@ function run_test(test_input, expected_phi, tolerance; args...)
         adv_input = advection_input("default", 1.0, 0.0, 0.0)
 		nrank_per_block = 0 # dummy value
 		irank = 0 # dummy value
-		comm = MPI.COMM_NULL # dummy value 
+		comm = MPI.COMM_NULL # dummy value
+        element_spacing_option = "uniform"
         input = grid_input("coord", test_input["z_ngrid"], test_input["z_nelement"], 
 						   test_input["z_nelement"], nrank_per_block, irank, 1.0,
                            test_input["z_discretization"], "", test_input["z_bc"],
-                           adv_input, comm)
+                           adv_input, comm, element_spacing_option)
         z, z_spectral = define_coordinate(input)
 
         # Cross comparison of all discretizations to same benchmark

From 50449bc99e9b28d87ecf581c8b9a944dc3db5594 Mon Sep 17 00:00:00 2001
From: mrhardman <michael.hardman@physics.ox.ac.uk>
Date: Wed, 13 Sep 2023 09:03:53 +0100
Subject: [PATCH 04/12] Update src/chebyshev.jl

Delete commented out code

Co-authored-by: John Omotani <john.omotani@ukaea.uk>
---
 src/chebyshev.jl | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/chebyshev.jl b/src/chebyshev.jl
index 854cdf7a5..b62eb5601 100644
--- a/src/chebyshev.jl
+++ b/src/chebyshev.jl
@@ -63,8 +63,6 @@ function scaled_chebyshev_grid(ngrid, nelement_local, n,
     @inbounds for j ∈ 1:nelement_local
         #wgts[imin[j]:imax[j]] .= sqrt.(1.0 .- reverse(chebyshev_grid)[k:ngrid].^2) * scale_factor
         # amount by which to shift the centre of this element from zero
-        #iel_global = j + irank*nelement_local
-        #shift = box_length*((float(iel_global)-0.5)/float(nelement_global) - 0.5)
         scale_factor = element_scale[j]
         shift = element_shift[j]
         # reverse the order of the original chebyshev_grid (ran from [1,-1])

From d76373c0f88f61a2cff20466caedb6f2bcec8ac5 Mon Sep 17 00:00:00 2001
From: Michael Hardman <michael.hardman@physics.ox.ac.uk>
Date: Wed, 13 Sep 2023 13:07:24 +0100
Subject: [PATCH 05/12] Address @johnomotani comments regarding commented out
 lines and documentation.

---
 src/chebyshev.jl       | 32 ++++++++++++++++++++------------
 src/coordinates.jl     | 14 ++------------
 src/file_io.jl         |  2 +-
 src/post_processing.jl |  3 +--
 4 files changed, 24 insertions(+), 27 deletions(-)

diff --git a/src/chebyshev.jl b/src/chebyshev.jl
index b62eb5601..2a235d39d 100644
--- a/src/chebyshev.jl
+++ b/src/chebyshev.jl
@@ -41,6 +41,22 @@ end
 
 """
 initialize chebyshev grid scaled to interval [-box_length/2, box_length/2]
+we no longer pass the box_length to this function, but instead pass precomputed
+arrays element_scale and element_shift that are needed to compute the grid.
+
+ngrid -- number of points per element (including boundary points)
+nelement_local -- number of elements in the local (distributed memory MPI) grid
+n -- total number of points in the local grid (excluding duplicate points)
+element_scale -- the scale factor in the transform from the coordinates 
+                 where the element limits are -1, 1 to the coordinate where
+                 the limits are Aj = coord.grid[imin[j]-1] and Bj = coord.grid[imax[j]]
+                 element_scale = 0.5*(Bj - Aj)
+element_shift -- the centre of the element in the extended grid coordinate
+                 element_shift = 0.5*(Aj + Bj)
+imin -- the array of minimum indices of each element on the extended grid.
+        By convention, the duplicated points are not included, so for element index j > 1
+        the lower boundary point is actually imin[j] - 1
+imax -- the array of maximum indices of each element on the extended grid.
 """
 function scaled_chebyshev_grid(ngrid, nelement_local, n,
 			element_scale, element_shift, imin, imax)
@@ -55,14 +71,12 @@ function scaled_chebyshev_grid(ngrid, nelement_local, n,
     # setup the scale factor by which the Chebyshev grid on [-1,1]
     # is to be multiplied to account for the full domain [-L/2,L/2]
     # and the splitting into nelement elements with ngrid grid points
-    #scale_factor = 0.5*box_length/float(nelement_global)
+    
     # account for the fact that the minimum index needed for the chebyshev_grid
     # within each element changes from 1 to 2 in going from the first element
     # to the remaining elements
     k = 1
     @inbounds for j ∈ 1:nelement_local
-        #wgts[imin[j]:imax[j]] .= sqrt.(1.0 .- reverse(chebyshev_grid)[k:ngrid].^2) * scale_factor
-        # amount by which to shift the centre of this element from zero
         scale_factor = element_scale[j]
         shift = element_shift[j]
         # reverse the order of the original chebyshev_grid (ran from [1,-1])
@@ -88,11 +102,9 @@ function elementwise_derivative!(coord, ff, chebyshev::chebyshev_info)
     # check array bounds
     @boundscheck nelement == size(chebyshev.f,2) || throw(BoundsError(chebyshev.f))
     @boundscheck nelement == size(df,2) && coord.ngrid == size(df,1) || throw(BoundsError(df))
-    # note that one must multiply by 2*nelement/L to get derivative
-    # in scaled coordinate
-    #scale_factor = 2.0*float(coord.nelement_global)/coord.L
-	# scale factor is (length of a single element/2)^{-1}
-	
+    # note that one must multiply by a coordinate transform factor 1/element_scale[j] 
+    # for each element j to get derivative on the extended grid
+    
     # variable k will be used to avoid double counting of overlapping point
     # at element boundaries (see below for further explanation)
     k = 0
@@ -335,13 +347,9 @@ function clenshaw_curtis_weights(ngrid, nelement_local, n, imin, imax, element_s
         # scale to account for modified domain (not [-1,1])
         wgts[1:ngrid] = w*element_scale[1]
         if nelement_local > 1
-            # account for double-counting of points at inner element boundaries
-            #wgts[ngrid] *= 2
             for j ∈ 2:nelement_local
                 wgts[imin[j]-1:imax[j]] .+= w*element_scale[j]
             end
-            # remove double-counting of outer element boundary for last element
-            #wgts[n] *= 0.5
         end
     end
     return wgts
diff --git a/src/coordinates.jl b/src/coordinates.jl
index b44a7b2ec..6ad0164b7 100644
--- a/src/coordinates.jl
+++ b/src/coordinates.jl
@@ -4,9 +4,7 @@ module coordinates
 
 export define_coordinate, write_coordinate
 export equally_spaced_grid
-# testing
 export set_element_boundaries
-export init_grid
 
 using ..type_definitions: mk_float, mk_int
 using ..array_allocation: allocate_float, allocate_int
@@ -91,8 +89,6 @@ struct coordinate
     element_scale::Array{mk_float,1}
     # shift for each element
     element_shift::Array{mk_float,1}
-    # boundaries for each element
-    #element_boundaries::Array{mk_float,1}
     # option used to set up element spacing
     element_spacing_option::String
 end
@@ -163,7 +159,7 @@ function define_coordinate(input, parallel_io::Bool=false)
         cell_width, igrid, ielement, imin, imax, input.discretization, input.fd_option,
         input.bc, wgts, uniform_grid, duniform_dgrid, scratch, copy(scratch), copy(scratch),
         scratch_2d, copy(scratch_2d), advection, send_buffer, receive_buffer, input.comm,
-        local_io_range, global_io_range, element_scale, element_shift, input.element_spacing_option)#, element_boundaries)
+        local_io_range, global_io_range, element_scale, element_shift, input.element_spacing_option)
 
     if input.discretization == "chebyshev_pseudospectral" && coord.n > 1
         # create arrays needed for explicit Chebyshev pseudospectral treatment in this
@@ -188,15 +184,9 @@ function set_element_boundaries(nelement_global, L, element_spacing_option)
         # number of boundaries of sqrt grid
         nsqrt = floor(mk_int,(nelement_global)/2) + 1
         if nelement_global%2 > 0 # odd
-            #fac = 0.0
-            #for j in 2:nsqrt
-            #    fac += 2.0*(((j-2)/(nsqrt-1))^2 - ((j-1)/(nsqrt-1))^2 )
-            #end
-            x = ((nsqrt-2)/(nsqrt-1))^2
             if nsqrt < 3
                 fac = 2.0/3.0
             else
-                #fac = 0.5*(sqrt( 1.0 + 4.0*x) - 1.0)/x
                 fac = 1.0/( 3.0/2.0 - 0.5*((nsqrt-2)/(nsqrt-1))^2)
             end
         else
@@ -210,7 +200,7 @@ function set_element_boundaries(nelement_global, L, element_spacing_option)
             element_boundaries[(nelement_global+1)+ 1 - j] = (L/2.0) - fac*(L/2.0)*((j-1)/(nsqrt-1))^2
         end
         
-    elseif element_spacing_option == "uniform" || nelement_global < 4 # uniform spacing 
+    elseif element_spacing_option == "uniform" || (element_spacing_option == "sqrt" && nelement_global < 4) # uniform spacing 
         for j in 1:nelement_global+1
             element_boundaries[j] = L*((j-1)/(nelement_global) - 0.5)
         end
diff --git a/src/file_io.jl b/src/file_io.jl
index 5dd229843..d77ae35dd 100644
--- a/src/file_io.jl
+++ b/src/file_io.jl
@@ -521,7 +521,7 @@ function define_io_coordinate!(parent, coord, coord_name, description, parallel_
         write_single_value!(group, "bc", coord.bc; parallel_io=parallel_io,
                             description="boundary condition for $coord_name")
 
-        # write the boundary condition for the coordinate
+        # write the element spacing option for the coordinate
         write_single_value!(group, "element_spacing_option", coord.element_spacing_option; parallel_io=parallel_io,
                             description="element_spacing_option for $coord_name")
 
diff --git a/src/post_processing.jl b/src/post_processing.jl
index d2d1f7ef3..5b3857526 100644
--- a/src/post_processing.jl
+++ b/src/post_processing.jl
@@ -196,7 +196,6 @@ end
 function construct_global_zr_coords(r_local, z_local)
 
     function make_global_input(coord_local)
-        #element_spacing_option = "" # dummy value
         return grid_input(coord_local.name, coord_local.ngrid,
             coord_local.nelement_global, coord_local.nelement_global, 1, 0, coord_local.L,
             coord_local.discretization, coord_local.fd_option, coord_local.bc,
@@ -1055,7 +1054,7 @@ function analyze_and_plot_data(prefix...; run_index=nothing)
         composition, species, collisions,
         geometry, drive_input, num_diss_params, manufactured_solns_input = input
 
-    if !is_1D1V || true
+    if !is_1D1V
         # make plots and animations of the phi, Ez and Er
         plot_charged_moments_2D(density, parallel_flow, parallel_pressure, time,
                                 z_global.grid, r_global.grid, iz0, ir0, n_ion_species,

From d81cf24dfd872a99ea3483e2917f84d333667f00 Mon Sep 17 00:00:00 2001
From: mrhardman <michael.hardman@physics.ox.ac.uk>
Date: Wed, 13 Sep 2023 16:57:09 +0100
Subject: [PATCH 06/12] Update src/chebyshev.jl

Remove out of date comment

Co-authored-by: John Omotani <john.omotani@ukaea.uk>
---
 src/chebyshev.jl | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/src/chebyshev.jl b/src/chebyshev.jl
index 2a235d39d..30209ec7b 100644
--- a/src/chebyshev.jl
+++ b/src/chebyshev.jl
@@ -68,9 +68,6 @@ function scaled_chebyshev_grid(ngrid, nelement_local, n,
     chebyshev_grid = chebyshevpoints(ngrid)
     # create array for the full grid
     grid = allocate_float(n)
-    # setup the scale factor by which the Chebyshev grid on [-1,1]
-    # is to be multiplied to account for the full domain [-L/2,L/2]
-    # and the splitting into nelement elements with ngrid grid points
     
     # account for the fact that the minimum index needed for the chebyshev_grid
     # within each element changes from 1 to 2 in going from the first element

From e2a2224db44f46fc037fc9c6f0eb0e88ac334c6c Mon Sep 17 00:00:00 2001
From: Michael Hardman <michael.hardman@physics.ox.ac.uk>
Date: Mon, 25 Sep 2023 09:23:35 +0100
Subject: [PATCH 07/12] Add "sqrt" grid to fundamental theorem of calculus
 tests with large error tolerance = 1.0e-2.

---
 test/calculus_tests.jl | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/test/calculus_tests.jl b/test/calculus_tests.jl
index eab6e882e..2adf15a07 100644
--- a/test/calculus_tests.jl
+++ b/test/calculus_tests.jl
@@ -14,7 +14,7 @@ function runtests()
         println("calculus tests")
         @testset "fundamental theorem of calculus" begin
             @testset "$discretization $ngrid $nelement" for
-                    discretization ∈ ("finite_difference", "chebyshev_pseudospectral"),
+                    (discretization, element_spacing_option, etol) ∈ (("finite_difference", "uniform", 1.0e-15), ("chebyshev_pseudospectral", "uniform", 1.0e-15), ("chebyshev_pseudospectral", "sqrt", 1.0e-2)),
                     ngrid ∈ (5,6,7,8,9,10), nelement ∈ (1, 2, 3, 4, 5)
 
                 if discretization == "finite_difference" && (ngrid - 1) * nelement % 2 == 1
@@ -26,7 +26,6 @@ function runtests()
                     continue
                 end
 
-                etol = 1.0e-15
                 # define inputs needed for the test
                 L = 6.0
                 bc = "periodic"
@@ -39,7 +38,7 @@ function runtests()
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
 				comm = MPI.COMM_NULL # dummy value 
-                element_spacing_option = "uniform" # dummy value
+                #element_spacing_option = "uniform" # dummy value
 				input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
                     discretization, fd_option, bc, adv_input, comm,
@@ -676,7 +675,7 @@ function runtests()
         end
 
         @testset "Chebyshev pseudospectral derivatives (4 argument), Neumann" verbose=false begin
-            @testset "$nelement $ngrid" for bc ∈ ("constant", "zero"),
+            @testset "$nelement $ngrid" for bc ∈ ("constant", "zero"), element_spacing_option ∈ ("uniform", "sqrt"),
                     nelement ∈ (1:5), ngrid ∈ (3:33)
 
                 # define inputs needed for the test
@@ -691,7 +690,7 @@ function runtests()
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
 				comm = MPI.COMM_NULL # dummy value
-                element_spacing_option = "uniform"
+                #element_spacing_option = "uniform"
 				input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
                     "chebyshev_pseudospectral", fd_option, bc, adv_input, comm,

From a61ff3776c361fba5cd70d7a298f0dfd79610163 Mon Sep 17 00:00:00 2001
From: Michael Hardman <michael.hardman@physics.ox.ac.uk>
Date: Mon, 25 Sep 2023 10:15:14 +0100
Subject: [PATCH 08/12] Include sqrt grid in calculus tests as suggested by
 @johnomotani.

---
 test/calculus_tests.jl | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/test/calculus_tests.jl b/test/calculus_tests.jl
index 2adf15a07..bdcd10453 100644
--- a/test/calculus_tests.jl
+++ b/test/calculus_tests.jl
@@ -690,8 +690,7 @@ function runtests()
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
 				comm = MPI.COMM_NULL # dummy value
-                #element_spacing_option = "uniform"
-				input = grid_input("coord", ngrid, nelement,
+                input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
                     "chebyshev_pseudospectral", fd_option, bc, adv_input, comm,
                     element_spacing_option)
@@ -727,7 +726,7 @@ function runtests()
         end
 
         @testset "Chebyshev pseudospectral derivatives upwinding (5 argument), Neumann" verbose=false begin
-            @testset "$nelement $ngrid" for bc ∈ ("constant", "zero"),
+            @testset "$nelement $ngrid" for bc ∈ ("constant", "zero"), element_spacing_option ∈ ("uniform", "sqrt"),
                     nelement ∈ (1:5), ngrid ∈ (3:33)
 
                 # define inputs needed for the test
@@ -742,8 +741,7 @@ function runtests()
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
 				comm = MPI.COMM_NULL # dummy value
-                element_spacing_option = "uniform"
-				input = grid_input("coord", ngrid, nelement,
+                input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
                     "chebyshev_pseudospectral", fd_option, bc, adv_input, comm,
                     element_spacing_option)

From e24ef9c789642c3a0f5c9f435dbed99440f30f72 Mon Sep 17 00:00:00 2001
From: Michael Hardman <michael.hardman@physics.ox.ac.uk>
Date: Mon, 25 Sep 2023 10:15:40 +0100
Subject: [PATCH 09/12] Include sqrt grid in interpolation tests as suggested
 by @johnomotani.

---
 test/interpolation_tests.jl | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl
index 33cd0fbaa..8885224cc 100644
--- a/test/interpolation_tests.jl
+++ b/test/interpolation_tests.jl
@@ -27,8 +27,9 @@ 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, rtol) ∈
-                    (("finite_difference", 1.e-5), ("chebyshev_pseudospectral", 1.e-8)),
+                (discretization, element_spacing_option, rtol) ∈
+                    (("finite_difference", "uniform", 1.e-5), ("chebyshev_pseudospectral", "uniform", 1.e-8),
+                    ("chebyshev_pseudospectral", "sqrt", 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
@@ -36,7 +37,7 @@ function runtests()
 			nrank_per_block = 0 # dummy value
 			irank = 0 # dummy value
 			comm = MPI.COMM_NULL # dummy value
-            element_spacing_option = "uniform"
+            #element_spacing_option = "uniform"
 			input = grid_input("coord", ngrid, nelement,
 				nelement_local, nrank_per_block, irank, L,
 				discretization, fd_option, bc, adv_input, comm,

From 05b7bfecc2acac68871fe62e8f7dd0038b23aae4 Mon Sep 17 00:00:00 2001
From: Michael Hardman <michael.hardman@physics.ox.ac.uk>
Date: Mon, 25 Sep 2023 10:17:02 +0100
Subject: [PATCH 10/12] Include sqrt grid in wall boundary condition tests as
 suggested by @johnomotani. Note that the minimum number of elements for a
 nontrivial test are 5 and 6 -- the "sqrt" element spacing option has
 different behaviour for odd and even number of elements, so both should be
 tested.

---
 test/wall_bc_tests.jl | 60 +++++++++++++++++++++++++++++++++++++++----
 1 file changed, 55 insertions(+), 5 deletions(-)

diff --git a/test/wall_bc_tests.jl b/test/wall_bc_tests.jl
index 9b51063e2..2bcda81b0 100644
--- a/test/wall_bc_tests.jl
+++ b/test/wall_bc_tests.jl
@@ -82,6 +82,7 @@ test_input_finite_difference = Dict("n_ion_species" => 1,
                                     "z_nelement" => 1,
                                     "z_bc" => "wall",
                                     "z_discretization" => "finite_difference",
+                                    "z_element_spacing_option" => "uniform",
                                     "vpa_ngrid" => 400,
                                     "vpa_nelement" => 1,
                                     "vpa_L" => 8.0,
@@ -98,6 +99,32 @@ test_input_chebyshev = merge(test_input_finite_difference,
                                   "z_discretization" => "chebyshev_pseudospectral",
                                   "z_ngrid" => 9,
                                   "z_nelement" => 2,
+                                  "z_element_spacing_option" => "uniform",
+                                  "vpa_discretization" => "chebyshev_pseudospectral",
+                                  "vpa_ngrid" => 17,
+                                  "vpa_nelement" => 10,
+                                  "vz_discretization" => "chebyshev_pseudospectral",
+                                  "vz_ngrid" => 17,
+                                  "vz_nelement" => 10))
+                                  
+test_input_chebyshev_sqrt_grid_odd = merge(test_input_finite_difference,
+                             Dict("run_name" => "chebyshev_pseudospectral",
+                                  "z_discretization" => "chebyshev_pseudospectral",
+                                  "z_ngrid" => 9,
+                                  "z_nelement" => 5, # minimum nontrival nelement (odd)
+                                  "z_element_spacing_option" => "sqrt",
+                                  "vpa_discretization" => "chebyshev_pseudospectral",
+                                  "vpa_ngrid" => 17,
+                                  "vpa_nelement" => 10,
+                                  "vz_discretization" => "chebyshev_pseudospectral",
+                                  "vz_ngrid" => 17,
+                                  "vz_nelement" => 10))
+test_input_chebyshev_sqrt_grid_even = merge(test_input_finite_difference,
+                             Dict("run_name" => "chebyshev_pseudospectral",
+                                  "z_discretization" => "chebyshev_pseudospectral",
+                                  "z_ngrid" => 9,
+                                  "z_nelement" => 6, # minimum nontrival nelement (even)
+                                  "z_element_spacing_option" => "sqrt",
                                   "vpa_discretization" => "chebyshev_pseudospectral",
                                   "vpa_ngrid" => 17,
                                   "vpa_nelement" => 10,
@@ -125,7 +152,7 @@ function run_test(test_input, expected_phi, tolerance; args...)
     # update the default inputs
 
     # Convert keyword arguments to a unique name
-    name = test_input["run_name"]
+    name = test_input["run_name"] * ", with element spacing: " * test_input["z_element_spacing_option"]
     if length(args) > 0
         name = string(name, "_", (string(k, "-", v, "_") for (k, v) in args)...)
 
@@ -196,12 +223,15 @@ function run_test(test_input, expected_phi, tolerance; args...)
         input = grid_input("coord", test_input["z_ngrid"], test_input["z_nelement"], 
 						   test_input["z_nelement"], nrank_per_block, irank, 1.0,
                            test_input["z_discretization"], "", test_input["z_bc"],
-                           adv_input, comm, element_spacing_option)
+                           adv_input, comm, test_input["z_element_spacing_option"])
         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)
-        @test isapprox(phi_interp, cross_compare_phi, rtol=tolerance, atol=1.e-15)
+        if test_input["z_element_spacing_option"] == "uniform" 
+            # only support this test for uniform element spacing
+            phi_interp = interpolate_to_grid_z(cross_compare_points, phi[:, end], z, z_spectral)
+            @test isapprox(phi_interp, cross_compare_phi, rtol=tolerance, atol=1.e-15)
+        end
     end
 end
 
@@ -214,12 +244,32 @@ function runtests()
             run_test(test_input_finite_difference, nothing, 2.e-3)
         end
 
-        @testset "Chebyshev" begin
+        @testset "Chebyshev uniform" begin
             run_test(test_input_chebyshev,
                      [-1.1689445031600718, -0.7479504438063098, -0.6947559936893813,
                       -0.6917252442591313, -0.7180152498764835, -0.9980114095597415],
                      2.e-3)
         end
+        
+        @testset "Chebyshev sqrt grid odd" begin
+            run_test(test_input_chebyshev_sqrt_grid_odd,
+                     [-1.2047298885671576, -0.9431378294506091, -0.8084332392927167,
+                     -0.7812620422650213, -0.7233303514000929, -0.7003878610612269,
+                     -0.69572751349158, -0.6933148921301019, -0.6992503992521327,
+                     -0.7115787972775218, -0.7596015032228407, -0.795776514029509,
+                     -0.876303297135126, -1.1471244425913258],
+                     2.e-3)
+        end
+        @testset "Chebyshev sqrt grid even" begin
+            run_test(test_input_chebyshev_sqrt_grid_even,
+                     [-1.213617049279473, -1.0054529928344382, -0.871444761913497,
+                     -0.836017699317097, -0.7552110924643832, -0.7264644073096705,
+                     -0.7149147366621806, -0.6950077192395091, -0.6923364889119271,
+                     -0.6950077192395089, -0.7149147366621814, -0.7264644073096692,
+                     -0.7552110924643836, -0.8360176993170979, -0.8714447619134948,
+                     -1.0054529928344376, -1.2136170492794727],
+                     2.e-3)
+        end
     end
 end
 

From 178178a5f6351968665974db5c4f7f0bdb144c52 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Mon, 25 Sep 2023 11:48:20 +0100
Subject: [PATCH 11/12] Comment why sqrt spacing skips cross-compare test

[skip ci]
---
 test/calculus_tests.jl | 1 -
 test/wall_bc_tests.jl  | 4 +++-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/test/calculus_tests.jl b/test/calculus_tests.jl
index bdcd10453..ae198003c 100644
--- a/test/calculus_tests.jl
+++ b/test/calculus_tests.jl
@@ -38,7 +38,6 @@ function runtests()
 				nrank_per_block = 0 # dummy value
 				irank = 0 # dummy value
 				comm = MPI.COMM_NULL # dummy value 
-                #element_spacing_option = "uniform" # dummy value
 				input = grid_input("coord", ngrid, nelement,
                     nelement_local, nrank_per_block, irank, L,
                     discretization, fd_option, bc, adv_input, comm,
diff --git a/test/wall_bc_tests.jl b/test/wall_bc_tests.jl
index 2bcda81b0..c0b563393 100644
--- a/test/wall_bc_tests.jl
+++ b/test/wall_bc_tests.jl
@@ -228,7 +228,9 @@ function run_test(test_input, expected_phi, tolerance; args...)
 
         # Cross comparison of all discretizations to same benchmark
         if test_input["z_element_spacing_option"] == "uniform" 
-            # only support this test for uniform element spacing
+            # Only support this test for uniform element spacing.
+            # phi is better resolved by "sqrt" spacing grid, so disagrees with benchmark data from
+            # simulation with uniform element spacing.
             phi_interp = interpolate_to_grid_z(cross_compare_points, phi[:, end], z, z_spectral)
             @test isapprox(phi_interp, cross_compare_phi, rtol=tolerance, atol=1.e-15)
         end

From 3018b2b4331ffca190f28d11ec1aedd1c6dca202 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Mon, 25 Sep 2023 11:49:58 +0100
Subject: [PATCH 12/12] Increase CI timeout for 'test' workflow

Extra tests added mean the CI was timing out frequently.
---
 .github/workflows/test.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index a55d73248..de2c2690e 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -10,7 +10,7 @@ jobs:
       matrix:
         os: [ubuntu-latest, macOS-latest]
       fail-fast: false
-    timeout-minutes: 30
+    timeout-minutes: 40
 
     steps:
       - uses: actions/checkout@v2