Merge pull request #128 from mabarnes/restart-interpolation

Allow interpolation between grids when restarting

.github/workflows/longtest.yml

@@ -14,7 +14,7 @@ jobs:
       matrix:
         os: [ubuntu-latest, macOS-latest]
       fail-fast: false
-    timeout-minutes: 75
+    timeout-minutes: 90
 
     steps:
       - uses: actions/checkout@v2

.github/workflows/parallel_test.yml

@@ -12,7 +12,7 @@ jobs:
         include:
           - julia_version: '1.8'
       fail-fast: false
-    timeout-minutes: 90
+    timeout-minutes: 120
 
     steps:
       - uses: actions/checkout@v2

.github/workflows/test.yml

@@ -10,7 +10,7 @@ jobs:
       matrix:
         os: [ubuntu-latest, macOS-latest]
       fail-fast: false
-    timeout-minutes: 40
+    timeout-minutes: 50
 
     steps:
       - uses: actions/checkout@v2

docs/src/constants.md

@@ -0,0 +1,6 @@
+`constants`
+===========
+
+```@autodocs
+Modules = [moment_kinetics.constants]
+```

docs/src/krook_collisions.md

@@ -0,0 +1,6 @@
+`krook_collisions`
+==================
+
+```@autodocs
+Modules = [moment_kinetics.krook_collisions]
+```

docs/src/reference_parameters.md

@@ -0,0 +1,6 @@
+`reference_parameters`
+======================
+
+```@autodocs
+Modules = [moment_kinetics.reference_parameters]
+```

src/calculus.jl

@@ -6,7 +6,8 @@ export derivative!, second_derivative!
 export reconcile_element_boundaries_MPI!
 export integral
 
-using ..moment_kinetics_structs: chebyshev_info
+using ..moment_kinetics_structs: discretization_info, null_spatial_dimension_info,
+                                 null_velocity_dimension_info
 using ..type_definitions: mk_float, mk_int
 using MPI
 using ..communication: block_rank
@@ -42,7 +43,7 @@ function elementwise_second_derivative! end
 
 Upwinding derivative.
 """
-function derivative!(df, f, coord, adv_fac, spectral::Union{Bool,<:chebyshev_info})
+function derivative!(df, f, coord, adv_fac, spectral::discretization_info)
     # get the derivative at each grid point within each element and store in
     # coord.scratch_2d
     elementwise_derivative!(coord, f, adv_fac, spectral)
@@ -65,16 +66,11 @@ function derivative!(df, f, coord, spectral)
     derivative_elements_to_full_grid!(df, coord.scratch_2d, coord)
 end
 
-function second_derivative!(df, f, coord, spectral::Bool)
-    # Finite difference version must use an appropriate second derivative stencil, not
-    # apply the 1st derivative twice as for the spectral element method
-
-    # get the derivative at each grid point within each element and store in
-    # coord.scratch_2d
-    elementwise_second_derivative!(coord, f, spectral)
-    # map the derivative from the elem;ntal grid to the full grid;
-    # at element boundaries, use the average of the derivatives from neighboring elements.
-    derivative_elements_to_full_grid!(df, coord.scratch_2d, coord)
+# Special versions for 'null' coordinates with only one point
+function derivative!(df, f, coord, spectral::Union{null_spatial_dimension_info,
+                                                   null_velocity_dimension_info})
+    df .= 0.0
+    return nothing
 end
 
 function second_derivative!(d2f, f, Q, coord, spectral)

src/chebyshev.jl

@@ -15,7 +15,27 @@ using ..array_allocation: allocate_float, allocate_complex
 using ..clenshaw_curtis: clenshawcurtisweights
 import ..calculus: elementwise_derivative!
 import ..interpolation: interpolate_to_grid_1d!
-using ..moment_kinetics_structs: chebyshev_info
+using ..moment_kinetics_structs: discretization_info
+
+"""
+Chebyshev pseudospectral discretization
+"""
+struct chebyshev_info{TForward <: FFTW.cFFTWPlan, TBackward <: AbstractFFTs.ScaledPlan} <: discretization_info
+    # fext is an array for storing f(z) on the extended domain needed
+    # to perform complex-to-complex FFT using the fact that f(theta) is even in theta
+    fext::Array{Complex{mk_float},1}
+    # Chebyshev spectral coefficients of distribution function f
+    # first dimension contains location within element
+    # second dimension indicates the element
+    f::Array{mk_float,2}
+    # Chebyshev spectral coefficients of derivative of f
+    df::Array{mk_float,1}
+    # plan for the complex-to-complex, in-place, forward Fourier transform on Chebyshev-Gauss-Lobatto grid
+    forward::TForward
+    # plan for the complex-to-complex, in-place, backward Fourier transform on Chebyshev-Gauss-Lobatto grid
+    #backward_transform::FFTW.cFFTWPlan
+    backward::TBackward
+end
 
 """
 create arrays needed for explicit Chebyshev pseudospectral treatment

src/coordinates.jl

@@ -10,8 +10,10 @@ using ..type_definitions: mk_float, mk_int
 using ..array_allocation: allocate_float, allocate_int
 using ..calculus: derivative!
 using ..chebyshev: scaled_chebyshev_grid, setup_chebyshev_pseudospectral
+using ..finite_differences: finite_difference_info
 using ..quadrature: composite_simpson_weights
 using ..input_structs: advection_input
+using ..moment_kinetics_structs: null_spatial_dimension_info, null_velocity_dimension_info
 
 using MPI
 
@@ -161,16 +163,23 @@ function define_coordinate(input, parallel_io::Bool=false)
         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)
 
-    if input.discretization == "chebyshev_pseudospectral" && coord.n > 1
+    if coord.n == 1 && occursin("v", coord.name)
+        spectral = null_velocity_dimension_info()
+        coord.duniform_dgrid .= 1.0
+    elseif coord.n == 1
+        spectral = null_spatial_dimension_info()
+        coord.duniform_dgrid .= 1.0
+    elseif input.discretization == "chebyshev_pseudospectral"
         # 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
+        # finite_difference_info is just a type so that derivative methods, etc., dispatch
+        # to the finite difference versions, it does not contain any information.
+        spectral = finite_difference_info()
         coord.duniform_dgrid .= 1.0
     end
 

src/finite_differences.jl

@@ -3,7 +3,14 @@
 module finite_differences
 
 using ..type_definitions: mk_float
-import ..calculus: elementwise_derivative!, elementwise_second_derivative!
+import ..calculus: elementwise_derivative!, elementwise_second_derivative!,
+                   second_derivative!
+using ..moment_kinetics_structs: discretization_info
+
+"""
+Finite difference discretization
+"""
+struct finite_difference_info <: discretization_info end
 
 """
 """
@@ -26,38 +33,55 @@ function fd_check_option(option, ngrid)
 end
 
 """
-    elementwise_derivative!(coord, f, adv_fac, not_spectral::Bool)
+    elementwise_derivative!(coord, f, adv_fac, not_spectral::finite_difference_info)
 
 Calculate the derivative of f using finite differences, with particular scheme
 specified by coord.fd_option; result stored in coord.scratch_2d.
 """
-function elementwise_derivative!(coord, f, adv_fac, not_spectral::Bool)
+function elementwise_derivative!(coord, f, adv_fac, not_spectral::finite_difference_info)
     return derivative_finite_difference!(coord.scratch_2d, f, coord.cell_width, adv_fac,
         coord.bc, coord.fd_option, coord.igrid, coord.ielement)
 end
 
 """
-    elementwise_derivative!(coord, f, not_spectral::Bool)
+    elementwise_derivative!(coord, f, not_spectral::finite_difference_info)
 
 Calculate the derivative of f using 4th order centered finite differences; result stored
 in coord.scratch_2d.
 """
-function elementwise_derivative!(coord, f, not_spectral::Bool)
+function elementwise_derivative!(coord, f, not_spectral::finite_difference_info)
     return derivative_finite_difference!(coord.scratch_2d, f, coord.cell_width,
         coord.bc, "fourth_order_centered", coord.igrid, coord.ielement)
 end
 
 """
-    elementwise_second_derivative!(coord, f, not_spectral::Bool)
+    elementwise_second_derivative!(coord, f, not_spectral::finite_difference_info)
 
 Calculate the second derivative of f using 2nd order centered finite differences; result
 stored in coord.scratch_2d.
 """
-function elementwise_second_derivative!(coord, f, not_spectral::Bool)
+function elementwise_second_derivative!(coord, f, not_spectral::finite_difference_info)
     return second_derivative_finite_difference!(coord.scratch_2d, f, coord.cell_width,
         coord.bc, coord.igrid, coord.ielement)
 end
 
+function second_derivative!(df, f, Q, coord, spectral::finite_difference_info)
+    # Finite difference version must use an appropriate second derivative stencil, not
+    # apply the 1st derivative twice as for the spectral element method
+
+    if !all(Q .== 1.0)
+        error("Finite difference implementation of second derivative does not support "
+              * "Q!=1.")
+    end
+
+    # get the derivative at each grid point within each element and store in
+    # coord.scratch_2d
+    elementwise_second_derivative!(coord, f, spectral)
+    # map the derivative from the elem;ntal grid to the full grid;
+    # at element boundaries, use the average of the derivatives from neighboring elements.
+    derivative_elements_to_full_grid!(df, coord.scratch_2d, coord)
+end
+
 """
 """
 function derivative_finite_difference!(df, f, del, adv_fac, bc, fd_option, igrid, ielement)

src/hermite_spline_interpolation.jl

@@ -1,6 +1,7 @@
 module hermite_spline_interpolation
 
 using ..calculus: derivative!
+using ..finite_differences: finite_difference_info
 import ..interpolation: interpolate_to_grid_1d!
 
 """
@@ -17,11 +18,12 @@ f : Array{mk_float}
     Field to be interpolated
 coord : coordinate
     `coordinate` struct giving the coordinate along which f varies
-not_spectral : Bool
-    A Bool argument here indicates that the coordinate is not spectral-element
-    discretized, i.e. it is on a uniform ('finite difference') grid.
+not_spectral : finite_difference_info
+    A finite_difference_info argument here indicates that the coordinate is not
+    spectral-element discretized, i.e. it is on a uniform ('finite difference') grid.
 """
-function interpolate_to_grid_1d!(result, new_grid, f, coord, not_spectral::Bool)
+function interpolate_to_grid_1d!(result, new_grid, f, coord,
+                                 not_spectral::finite_difference_info)
     x = coord.grid
     n_new = length(new_grid)
 

src/interpolation.jl

@@ -8,6 +8,7 @@ module interpolation
 export interpolate_to_grid_z
 
 using ..array_allocation: allocate_float
+using ..moment_kinetics_structs: null_spatial_dimension_info, null_velocity_dimension_info
 using ..type_definitions: mk_float
 
 """
@@ -23,11 +24,33 @@ f : Array{mk_float}
     Field to be interpolated
 coord : coordinate
     `coordinate` struct giving the coordinate along which f varies
-spectral : Bool or chebyshev_info
+spectral : discretization_info
     struct containing information for discretization, whose type determines which method
     is used.
 """
-function interpolate_to_grid_1d!() end
+function interpolate_to_grid_1d! end
+
+function interpolate_to_grid_1d!(result, new_grid, f, coord,
+                                 spectral::null_spatial_dimension_info)
+    # There is only one point in the 'old grid' represented by coord (as indicated by the
+    # type of the `spectral` argument), and we are interpolating in a spatial dimension.
+    # Assume that the variable should be taken to be constant in this dimension to
+    # 'interpolate'.
+    result .= f[1]
+
+    return nothing
+end
+
+function interpolate_to_grid_1d!(result, new_grid, f, coord,
+                                 spectral::null_velocity_dimension_info)
+    # There is only one point in the 'old grid' represented by coord (as indicated by the
+    # type of the `spectral` argument), and we are interpolating in a velocity space
+    # dimension. Assume that the profile 'should be' a Maxwellian over the new grid, with
+    # a width of 1 in units of the reference speed.
+    @. result = f[1] * exp(-new_grid^2)
+
+    return nothing
+end
 
 """
 Interpolation from a regular grid to a 1d grid with arbitrary spacing