set p in alg_cache

SciML · Aug 14, 2024 · 025900d · 025900d
1 parent 791f557
commit 025900d
Show file tree

Hide file tree

Showing 4 changed files with 35 additions and 31 deletions.
diff --git a/lib/OrdinaryDiffEqDifferentiation/src/linsolve_utils.jl b/lib/OrdinaryDiffEqDifferentiation/src/linsolve_utils.jl
@@ -11,8 +11,12 @@ function dolinsolve(integrator, linsolve; A = nothing, linu = nothing, b = nothi
     _alg = unwrap_alg(integrator, true)
     if !isnothing(A)
         (;du, u, p, t) = integrator
-        p = isnothing(integrator) ? nothing : (du, u, p, t)
-        reinit!(linsolve; A, p)
+        if isnothing(integrator)
+            reinit!(linsolve; A)
+        else
+            p = (du, u, p, t)
+            reinit!(linsolve; A, p)
+        end
     end
 
     linres = solve!(linsolve; reltol)

diff --git a/lib/OrdinaryDiffEqNonlinearSolve/src/utils.jl b/lib/OrdinaryDiffEqNonlinearSolve/src/utils.jl
@@ -191,12 +191,12 @@ function build_nlsolver(
             end
             jac_config = build_jac_config(alg, nf, uf, du1, uprev, u, ztmp, dz)
         end
-        linprob = LinearProblem(W, _vec(k); u0 = _vec(dz))
+        linprob = LinearProblem(W, _vec(k), (isdae ? du1 : nothing,u,p,t); u0 = _vec(dz))
         Pl, Pr = wrapprecs(
             alg.precs(W, nothing, u, p, t, nothing, nothing, nothing,
                 nothing)...,
             weight, dz)
-        linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
+        linsolve = init(linprob, alg.linsolve, (isdae ? du1 : nothing,u,p,t); alias_A = true, alias_b = true,
             Pl = Pl, Pr = Pr,
             assumptions = LinearSolve.OperatorAssumptions(true))
 

diff --git a/lib/OrdinaryDiffEqRosenbrock/src/generic_rosenbrock.jl b/lib/OrdinaryDiffEqRosenbrock/src/generic_rosenbrock.jl
@@ -247,7 +247,7 @@ function gen_algcache(cacheexpr::Expr,constcachename::Symbol,algname::Symbol,tab
             tf = TimeGradientWrapper(f,uprev,p)
             uf = UJacobianWrapper(f,t,p)
             linsolve_tmp = zero(rate_prototype)
-            linprob = LinearProblem(W,_vec(linsolve_tmp); u0=_vec(tmp))
+            linprob = LinearProblem(W,_vec(linsolve_tmp), (nothing, u, p, t); u0=_vec(tmp))
             linsolve = init(linprob,alg.linsolve,alias_A=true,alias_b=true,
                             Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))),
                             Pr = Diagonal(_vec(weight)))
@@ -995,7 +995,7 @@ references = """
 """,
 "Rodas3",
 references = """
-- Steinebach G. Construction of Rosenbrock–Wanner method Rodas5P and numerical benchmarks 
+- Steinebach G. Construction of Rosenbrock–Wanner method Rodas5P and numerical benchmarks
   within the Julia Differential Equations package.
   In: BIT Numerical Mathematics, 63(2), 2023
 """,
@@ -1056,7 +1056,7 @@ lower if not corrected).
 """,
 "Rodas4P",
 references = """
-- Steinebach G., Rodas23W / Rodas32P - a Rosenbrock-type method for DAEs with additional error estimate 
+- Steinebach G., Rodas23W / Rodas32P - a Rosenbrock-type method for DAEs with additional error estimate
   for dense output and Julia implementation,
   In progress.
 """,
@@ -1070,7 +1070,7 @@ of Roadas4P and in case of inexact Jacobians a second order W method.
 """,
 "Rodas4P2",
 references = """
-- Steinebach G., Rodas23W / Rodas32P - a Rosenbrock-type method for DAEs with additional error estimate 
+- Steinebach G., Rodas23W / Rodas32P - a Rosenbrock-type method for DAEs with additional error estimate
   for dense output and Julia implementation,
   In progress.
 """,
@@ -1095,7 +1095,7 @@ Has improved stability in the adaptive time stepping embedding.
 """,
 "Rodas5P",
 references = """
-- Steinebach G. Construction of Rosenbrock–Wanner method Rodas5P and numerical benchmarks 
+- Steinebach G. Construction of Rosenbrock–Wanner method Rodas5P and numerical benchmarks
   within the Julia Differential Equations package.
   In: BIT Numerical Mathematics, 63(2), 2023
 """,
@@ -1225,10 +1225,10 @@ references = """
     ROS2PRTableau()
 
 2nd order stiffly accurate Rosenbrock method with 3 internal stages with (Rinf=0).
-For problems with medium stiffness the convergence behaviour is very poor and it is recommended to use 
+For problems with medium stiffness the convergence behaviour is very poor and it is recommended to use
 [`ROS2S`](@ref) instead.
 
-Rang, Joachim (2014): The Prothero and Robinson example: 
+Rang, Joachim (2014): The Prothero and Robinson example:
 Convergence studies for Runge-Kutta and Rosenbrock-Wanner methods. https://doi.org/10.24355/dbbs.084-201408121139-0
 """
 function ROS2PRTableau() # 2nd order
@@ -1248,7 +1248,7 @@ end
 @doc rosenbrock_docstring(
 """
 2nd order stiffly accurate Rosenbrock method with 3 internal stages with (Rinf=0).
-For problems with medium stiffness the convergence behaviour is very poor and it is recommended to use 
+For problems with medium stiffness the convergence behaviour is very poor and it is recommended to use
 [`ROS2S`](@ref) instead.
 """,
 "ROS2PR",
@@ -1267,7 +1267,7 @@ ROS2PR
 2nd order stiffly accurate Rosenbrock-Wanner W-method with 3 internal stages with B_PR consistent of order 2 with (Rinf=0).
 More Information at https://doi.org/10.24355/dbbs.084-201408121139-0
 
-Rang, Joachim (2014): The Prothero and Robinson example: 
+Rang, Joachim (2014): The Prothero and Robinson example:
 Convergence studies for Runge-Kutta and Rosenbrock-Wanner methods. https://doi.org/10.24355/dbbs.084-201408121139-0
 """
 function ROS2STableau() # 2nd order
@@ -1335,7 +1335,7 @@ references = """
 
 3nd order stiffly accurate Rosenbrock-Wanner method with 3 internal stages with B_PR consistent of order 3, which is strongly A-stable with Rinf~=-0.73.
 
-Rang, Joachim (2014): The Prothero and Robinson example: 
+Rang, Joachim (2014): The Prothero and Robinson example:
 Convergence studies for Runge-Kutta and Rosenbrock-Wanner methods. https://doi.org/10.24355/dbbs.084-201408121139-0
 """
 function ROS3PRTableau() # 3rd order
@@ -1371,7 +1371,7 @@ references = """
 3nd order stiffly accurate Rosenbrock method with 3 internal stages with B_PR consistent of order 3, which is strongly A-stable with Rinf~=-0.73
 Convergence with order 4 for the stiff case, but has a poor accuracy.
 
-Rang, Joachim (2014): The Prothero and Robinson example: 
+Rang, Joachim (2014): The Prothero and Robinson example:
 Convergence studies for Runge-Kutta and Rosenbrock-Wanner methods. https://doi.org/10.24355/dbbs.084-201408121139-0
 """
 function Scholz4_7Tableau() # 3rd order
@@ -1599,7 +1599,7 @@ references = """
 B_PR consistent of order 2 with Rinf=0.
 The order of convergence decreases if medium stiff problems are considered, but it has good results for very stiff cases.
 
-Rang, Joachim (2014): The Prothero and Robinson example: 
+Rang, Joachim (2014): The Prothero and Robinson example:
 Convergence studies for Runge-Kutta and Rosenbrock-Wanner methods. https://doi.org/10.24355/dbbs.084-201408121139-0
 """
 function ROS3PRLTableau() # 3rd order
@@ -1639,7 +1639,7 @@ references = """
 B_PR consistent of order 3.
 The order of convergence does NOT decreases if medium stiff problems are considered as it does for [`ROS3PRL`](@ref).
 
-Rang, Joachim (2014): The Prothero and Robinson example: 
+Rang, Joachim (2014): The Prothero and Robinson example:
 Convergence studies for Runge-Kutta and Rosenbrock-Wanner methods. https://doi.org/10.24355/dbbs.084-201408121139-0
 """
 function ROS3PRL2Tableau() # 3rd order
@@ -1677,7 +1677,7 @@ references = """
 
 4rd order L-stable Rosenbrock-Krylov method with 4 internal stages,
 with a 3rd order embedded method which is strongly A-stable with Rinf~=0.55. (when using exact Jacobians)
-Tranquilli, Paul and Sandu, Adrian (2014): Rosenbrock--Krylov Methods for Large Systems of Differential Equations 
+Tranquilli, Paul and Sandu, Adrian (2014): Rosenbrock--Krylov Methods for Large Systems of Differential Equations
 https://doi.org/10.1137/130923336
 """
 function ROK4aTableau() # 4rd order
@@ -1703,7 +1703,7 @@ with a 3rd order embedded method which is strongly A-stable with Rinf~=0.55. (wh
 """,
 "ROK4a",
 references = """
-- Tranquilli, Paul and Sandu, Adrian (2014): 
+- Tranquilli, Paul and Sandu, Adrian (2014):
   Rosenbrock--Krylov Methods for Large Systems of Differential Equations
   https://doi.org/10.1137/130923336
 """) ROK4a

diff --git a/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl b/lib/OrdinaryDiffEqRosenbrock/src/rosenbrock_caches.jl
@@ -94,7 +94,7 @@ function alg_cache(alg::Rosenbrock23, u, rate_prototype, ::Type{uEltypeNoUnits},
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
 
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
     linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
 
@@ -135,7 +135,7 @@ function alg_cache(alg::Rosenbrock32, u, rate_prototype, ::Type{uEltypeNoUnits},
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
 
     linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
@@ -280,7 +280,7 @@ function alg_cache(alg::ROS3P, u, rate_prototype, ::Type{uEltypeNoUnits},
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
     linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     grad_config = build_grad_config(alg, f, tf, du1, t)
@@ -362,7 +362,7 @@ function alg_cache(alg::Rodas3, u, rate_prototype, ::Type{uEltypeNoUnits},
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
     linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     grad_config = build_grad_config(alg, f, tf, du1, t)
@@ -551,7 +551,7 @@ function alg_cache(alg::Rodas23W, u, rate_prototype, ::Type{uEltypeNoUnits},
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
     linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     grad_config = build_grad_config(alg, f, tf, du1, t)
@@ -591,7 +591,7 @@ function alg_cache(alg::Rodas3P, u, rate_prototype, ::Type{uEltypeNoUnits},
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
     linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     grad_config = build_grad_config(alg, f, tf, du1, t)
@@ -710,7 +710,7 @@ function alg_cache(alg::Rodas4, u, rate_prototype, ::Type{uEltypeNoUnits},
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
     linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     grad_config = build_grad_config(alg, f, tf, du1, t)
@@ -766,7 +766,7 @@ function alg_cache(alg::Rodas42, u, rate_prototype, ::Type{uEltypeNoUnits},
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
     linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     grad_config = build_grad_config(alg, f, tf, du1, t)
@@ -822,7 +822,7 @@ function alg_cache(alg::Rodas4P, u, rate_prototype, ::Type{uEltypeNoUnits},
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
     linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     grad_config = build_grad_config(alg, f, tf, du1, t)
@@ -878,7 +878,7 @@ function alg_cache(alg::Rodas4P2, u, rate_prototype, ::Type{uEltypeNoUnits},
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
     linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     grad_config = build_grad_config(alg, f, tf, du1, t)
@@ -991,7 +991,7 @@ function alg_cache(alg::Rodas5, u, rate_prototype, ::Type{uEltypeNoUnits},
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
     linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     grad_config = build_grad_config(alg, f, tf, du1, t)
@@ -1051,7 +1051,7 @@ function alg_cache(
     tf = TimeGradientWrapper(f, uprev, p)
     uf = UJacobianWrapper(f, t, p)
     linsolve_tmp = zero(rate_prototype)
-    linprob = LinearProblem(W, _vec(linsolve_tmp); u0 = _vec(tmp))
+    linprob = LinearProblem(W, _vec(linsolve_tmp), (nothing,u,p,t); u0 = _vec(tmp))
     linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
         assumptions = LinearSolve.OperatorAssumptions(true))
     grad_config = build_grad_config(alg, f, tf, du1, t)