fix utility tests

Author: ChrisRackauckas

Project.toml

@@ -41,7 +41,7 @@ FiniteDiff = "2"
 ForwardDiff = "0.10.3"
 MuladdMacro = "0.2.1"
 NLsolve = "4"
-OrdinaryDiffEq = "6.2"
+OrdinaryDiffEq = "6.4"
 RandomNumbers = "1.5.3"
 RecursiveArrayTools = "2"
 Reexport = "0.2, 1.0"

src/algorithms.jl

@@ -635,7 +635,7 @@ struct ImplicitEM{CS,AD,F,F2,P,FDT,ST,CJ,T2,Controller} <: StochasticDiffEqNewto
   symplectic::Bool
 end
 ImplicitEM(;chunk_size=0,autodiff=true,diff_type=Val{:central},
-                          standardtag = Val{true}(),concrete_jac = Val{true}(),
+                          standardtag = Val{true}(),concrete_jac = nothing,
                           precs = OrdinaryDiffEq.DEFAULT_PRECS,
                           linsolve=nothing,nlsolve=NLNewton(),
                           extrapolant=:constant,
@@ -670,7 +670,7 @@ struct ImplicitEulerHeun{CS,AD,F,P,FDT,ST,CJ,N,T2,Controller} <: StochasticDiffE
   symplectic::Bool
 end
 ImplicitEulerHeun(;chunk_size=0,autodiff=true,diff_type=Val{:central},
-                          standardtag = Val{true}(),concrete_jac = Val{true}(),
+                          standardtag = Val{true}(),concrete_jac = nothing,
                           precs = OrdinaryDiffEq.DEFAULT_PRECS,
                           linsolve=nothing,nlsolve=NLNewton(),
                           extrapolant=:constant,
@@ -706,7 +706,7 @@ struct ImplicitRKMil{CS,AD,F,P,FDT,ST,CJ,N,T2,Controller,interpretation} <: Stoc
   symplectic::Bool
 end
 ImplicitRKMil(;chunk_size=0,autodiff=true,diff_type=Val{:central},
-                          standardtag = Val{true}(),concrete_jac = Val{true}(),
+                          standardtag = Val{true}(),concrete_jac = nothing,
                           precs = OrdinaryDiffEq.DEFAULT_PRECS,
                           linsolve=nothing,nlsolve=NLNewton(),
                           extrapolant=:constant,
@@ -741,7 +741,7 @@ struct ISSEM{CS,AD,F,P,FDT,ST,CJ,N,T2,Controller} <: StochasticDiffEqNewtonAdapt
   symplectic::Bool
 end
 ISSEM(;chunk_size=0,autodiff=true,diff_type=Val{:central},
-                       standardtag = Val{true}(),concrete_jac = Val{true}(),
+                       standardtag = Val{true}(),concrete_jac = nothing,
                        precs = OrdinaryDiffEq.DEFAULT_PRECS,
                        linsolve=nothing,nlsolve=NLNewton(),
                        extrapolant=:constant,
@@ -776,7 +776,7 @@ struct ISSEulerHeun{CS,AD,F,P,FDT,ST,CJ,N,T2,Controller} <: StochasticDiffEqNewt
  symplectic::Bool
 end
 ISSEulerHeun(;chunk_size=0,autodiff=true,diff_type=Val{:central},
-                      standardtag = Val{true}(),concrete_jac = Val{true}(),
+                      standardtag = Val{true}(),concrete_jac = nothing,
                       precs = OrdinaryDiffEq.DEFAULT_PRECS,
                       linsolve=nothing,nlsolve=NLNewton(),
                       extrapolant=:constant,
@@ -808,7 +808,7 @@ struct SKenCarp{CS,AD,F,P,FDT,ST,CJ,N,T2,Controller} <: StochasticDiffEqNewtonAd
 end
 
 SKenCarp(;chunk_size=0,autodiff=true,diff_type=Val{:central},
-                   standardtag = Val{true}(),concrete_jac = Val{true}(),
+                   standardtag = Val{true}(),concrete_jac = nothing,
                    precs = OrdinaryDiffEq.DEFAULT_PRECS,
                    linsolve=nothing,nlsolve=NLNewton(),
                    smooth_est=true,extrapolant=:min_correct,

test/utility_tests.jl

@@ -1,29 +1,13 @@
-using StochasticDiffEq, LinearAlgebra, SparseArrays, Random, Test, DiffEqOperators
+using StochasticDiffEq, LinearAlgebra, SparseArrays, Random, LinearSolve, Test
 using StochasticDiffEq.OrdinaryDiffEq: WOperator, set_gamma!, calc_W!
 
 @testset "Derivative Utilities" begin
-  @testset "WOperator" begin
-    Random.seed!(123)
-    y = zeros(2); b = rand(2)
-    mm = rand(2, 2)
-    for _J in [rand(2, 2), Diagonal(rand(2))]
-      _Ws = [-mm + 2.0 * _J, -mm/2.0 + _J]
-      for inplace in (true, false), (_W, W_transform) in zip(_Ws, [false, true])
-        W = WOperator{inplace}(mm, 1.0, DiffEqArrayOperator(_J), b, transform=W_transform)
-        set_gamma!(W, 2.0)
-        @test convert(AbstractMatrix,W) ≈ _W
-        @test W * b ≈ _W * b
-        mul!(y, W, b); @test y ≈ _W * b
-      end
-    end
-  end
-
   @testset "calc_W!" begin
     A = [-1.0 0.0; 0.0 -0.5]; σ = [0.9 0.0; 0.0 0.8]
     mm = [2.0 0.0; 0.0 1.0]
     u0 = [1.0, 1.0]; tmp = zeros(2)
-    tspan = (0.0,1.0); dt = 0.01
-    concrete_W = -mm + dt * A
+    tspan = (0.0,1.0); dt = 0.01; dtgamma = 0.5dt
+    concrete_W = -mm + dtgamma * A
 
     # Out-of-place
     _f = (u,p,t) -> A*u; _g = (u,p,t) -> σ*u
@@ -33,7 +17,7 @@ using StochasticDiffEq.OrdinaryDiffEq: WOperator, set_gamma!, calc_W!
     prob = SDEProblem(fun, _g, u0, tspan)
     integrator = init(prob, ImplicitEM(theta=1); adaptive=false, dt=dt)
     W = integrator.cache.nlsolver.cache.W
-    calc_W!(W, integrator, integrator.cache.nlsolver, integrator.cache, dt, false)
+    calc_W!(W, integrator, integrator.cache.nlsolver, integrator.cache, dtgamma, false)
     @test convert(AbstractMatrix, W) ≈ concrete_W
     @test W \ u0 ≈ concrete_W \ u0
 
@@ -45,9 +29,17 @@ using StochasticDiffEq.OrdinaryDiffEq: WOperator, set_gamma!, calc_W!
     prob = SDEProblem(fun, _g, u0, tspan)
     integrator = init(prob, ImplicitEM(theta=1); adaptive=false, dt=dt)
     W = integrator.cache.nlsolver.cache.W
-    calc_W!(W, integrator, integrator.cache.nlsolver, integrator.cache, dt, false)
-    @test convert(AbstractMatrix, integrator.cache.nlsolver.cache.W) ≈ concrete_W
-    ldiv!(tmp, lu!(integrator.cache.nlsolver.cache.W), u0); @test tmp ≈ concrete_W \ u0
+    calc_W!(W, integrator, integrator.cache.nlsolver, integrator.cache, dtgamma, false)
+
+    # Did not update because it's an array operator
+    # We don't want to build Jacobians when we have operators!
+    @test convert(AbstractMatrix, integrator.cache.nlsolver.cache.W) != concrete_W
+    ldiv!(tmp, lu!(integrator.cache.nlsolver.cache.W), u0); @test tmp != concrete_W \ u0
+
+    # But jacobian2W! will update the cache
+    StochasticDiffEq.OrdinaryDiffEq.jacobian2W!(integrator.cache.nlsolver.cache.W._concrete_form, mm, dtgamma, integrator.cache.nlsolver.cache.W.J.A, false)
+    @test convert(AbstractMatrix, integrator.cache.nlsolver.cache.W) == concrete_W
+    ldiv!(tmp, lu!(integrator.cache.nlsolver.cache.W), u0); @test tmp == concrete_W \ u0
   end
 
   @testset "Implicit solver with lazy W" begin
@@ -69,7 +61,7 @@ using StochasticDiffEq.OrdinaryDiffEq: WOperator, set_gamma!, calc_W!
       Random.seed!(0); sol2 = solve(prob2, Alg(theta=1); adaptive=false, dt=0.01)
       @test sol1(1.0) ≈ sol2(1.0) rtol=1e-4
       Random.seed!(0); sol1_ip = solve(prob1_ip, Alg(theta=1); adaptive=false, dt=0.01)
-      Random.seed!(0); sol2_ip = solve(prob2_ip, Alg(theta=1,linsolve=LinSolveFactorize(lu)); adaptive=false, dt=0.01)
+      Random.seed!(0); sol2_ip = solve(prob2_ip, Alg(theta=1); adaptive=false, dt=0.01)
       @test sol1_ip(1.0) ≈ sol2_ip(1.0) rtol=1e-4
     end
 
@@ -85,7 +77,7 @@ using StochasticDiffEq.OrdinaryDiffEq: WOperator, set_gamma!, calc_W!
     Random.seed!(0); sol2 = solve(prob2, SKenCarp(); adaptive=false, dt=0.01)
     @test sol1(1.0) ≈ sol2(1.0) rtol=1e-4
     Random.seed!(0); sol1_ip = solve(prob1_ip, SKenCarp(); adaptive=false, dt=0.01)
-    Random.seed!(0); sol2_ip = solve(prob2_ip, SKenCarp(linsolve=LinSolveFactorize(lu)); adaptive=false, dt=0.01)
+    Random.seed!(0); sol2_ip = solve(prob2_ip, SKenCarp(); adaptive=false, dt=0.01)
     @test sol1_ip(1.0) ≈ sol2_ip(1.0) rtol=1e-3
   end
 end