From 426483cdf135bdd74989fbd40d717ecd5007ccc4 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Thu, 10 Oct 2024 15:28:40 -0400 Subject: [PATCH 1/4] don't use jacvec with concrete_jac --- lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl index abf994e5bc..879921902b 100644 --- a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl +++ b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl @@ -734,9 +734,7 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, else (u, p, t) -> _f(u, p, t) end - jacvec = JacVec(__f, copy(u), p, t; - autodiff = alg_autodiff(alg), tag = OrdinaryDiffEqTag()) - WOperator{IIP}(f.mass_matrix, dt, J, u, jacvec) + WOperator{IIP}(f.mass_matrix, dt, J, u, nothing) end else J = if !IIP && DiffEqBase.has_jac(f) From 8dba9c5b189d00dae4beb06415ddc592251f469b Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Thu, 10 Oct 2024 15:29:51 -0400 Subject: [PATCH 2/4] update comment --- lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl index 879921902b..7c62f56f97 100644 --- a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl +++ b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl @@ -719,7 +719,7 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, concrete_jac(alg) !== nothing && concrete_jac(alg) # The linear solver does not need a concrete Jacobian, but the user has # asked for one. This will happen when the Jacobian is used in the preconditioner - # Thus setup JacVec and a concrete J, using sparsity when possible + # or when jvp computation is expensive. Therefore, use concrete J, and sparsity when possible _f = islin ? (isode ? f.f : f.f1.f) : f J = if f.jac_prototype === nothing ArrayInterface.undefmatrix(u) From 83cc7cec9e23fa6f51793a544bac02dcae2fdd05 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Thu, 10 Oct 2024 15:57:27 -0400 Subject: [PATCH 3/4] use non SciMLOperator W_prototype --- lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl index 7c62f56f97..240559bae9 100644 --- a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl +++ b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl @@ -391,7 +391,7 @@ function do_newJW(integrator, alg, nlsolver, repeat_step)::NTuple{2, Bool} # TODO: add `isJcurrent` support for Rosenbrock solvers if !isnewton(nlsolver) isfreshJ = !(integrator.alg isa CompositeAlgorithm) && - (integrator.iter > 1 && errorfail && !integrator.u_modified) + (!integrator.u_modified) return !isfreshJ, true end isfirstcall(nlsolver) && return true, true @@ -684,9 +684,9 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, # TODO - if jvp given, make it SciMLOperators.FunctionOperator # TODO - make mass matrix a SciMLOperator so it can be updated with time. Default to IdentityOperator islin, isode = islinearfunction(f, alg) - if isdefined(f, :W_prototype) && (f.W_prototype isa AbstractSciMLOperator) - # We use W_prototype when it is provided as a SciMLOperator, and in this case we require jac_prototype to be a SciMLOperator too. - if !(f.jac_prototype isa AbstractSciMLOperator) + if isdefined(f, :W_prototype) && !isnothing(f.W_prototype) + # If W_prototype is a SciMLOperator, we require jac_prototype to be a SciMLOperator too. + if f.W_prototype isa AbstractSciMLOperator && !(f.jac_prototype isa AbstractSciMLOperator) error("SciMLOperator for W_prototype only supported when jac_prototype is a SciMLOperator, but got $(typeof(f.jac_prototype))") end W = f.W_prototype From 851e3aa19c37e657ad9362d3a2415a168ecffd27 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Thu, 10 Oct 2024 16:00:56 -0400 Subject: [PATCH 4/4] simplify logic slightly --- lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl index 240559bae9..6f98243a65 100644 --- a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl +++ b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl @@ -703,7 +703,7 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, # If factorization, then just use the jac_prototype J = similar(f.jac_prototype) W = similar(J) - elseif (IIP && (concrete_jac(alg) === nothing || !concrete_jac(alg)) && + elseif (IIP && (concrete_jac(alg) != true) && alg.linsolve !== nothing && !LinearSolve.needs_concrete_A(alg.linsolve)) # If the user has chosen GMRES but no sparse Jacobian, assume that the dense @@ -716,7 +716,7 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits}, J = jacvec W = WOperator{IIP}(f.mass_matrix, dt, J, u, jacvec) elseif alg.linsolve !== nothing && !LinearSolve.needs_concrete_A(alg.linsolve) || - concrete_jac(alg) !== nothing && concrete_jac(alg) + concrete_jac(alg) == true # The linear solver does not need a concrete Jacobian, but the user has # asked for one. This will happen when the Jacobian is used in the preconditioner # or when jvp computation is expensive. Therefore, use concrete J, and sparsity when possible