From 6127e95840e36556747dd3733dbd1addb9dd75db Mon Sep 17 00:00:00 2001 From: Sukhreen Sandhu <70529624+sukhreens@users.noreply.github.com> Date: Fri, 12 Nov 2021 21:49:08 -0800 Subject: [PATCH 1/2] Update paper.md 1 edit --- paper.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/paper.md b/paper.md index 8a1ca78a..b480ca54 100644 --- a/paper.md +++ b/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 From cafd53c6c4a58bc4c7121fae6b04e1e4edea27b5 Mon Sep 17 00:00:00 2001 From: Sukhreen Sandhu Date: Fri, 19 Nov 2021 14:38:04 -0800 Subject: [PATCH 2/2] changes to TestDeSolver --- src/models/L96.jl | 47 +++++++++++++-------------- test/Test.jl | 3 ++ test/TestDeSolvers.jl | 55 ++++++++++++++++++++------------ test/TestTimeSeriesGeneration.jl | 40 +++++++++++------------ test/runtests.jl | 9 +++++- 5 files changed, 88 insertions(+), 66 deletions(-) create mode 100644 test/Test.jl diff --git a/src/models/L96.jl b/src/models/L96.jl index 1eaebd36..48f9b0c5 100644 --- a/src/models/L96.jl +++ b/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,7 +131,7 @@ 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 @@ -139,13 +139,13 @@ end 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 - diff --git a/test/Test.jl b/test/Test.jl new file mode 100644 index 00000000..5a384d8d --- /dev/null +++ b/test/Test.jl @@ -0,0 +1,3 @@ +module Test + +end diff --git a/test/TestDeSolvers.jl b/test/TestDeSolvers.jl index f705ef38..cc4f4a67 100644 --- a/test/TestDeSolvers.jl +++ b/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 diff --git a/test/TestTimeSeriesGeneration.jl b/test/TestTimeSeriesGeneration.jl index 03eddcb7..6cdbff88 100644 --- a/test/TestTimeSeriesGeneration.jl +++ b/test/TestTimeSeriesGeneration.jl @@ -1,4 +1,4 @@ -####################################################################################################################### +#git ####################################################################################################################### module TestTimeSeriesGeneration ####################################################################################################################### # imports and exports @@ -8,13 +8,12 @@ using JLD using Random ####################################################################################################################### # Test generation of the L96 model time series - function testL96() - # define the model and the solver +# define the model and the solver dx_dt = L96.dx_dt step_model! = DeSolvers.em_step! - # set model and experimental parameters +# set model and experimental parameters F = 8.0 h = 0.001 nanl = 1000 @@ -24,11 +23,11 @@ function testL96() seed = 0 Random.seed!(seed) - # define the dx_params dict +# define the dx_params dict dx_params = Dict{String, Array{Float64}}("F" => [F]) fore_steps = convert(Int64, tanl/h) - # set the kwargs for the integration scheme +# set the kwargs for the integration scheme kwargs = Dict{String, Any}( "h" => h, "diffusion" => diffusion, @@ -36,13 +35,13 @@ function testL96() "dx_dt" => L96.dx_dt, ) - # set arbitrary initial condition +# set arbitrary initial condition xt = ones(sys_dim) - # pre-allocate storage for the time series observations +# pre-allocate storage for the time series observations tobs = Array{Float64}(undef,sys_dim, nanl) - # loop the experiment, taking observations at time length tanl +# loop the experiment, taking observations at time length tanl for i in 1:nanl for j in 1:fore_steps step_model!(xt, 0.0, kwargs) @@ -50,8 +49,8 @@ function testL96() tobs[:,i] = xt end - # define the file name for the experiment output - # dynamically based on experiment parameters +# define the file name for the experiment output +# dynamically based on experiment parameters fname = "time_series_data_seed_" * lpad(seed, 4, "0") * "_dim_" * lpad(sys_dim, 2, "0") * "_diff_" * rpad(diffusion, 5, "0") * @@ -61,7 +60,7 @@ function testL96() "_h_" * rpad(h, 5, "0") * ".jld" - # define the experimental data in a dictionary to write with JLD +# define the experimental data in a dictionary to write with JLD data = Dict{String, Any}( "h" => h, "diffusion" => diffusion, @@ -73,33 +72,32 @@ function testL96() ) path = "../data/time_series/" - # test to see if the data can be written to standard output directory +# test to see if the data can be written to standard output directory function write_file() try save(path * fname, data) - true + did_write = true catch # if not, set test case false - false + did_write = false end end - # test to see if the data can be read from standard output directory +# test to see if the data can be read from standard output directory function load_file() try tmp = load(path * fname) - true + did_read = true catch # if not, set test case false - false + did_read = false end end - return write_file() && load_file() - +write_file() +load_file() end - ####################################################################################################################### end diff --git a/test/runtests.jl b/test/runtests.jl index 61f3c8f3..eb68d33b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,9 +8,10 @@ using DataAssimilationBenchmarks.L96 using Test ######################################################################################################################## -# include test sub-modules +# include test sub-modules include("TestL96.jl") include("TestTimeSeriesGeneration.jl") +include("TestDeSolvers.jl") ######################################################################################################################## # Run tests @@ -25,6 +26,12 @@ end @test TestTimeSeriesGeneration.testL96() end +# test set 3: +@testset "Calculate Order Convergence" begin + @test TestDeSolvers.testEMExponential() + @test TestDeSolvers.testRKExponential() +end + ########################################################################################################################