Merge pull request #462 from isaacsas/add_vrj_test

add implicit solver VariableRateJump test

HISTORY.md

@@ -3,7 +3,8 @@
 ## JumpProcesses unreleased (master branch)
 
 ## 9.14
-  - Added the constant complexity next reaction method (CCNRM). 
+
+  - Added the constant complexity next reaction method (CCNRM).
 
 ## 9.13
 

Project.toml

@@ -49,6 +49,7 @@ julia = "1.10"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
@@ -57,4 +58,5 @@ StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["DiffEqCallbacks", "LinearAlgebra", "OrdinaryDiffEq", "SafeTestsets", "StableRNGs", "Statistics", "StochasticDiffEq", "Test", "FastBroadcast"]
+test = ["DiffEqCallbacks", "LinearAlgebra", "LinearSolve", "OrdinaryDiffEq", 
+    "SafeTestsets", "StableRNGs", "Statistics", "StochasticDiffEq", "Test", "FastBroadcast"]

test/table_test.jl

@@ -65,7 +65,7 @@ end
 ##### PRIORITY TIME TABLE TESTS FOR CCNRM
 mintime = 0.0;
 maxtime = 100.0;
-timestep = 1.5/16;
+timestep = 1.5 / 16;
 times = [2.0, 8.0, 13.0, 15.0, 74.0]
 
 ptt = DJ.PriorityTimeTable(times, mintime, timestep)
@@ -89,7 +89,7 @@ end
 @test DJ.getfirst(ptt) === (0, 0) # No more left.
 
 mintime = 66.0;
-timestep = 0.75/16;
+timestep = 0.75 / 16;
 DJ.rebuild!(ptt, mintime, timestep)
 @test ptt.groups[11].numpids == 2 # 73.5-74.25
 @test ptt.groups[18].numpids == 1

test/variable_rate.jl

@@ -1,5 +1,5 @@
 using DiffEqBase, JumpProcesses, OrdinaryDiffEq, StochasticDiffEq, Test
-using Random
+using Random, LinearSolve
 using StableRNGs
 rng = StableRNG(12345)
 
@@ -275,6 +275,16 @@ end
 # https://github.com/SciML/JumpProcesses.jl/issues/320
 # note that even with the seeded StableRNG this test is not 
 # deterministic for some reason.
+function getmean(Nsims, prob, alg, dt, tsave)
+    umean = zeros(length(tsave))
+    for i in 1:Nsims
+        sol = solve(prob, alg; saveat = dt)
+        umean .+= Array(sol(tsave; idxs = 1))
+    end
+    umean ./= Nsims
+    return umean
+end
+
 let
     rng = StableRNG(12345)
     b = 2.0
@@ -309,11 +319,8 @@ let
     sjm_prob = JumpProblem(ode_prob, b_jump, d_jump; rng)
     dt = 0.1
     tsave = range(tspan[1], tspan[2]; step = dt)
-    umean = zeros(length(tsave))
-    for i in 1:Nsims
-        sol = solve(sjm_prob, Tsit5(); saveat = dt)
-        umean .+= Array(sol(tsave; idxs = 1))
+    for alg in (Tsit5(), Rodas5P(linsolve = QRFactorization()))
+        umean = getmean(Nsims, sjm_prob, alg, dt, tsave)
+        @test all(abs.(umean .- n.(tsave)) .< 0.05 * n.(tsave))
     end
-    umean ./= Nsims
-    @test all(abs.(umean .- n.(tsave)) .< 0.05 * n.(tsave))
 end