try this

SciML · Sep 10, 2024 · aeafc10 · aeafc10
1 parent 81559b2
commit aeafc10
Show file tree

Hide file tree

Showing 2 changed files with 7 additions and 28 deletions.
diff --git a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl
@@ -335,7 +335,7 @@ function Base.convert(::Type{AbstractMatrix}, W::WOperator{IIP}) where {IIP}
     else
         # Non-allocating already updated
         #_W = W._concrete_form
-        #jacobian2W!(_W, W.mass_matrix, W.gamma, convert(AbstractMatrix, W.J), W.transform)
+        #jacobian2W!(_W, W.mass_matrix, W.gamma, W.J, W.transform)
     end
     return W._concrete_form
 end
@@ -392,10 +392,6 @@ function Base.:\(W::WOperator, x::Number)
     end
 end
 
-function LinearSolve.defaultalg(W::WOperator, b, assump::OperatorAssumptions{Bool})
-    LinearSolve.defaultalg(W.J, b, assump)
-end
-
 function LinearAlgebra.mul!(Y::AbstractVecOrMat, W::WOperator, B::AbstractVecOrMat)
     if W.transform
         # Compute mass_matrix * B
@@ -455,11 +451,7 @@ function do_newJW(integrator, alg, nlsolver, repeat_step)::NTuple{2, Bool}
     integrator.iter <= 1 && return true, true # at least one JW eval at the start
     repeat_step && return false, false
     islin, _ = islinearfunction(integrator)
-    if islin
-        # no further J eval when it's linear
-        # W eval needed if not adaptive
-        return false, !integrator.opts.adaptive
-    end
+    islin && return false, false # no further JW eval when it's linear
     !integrator.opts.adaptive && return true, true # Not adaptive will always refactorize
     errorfail = integrator.EEst > one(integrator.EEst)
     if alg isa DAEAlgorithm
@@ -708,22 +700,12 @@ function calc_W!(W, integrator, nlsolver::Union{Nothing, AbstractNLSolver}, cach
         else
             update_coefficients!(W, uprev, p, t; transform = W_transform, dtgamma)
         end
-        if W.J !== nothing && (W.J isa MatrixOperator || !(W.J isa AbstractSciMLOperator))
+        if W.J !== nothing && !(W.J isa AbstractSciMLOperator)
             islin, isode = islinearfunction(integrator)
-            if islin
-                J = isode ? f.f : f.f1.f
-            elseif new_jac
-                calc_J!(W.J, integrator, lcache, next_step)
-            end
-            if W.J isa MatrixOperator
-                J = J.A
-            end
-            if mass_matrix isa MatrixOperator
-                mass_matrix = mass_matrix.A
-            end
-            if new_W && !isdae
+            islin ? (J = isode ? f.f : f.f1.f) :
+            (new_jac && (calc_J!(W.J, integrator, lcache, next_step)))
+            new_W && !isdae &&
                 jacobian2W!(W._concrete_form, mass_matrix, dtgamma, J, W_transform)
-            end
         end
     elseif W isa AbstractSciMLOperator && !(W isa StaticWOperator)
         update_coefficients!(W, uprev, p, t; transform = W_transform, dtgamma)
@@ -929,7 +911,6 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits},
         else
             deepcopy(f.jac_prototype)
         end
-
         W = if alg isa DAEAlgorithm
             J
         elseif IIP

diff --git a/test/interface/utility_tests.jl b/test/interface/utility_tests.jl
@@ -60,9 +60,7 @@ end
         jac_prototype = MatrixOperator(similar(A);
             update_func! = (J, u, p, t) -> J .= t .*
                                                 A))
-
-    for Alg in [ImplicitEuler, Rosenbrock23, Rodas5]
-        println(Alg)
+    @testset "$Alg" for Alg in [ImplicitEuler, Rosenbrock23, Rodas5]
         sol1 = solve(ODEProblem(fun1, u0, tspan), Alg(); adaptive = false, dt = 0.01)
         sol2 = solve(ODEProblem(fun2, u0, tspan), Alg(); adaptive = false, dt = 0.01)
         @test sol1(1.0)≈sol2(1.0) atol=1e-3