Merge pull request #2079 from SciML/rosenbrock23

Alternative error control on Rosenbrock23 algebraic equations
SciML · Dec 14, 2023 · 146c873 · 146c873
2 parents 6561a2f + 3df1150
commit 146c873
Show file tree

Hide file tree

Showing 2 changed files with 37 additions and 4 deletions.
diff --git a/src/perform_step/rosenbrock_perform_step.jl b/src/perform_step/rosenbrock_perform_step.jl
@@ -117,11 +117,16 @@ end
             veck₃ = _vec(k₃)
             vectmp = _vec(tmp)
             @.. broadcast=false vectmp=ifelse(cache.algebraic_vars,
-                dto6 * (veck₁ - 2 * veck₂ + veck₃))
+                false, dto6 * (veck₁ - 2 * veck₂ + veck₃))
         end
         calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol,
             integrator.opts.reltol, integrator.opts.internalnorm, t)
         integrator.EEst = integrator.opts.internalnorm(atmp, t)
+
+        if mass_matrix !== I
+            @.. broadcast=false atmp = ifelse(cache.algebraic_vars,fsallast,false)/integrator.opts.abstol
+            integrator.EEst += integrator.opts.internalnorm(atmp, t)
+        end
     end
     cache.linsolve = linres.cache
 end
@@ -245,6 +250,19 @@ end
         calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol,
             integrator.opts.reltol, integrator.opts.internalnorm, t)
         integrator.EEst = integrator.opts.internalnorm(atmp, t)
+
+        if mass_matrix !== I
+            if integrator.opts.abstol isa AbstractVector
+                @inbounds @simd ivdep for i in 1:length(vectmp)
+                    vectmp[i] = ifelse(cache.algebraic_vars[i],fsallast[i],false)/integrator.opts.abstol[i]
+                end
+            else
+                @inbounds @simd ivdep for i in 1:length(vectmp)
+                    vectmp[i] = ifelse(cache.algebraic_vars[i],fsallast[i],false)/integrator.opts.abstol
+                end
+            end
+            integrator.EEst += integrator.opts.internalnorm(vectmp, t)
+        end
     end
     cache.linsolve = linres.cache
 end
@@ -337,6 +355,11 @@ end
         calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol,
             integrator.opts.reltol, integrator.opts.internalnorm, t)
         integrator.EEst = integrator.opts.internalnorm(atmp, t)
+
+        if mass_matrix !== I
+            @.. broadcast=false atmp = ifelse(cache.algebraic_vars,fsallast,false)/integrator.opts.abstol
+            integrator.EEst += integrator.opts.internalnorm(atmp, t)
+        end
     end
     cache.linsolve = linres.cache
 end
@@ -399,6 +422,11 @@ end
         atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol,
             integrator.opts.reltol, integrator.opts.internalnorm, t)
         integrator.EEst = integrator.opts.internalnorm(atmp, t)
+
+        if mass_matrix !== I
+            atmp = @. ifelse(!integrator.differential_vars,integrator.fsallast,false)./integrator.opts.abstol
+            integrator.EEst += integrator.opts.internalnorm(atmp, t)
+        end
     end
     integrator.k[1] = k₁
     integrator.k[2] = k₂
@@ -466,6 +494,11 @@ end
         atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol,
             integrator.opts.reltol, integrator.opts.internalnorm, t)
         integrator.EEst = integrator.opts.internalnorm(atmp, t)
+
+        if mass_matrix !== I
+            atmp = @. ifelse(!integrator.differential_vars,integrator.fsallast,false)./integrator.opts.abstol
+            integrator.EEst += integrator.opts.internalnorm(atmp, t)
+        end
     end
 
     integrator.k[1] = k₁

diff --git a/test/regression/hard_dae.jl b/test/regression/hard_dae.jl
@@ -216,13 +216,13 @@ cb = DiscreteCallback(condition, affect!)
 
 prob = ODEProblem(f, deepcopy(res.zero), (0, 20.0), deepcopy(p_inv))
 refsol = solve(prob, Rodas4(), saveat = 0.1, callback = cb, tstops = [1.0], reltol = 1e-12,
-    abstol = 1e-18)
+    abstol = 1e-17)
 
 for solver in (Rodas4, Rodas4P, Rodas5, Rodas5P, FBDF, QNDF, Rosenbrock23)
     @show solver
     prob = ODEProblem(f, deepcopy(res.zero), (0, 20.0), deepcopy(p_inv))
-    sol = solve(prob, solver(), saveat = 0.1, callback = cb, tstops = [1.0], reltol = 1e-12,
-        abstol = 1e-16)
+    sol = solve(prob, solver(), saveat = 0.1, callback = cb, tstops = [1.0], reltol = 1e-14,
+        abstol = 1e-14)
     @test sol.retcode == ReturnCode.Success
     @test sol.t[end] == 20.0
     @test maximum(sol - refsol) < 1e-11