Added derivative vector D0 for computing the derivative at vperp = 0 …

…(eventually) for imposing the regularity condition d F / d vperp. In normalised coordinates, D0 computes the derivative at x = -1, which is useful on the Radau grid because the grid runs in normalised units from (-1,1]. Note that like the other derivative matrices and vectors, an explicit scale length must be included. A test script is provided to check the functionality.

src/chebyshev.jl

@@ -39,6 +39,8 @@ struct chebyshev_base_info{TForward <: FFTW.cFFTWPlan, TBackward <: AbstractFFTs
     backward::TBackward
     # elementwise differentiation matrix (ngrid*ngrid)
     Dmat::Array{mk_float,2}
+    # elementwise differentiation vector (ngrid) for the point x = -1
+    D0::Array{mk_float,1}
 end
 
 struct chebyshev_info{TForward <: FFTW.cFFTWPlan, TBackward <: AbstractFFTs.ScaledPlan} <: discretization_info
@@ -72,9 +74,11 @@ function setup_chebyshev_pseudospectral_lobatto(coord)
     # create array for differentiation matrix 
     Dmat = allocate_float(coord.ngrid, coord.ngrid)
     cheb_derivative_matrix_elementwise!(Dmat,coord.ngrid)
+    D0 = allocate_float(coord.ngrid)
+    D0 .= Dmat[1,:]
     # return a structure containing the information needed to carry out
     # a 1D Chebyshev transform
-    return chebyshev_base_info(fext, fcheby, dcheby, forward_transform, backward_transform, Dmat)
+    return chebyshev_base_info(fext, fcheby, dcheby, forward_transform, backward_transform, Dmat, D0)
 end
 
 function setup_chebyshev_pseudospectral_radau(coord)
@@ -93,9 +97,11 @@ function setup_chebyshev_pseudospectral_radau(coord)
         # create array for differentiation matrix 
         Dmat = allocate_float(coord.ngrid, coord.ngrid)
         cheb_derivative_matrix_elementwise_radau_by_FFT!(Dmat, coord, fcheby, dcheby, fext, forward_transform)
+        D0 = allocate_float(coord.ngrid)
+        cheb_lower_endpoint_derivative_vector_elementwise_radau_by_FFT!(D0, coord, fcheby, dcheby, fext, forward_transform)
         # return a structure containing the information needed to carry out
         # a 1D Chebyshev transform
-        return chebyshev_base_info(fext, fcheby, dcheby, forward_transform, backward_transform, Dmat)
+        return chebyshev_base_info(fext, fcheby, dcheby, forward_transform, backward_transform, Dmat, D0)
 end
 
 """
@@ -780,6 +786,22 @@ function chebyshev_radau_forward_transform!(chebyf, fext, ff, transform, n)
         # inverse Chebyshev transform to get df/dcoord
         chebyshev_radau_backward_transform!(df, cheby_fext, cheby_df, forward, coord.ngrid)
     end
+    function chebyshev_radau_derivative_lower_endpoint(ff, cheby_f, cheby_df, cheby_fext, forward, coord)
+        # calculate the Chebyshev coefficients of the real-space function ff and return
+        # as cheby_f
+        chebyshev_radau_forward_transform!(cheby_f, cheby_fext, ff, forward, coord.ngrid)
+        # calculate the Chebyshev coefficients of the derivative of ff with respect to coord.grid
+        chebyshev_spectral_derivative!(cheby_df, cheby_f)
+        # form the derivative at x = - 1 using that T_n(-1) = (-1)^n
+        # and converting the normalisation factors to undo the normalisation in the FFT
+        # df = d0 + sum_n=1 (-1)^n d_n/2 with d_n the coeffs
+        # of the Cheb derivative in the Fourier representation
+        df = cheby_df[1]
+        for i in 2:coord.ngrid
+            df += ((-1)^(i-1))*0.5*cheby_df[i]
+        end
+        return df
+    end
 
 
 """
@@ -888,5 +910,22 @@ https://people.maths.ox.ac.uk/trefethen/pdetext.html
             D[j,j] = -sum(D[j,:])
         end
     end
+
+    function cheb_lower_endpoint_derivative_vector_elementwise_radau_by_FFT!(D::Array{Float64,1}, coord, f, df, fext, forward)
+        ff_buffer = Array{Float64,1}(undef,coord.ngrid)
+        df_buffer = Array{Float64,1}(undef,coord.ngrid)
+        # use response matrix approach to calculate derivative vector D 
+        for j in 1:coord.ngrid 
+            ff_buffer .= 0.0 
+            ff_buffer[j] = 1.0
+            @views df_buffer = chebyshev_radau_derivative_lower_endpoint(ff_buffer[:],
+                f[:,1], df, fext, forward, coord)
+            D[j] = df_buffer # assign appropriate value of derivative vector 
+        end
+        # correct diagonal elements to gurantee numerical stability
+        # gives D*[1.0, 1.0, ... 1.0] = [0.0, 0.0, ... 0.0]
+        D[1] = 0.0
+        D[1] = -sum(D[:])
+    end
 
 end

test_scripts/chebyshev_radau_test.jl

@@ -0,0 +1,107 @@
+export chebyshevradau_test
+
+using Printf
+using Plots
+using LaTeXStrings
+using MPI
+using Measures
+
+import moment_kinetics
+using moment_kinetics.chebyshev
+using moment_kinetics.input_structs: grid_input, advection_input
+using moment_kinetics.coordinates: define_coordinate
+using moment_kinetics.calculus: derivative!
+
+
+function print_matrix(matrix,name,n,m)
+    println("\n ",name," \n")
+    for i in 1:n
+        for j in 1:m
+            @printf("%.3f ", matrix[i,j])
+        end
+        println("")
+    end
+    println("\n")
+end
+
+function print_vector(vector,name,m)
+    println("\n ",name," \n")
+    for j in 1:m
+        @printf("%.3f ", vector[j])
+    end
+    println("")
+    println("\n")
+end 
+"""
+Test for derivative vector D0 that is used to compute 
+the numerical derivative on the Chebyshev-Radau elements
+at the lower endpoint of the domain (-1,1] in the normalised
+coordinate x. Here in the tests the shifted coordinate y 
+is used with the vperp label so that the grid runs from (0,L].
+"""
+function chebyshevradau_test(; ngrid=5, L_in=3.0)
+
+    # elemental grid tests 
+    #ngrid = 17
+    #nelement = 4
+    y_ngrid = ngrid #number of points per element 
+    y_nelement_local = 1 # number of elements per rank
+    y_nelement_global = y_nelement_local # total number of elements 
+    y_L = L_in
+    bc = "zero" 
+    discretization = "gausslegendre_pseudospectral"
+    # fd_option and adv_input not actually used so given values unimportant
+    fd_option = "fourth_order_centered"
+    cheb_option = "matrix"
+    adv_input = advection_input("default", 1.0, 0.0, 0.0)
+    nrank = 1
+    irank = 0#1
+    comm = MPI.COMM_NULL
+    element_spacing_option = "uniform"
+    # create the 'input' struct containing input info needed to create a
+    # coordinate
+    y_name = "vperp" # to use radau grid
+    y_input = grid_input(y_name, y_ngrid, y_nelement_global, y_nelement_local, 
+            nrank, irank, y_L, discretization, fd_option, cheb_option, bc, adv_input,comm,element_spacing_option)
+    y, y_spectral = define_coordinate(y_input)
+
+    Dmat = y_spectral.radau.Dmat
+    print_matrix(Dmat,"Radau Dmat",y.ngrid,y.ngrid)
+    D0 = y_spectral.radau.D0
+    print_vector(D0,"Radau D0",y.ngrid)
+
+    ff_err = Array{Float64,1}(undef,y.n)
+    ff = Array{Float64,1}(undef,y.n)
+    for iy in 1:y.n
+        ff[iy] = exp(-4.0*y.grid[iy])
+    end
+    df_exact = -4.0
+    df_num = sum(D0.*ff)/y.element_scale[1]
+    df_err = abs(df_num - df_exact)
+    println("f(y) = exp(-4 y) test")
+    println("exact df: ",df_exact," num df: ",df_num," abs(err): ",df_err) 
+
+    for iy in 1:y.n
+        ff[iy] = sin(y.grid[iy])
+    end
+    df_exact = 1.0
+    df_num = sum(D0.*ff)/y.element_scale[1]
+    df_err = abs(df_num - df_exact)
+    println("f(y) = sin(y) test")
+    println("exact df: ",df_exact," num df: ",df_num," abs(err): ",df_err) 
+    for iy in 1:y.n
+        ff[iy] = y.grid[iy] + (y.grid[iy])^2 + (y.grid[iy])^3
+    end
+    df_exact = 1.0
+    df_num = sum(D0.*ff)/y.element_scale[1]
+    df_err = abs(df_num - df_exact)
+    println("f(y) = y + y^2 + y^3 test")
+    println("exact df: ",df_exact," num df: ",df_num," abs(err): ",df_err) 
+end
+
+if abspath(PROGRAM_FILE) == @__FILE__
+    using Pkg
+    Pkg.activate(".")
+
+    chebyshevradau_test()
+end