Update paper.md

paper.md

@@ -158,7 +158,7 @@ Sandhu supported development of the package by building and validating code test
 the development of the package structure and organization.  This work
 was supported by the University of Nevada, Reno, Office of Undergraduate Research's
 Pack Research Experience Program (PREP) which supported Sukhreen Sandhu as a research assistant.
-This work benefitted from the DAPPER library which was referenced at times for the development
+This work benefited from the DAPPER library which was referenced at times for the development
 of DA schemes.
 
 # References

src/models/L96.jl

@@ -42,14 +42,14 @@ function dx_dt(x::T, t::Float64, dx_params::ParamDict) where {T <: VecA}
     F = dx_params["F"][1]::Float64
     x_dim = length(x)
     dx = Vector{Float64}(undef, x_dim)
-	
+
     for j in 1:x_dim
         # index j minus 2, modulo the system dimension
         j_m_2 = mod_indx!(j - 2, x_dim)
-        
+
         # index j minus 1, modulo the system dimension
         j_m_1 = mod_indx!(j - 1, x_dim)
-        
+
         # index j plus 1, modulo the system dimension
         j_p_1 = mod_indx!(j + 1, x_dim)
 
@@ -85,7 +85,7 @@ end
 
 function jacobian(x::Vector{Float64}, t::Float64, dx_params::ParamDict)
     """"This computes the Jacobian of the Lorenz 96, for arbitrary dimension, equation about the state x.
-    
+
     Note that this has been designed to load the entries in a standard zeros array but return a sparse array after
     to make a compromise between memory and computational resources."""
 
@@ -97,16 +97,16 @@ function jacobian(x::Vector{Float64}, t::Float64, dx_params::ParamDict)
 
         # index j minus 2, modulo the system dimension
         j_m_2 = mod_indx!(j - 2, x_dim)
-        
+
         # index j minus 1, modulo the system dimension
         j_m_1 = mod_indx!(j - 1, x_dim)
-        
+
         # index j plus 1, modulo the system dimension
         j_p_1 = mod_indx!(j + 1, x_dim)
-        
+
         # index j plus 2, modulo the system dimension
         j_p_2 = mod_indx!(j + 2, x_dim)
-        
+
         # load the jacobian entries in corresponding rows
         dxF[j_p_2, j] = -x[j_p_1]
         dxF[j_p_1, j] = x[j_p_2] - x[j_m_1]
@@ -117,7 +117,7 @@ function jacobian(x::Vector{Float64}, t::Float64, dx_params::ParamDict)
 end
 
 ########################################################################################################################
-# Auxiliary functions for the 2nd order Taylor-Stratonovich expansion below. These need to be computed once, only as 
+# Auxiliary functions for the 2nd order Taylor-Stratonovich expansion below. These need to be computed once, only as
 # a function of the order of truncation of the fourier series, p for full timeseries.
 
 function ρ(p::Int64)
@@ -131,21 +131,21 @@ end
 ########################################################################################################################
 # 2nd order strong taylor SDE step
 # This method is derived in
-# Grudzien, C. et al.: On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark 
+# Grudzien, C. et al.: On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark
 # twin experiments, Geosci. Model Dev., 13, 1903–1924, https://doi.org/10.5194/gmd-13-1903-2020, 2020.
 # This depends on rho and alpha as above
 # NOTE: this still needs to be debugged
 
 function l96s_tay2_step!(x::Vector{Float64}, t::Float64, kwargs::Dict{String,Any})
     """One step of integration rule for l96 second order taylor rule
 
-    The rho and alpha are to be computed by the auxiliary functions, depending only on p, and supplied for all steps.  
-    This is the general formulation which includes, eg. dependence on the truncation of terms in the auxilliary 
+    The rho and alpha are to be computed by the auxiliary functions, depending only on p, and supplied for all steps.
+    This is the general formulation which includes, eg. dependence on the truncation of terms in the auxilliary
     function C with respect to the parameter p.  In general, truncation at p=1 is all that is necessary for order
-    2.0 convergence, and in this case C below is identically equal to zero.  This auxilliary function can be removed 
+    2.0 convergence, and in this case C below is identically equal to zero.  This auxilliary function can be removed
     (and is removed) in other implementations for simplicity."""
 
-    # Infer model and parameters 
+    # Infer model and parameters
     sys_dim = length(x)
     dx_params = kwargs["dx_params"]::ParamDict
     h = kwargs["h"]::Float64
@@ -159,29 +159,29 @@ function l96s_tay2_step!(x::Vector{Float64}, t::Float64, kwargs::Dict{String,Any
     Jac_x = jacobian(x, 0.0, dx_params)
 
     ## random variables
-    # Vectors ξ, μ, ϕ are sys_dim X 1 vectors of iid standard normal variables, 
+    # Vectors ξ, μ, ϕ are sys_dim X 1 vectors of iid standard normal variables,
     # ζ and η are sys_dim X p matrices of iid standard normal variables. Functional relationships describe each
     # variable W_j as the transformation of ξ_j to be of variace given by the length of the time step h. The functions
     # of random Fourier coefficients a_i, b_i are given in terms μ / η and ϕ / ζ respectively.
-    
+
     # draw standard normal samples
     rndm = rand(Normal(), sys_dim, 2*p + 3)
     ξ = rndm[:, 1]
-    
+
     μ = rndm[:, 2]
     ϕ = rndm[:, 3]
 
     ζ = rndm[:, 4: p+3]
     η = rndm[:, p+4: end]
-    
+
     ### define the auxiliary functions of random fourier coefficients, a and b
-    
+
     # denominators for the a series
     denoms = repeat((1.0 ./ Vector{Float64}(1:p)), 1, sys_dim)
 
     # vector of sums defining a terms
     a = -2.0 * sqrt(h * ρ) * μ - sqrt(2.0*h) * sum(ζ' .* denoms, dims=1)' / π
-    
+
     # denominators for the b series
     denoms = repeat((1.0 ./ Vector{Float64}(1:p).^2.0), 1, sys_dim)
 
@@ -190,12 +190,12 @@ function l96s_tay2_step!(x::Vector{Float64}, t::Float64, kwargs::Dict{String,Any
 
     # vector of first order Stratonovich integrals
     J_pdelta = (h / 2.0) * (sqrt(h) * ξ + a)
-    
+
     ### auxiliary functions for higher order stratonovich integrals ###
     # C function is optional for higher precision but does not change order of convergence to leave out
     function C(l, j)
         if p == 1
-            return 0.0 
+            return 0.0
         end
         c = zeros(p, p)
         # define the coefficient as a sum of matrix entries where r and k do not agree
@@ -230,7 +230,7 @@ function l96s_tay2_step!(x::Vector{Float64}, t::Float64, kwargs::Dict{String,Any
 
     # the final vectorized step forward is given as
     x .= collect(Iterators.flatten(
-            x + dx * h + h^2.0 * 0.5 * Jac_x * dx +  # deterministic taylor step 
+            x + dx * h + h^2.0 * 0.5 * Jac_x * dx +  # deterministic taylor step
             diffusion * sqrt(h) * ξ +                # stochastic euler step
             diffusion * Jac_x * J_pdelta +           # stochastic first order taylor step
             diffusion^2.0 * (Ψ_plus - Ψ_minus)       # stochastic second order taylor step
@@ -240,4 +240,3 @@ end
 ########################################################################################################################
 
 end
-

test/Test.jl

@@ -0,0 +1,3 @@
+module Test
+
+end

test/TestDeSolvers.jl

@@ -7,15 +7,18 @@ using DataAssimilationBenchmarks.DeSolvers
 using DataAssimilationBenchmarks.L96
 using LsqFit
 using Test
-
 ########################################################################################################################
 ########################################################################################################################
 # Define test function for analytical verification of discretization error
+VecA = Union{Vector{Float64}, SubArray{Float64, 1}}
+ParamDict = Union{Dict{String, Array{Float64}}, Dict{String, Vector{Float64}}}
 
-function exponentialODE(x, t, dx_params)
+function exponentialODE(x::T, t::Float64, dx_params::ParamDict) where {T <: VecA}
     # fill in exponential function here in the form of the other time derivatives
     # like L96.dx_dt, where the function is in terms of dx/dt=e^(t) so that we can
     # compute the answer analytically for comparision
+
+    return [exp(t)]
 end
 
 
@@ -25,13 +28,13 @@ end
 function expDiscretizationError(step_model!, h)
     # continuous time length of the integration
     tanl = 0.1
-    
+
     # discrete integration steps
     fore_steps = convert(Int64, tanl/h)
     time_steps = LinRange(0, tanl, fore_steps + 1)
 
     # initial data for the exponential function
-    x = 1.0
+    x = [1.0]
 
     # set the kwargs for the integration scheme
     # with empty values for the uneccessary parameters
@@ -43,12 +46,12 @@ function expDiscretizationError(step_model!, h)
             "dx_params" => dx_params,
             "dx_dt" => exponentialODE,
             )
-
     for i in 1:fore_steps
-        step_model!(x, time_steps[i], kwargs) 
+        step_model!(x, time_steps[i], kwargs)
     end
 
-    abs(x - exp(tanl))
+    abs(x[1] - exp(tanl))
+
 end
 
 
@@ -57,15 +60,24 @@ end
 
 function calculateOrderConvergence(step_model!)
     # set step sizes in increasing order for log-10 log-10 analysis
-    h_range = .1 .^[3, 2, 1]
+    h_range = [0.0001, 0.001, 0.01]
     error_range = Vector{Float64}(undef, 3)
-    
+
     # loop the discretization and calculate the errors
     for i in 1:length(h_range)
         error_range[i] = expDiscretizationError(step_model!, h_range[i])
     end
 
-    # Set the least squares estimate for the line in log-10 log-10 scale using LsqFit 
+    h_range_log10 = log10.(h_range)
+    error_range_log10 = log10.(error_range)
+
+    function model_lsq_squares(x,p)
+        @.p[1] + p[2]*x
+    end
+
+    fit = curve_fit(model_lsq_squares, h_range_log10, error_range_log10, [1.0,1.0])
+    return coef(fit)
+    # Set the least squares estimate for the line in log-10 log-10 scale using LsqFit
     # where we have the h_range in the x-axis and the error_range in the
     # y-axis.  The order of convergence is defined by the slope of this line
     # whereas the intercept is a constant factor (not necessary to report here).
@@ -74,16 +86,17 @@ function calculateOrderConvergence(step_model!)
 end
 
 
+
 ########################################################################################################################
 # Wrapper function to be supplied to runtests
 
 function testEMExponential()
-    slope = calculateOrderConvergence(em_step!)
-    
-    if abs(slope - 1.0) > 0.1
-        false
+    coef = calculateOrderConvergence(em_step!)
+
+    if abs(coef[2] - 1.0) > 0.1
+        return false
     else
-        true
+        return true
     end
 end
 
@@ -92,16 +105,18 @@ end
 # Wrapper function to be supplied to runtests
 
 function testRKExponential()
-    slope = calculateOrderConvergence(rk_step!)
-    
-    if abs(slope - 4.0) > 0.1
-        false
+    coef = calculateOrderConvergence(rk4_step!)
+
+    if abs(coef[2] - 4.0) > 0.1
+        return false
     else
-        true
+        return true
     end
 end
 
 
+testRKExponential()
+
 ########################################################################################################################
 
 end