From aeafc103b5ee157d2ce4e8c300927eb705983512 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Tue, 10 Sep 2024 15:44:39 -0400 Subject: [PATCH] try this --- .../src/derivative_utils.jl | 31 ++++--------------- test/interface/utility_tests.jl | 4 +-- 2 files changed, 7 insertions(+), 28 deletions(-) diff --git a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl index 6c30602ea7..b9e3695b69 100644 --- 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 index d0dbf2427b..f95e3aebd5 100644 --- 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