Trim unused OrdinaryDiffEq code

Author: charleskawczynski

examples/hybrid/ode_config.jl

@@ -1,5 +1,4 @@
 import SciMLBase
-import OrdinaryDiffEq as ODE
 import ClimaTimeSteppers as CTS
 
 is_explicit_CTS_algo_type(alg_or_tableau) =
@@ -8,35 +7,13 @@ is_explicit_CTS_algo_type(alg_or_tableau) =
 is_imex_CTS_algo_type(alg_or_tableau) =
     alg_or_tableau <: CTS.IMEXARKAlgorithmName
 
-is_implicit_type(::typeof(ODE.IMEXEuler)) = true
-is_implicit_type(alg_or_tableau) =
-    alg_or_tableau <: Union{
-        ODE.OrdinaryDiffEqImplicitAlgorithm,
-        ODE.OrdinaryDiffEqAdaptiveImplicitAlgorithm,
-    } || is_imex_CTS_algo_type(alg_or_tableau)
-
-is_ordinary_diffeq_newton(::typeof(ODE.IMEXEuler)) = true
-is_ordinary_diffeq_newton(alg_or_tableau) =
-    alg_or_tableau <: Union{
-        ODE.OrdinaryDiffEqNewtonAlgorithm,
-        ODE.OrdinaryDiffEqNewtonAdaptiveAlgorithm,
-    }
-
 is_imex_CTS_algo(::CTS.IMEXAlgorithm) = true
 is_imex_CTS_algo(::SciMLBase.AbstractODEAlgorithm) = false
 
-is_implicit(::ODE.OrdinaryDiffEqImplicitAlgorithm) = true
-is_implicit(::ODE.OrdinaryDiffEqAdaptiveImplicitAlgorithm) = true
-is_implicit(ode_algo) = is_imex_CTS_algo(ode_algo)
-
-is_rosenbrock(::ODE.Rosenbrock23) = true
-is_rosenbrock(::ODE.Rosenbrock32) = true
-is_rosenbrock(::SciMLBase.AbstractODEAlgorithm) = false
-use_transform(ode_algo) =
-    !(is_imex_CTS_algo(ode_algo) || is_rosenbrock(ode_algo))
+use_transform(ode_algo) = !is_imex_CTS_algo(ode_algo)
 
 function jac_kwargs(ode_algo, Y, jacobi_flags)
-    if is_implicit(ode_algo)
+    if is_imex_CTS_algo(ode_algo)
         W = SchurComplementW(Y, use_transform(ode_algo), jacobi_flags)
         if use_transform(ode_algo)
             return (; jac_prototype = W, Wfact_t = Wfact!)
@@ -57,32 +34,17 @@ function ode_configuration(
         ode_name = split(ode_name, ".")[end]
     end
     ode_sym = Symbol(ode_name)
-    alg_or_tableau = if hasproperty(ODE, ode_sym)
-        @warn "apply_limiter flag is ignored for OrdinaryDiffEq algorithms"
-        getproperty(ODE, ode_sym)
-    else
-        getproperty(CTS, ode_sym)
-    end
+    alg_or_tableau = getproperty(CTS, ode_sym)
     @info "Using ODE config: `$alg_or_tableau`"
 
     if is_explicit_CTS_algo_type(alg_or_tableau)
         return CTS.ExplicitAlgorithm(alg_or_tableau())
-    elseif !is_implicit_type(alg_or_tableau)
+    elseif !is_imex_CTS_algo_type(alg_or_tableau)
         return alg_or_tableau()
-    elseif is_ordinary_diffeq_newton(alg_or_tableau)
-        if max_newton_iters == 1
-            error("OridinaryDiffEq requires at least 2 Newton iterations")
-        end
-        # κ like a relative tolerance; its default value in ODE is 0.01
-        nlsolve = ODE.NLNewton(;
-            κ = max_newton_iters == 2 ? Inf : 0.01,
-            max_iter = max_newton_iters,
-        )
-        return alg_or_tableau(; linsolve = linsolve!, nlsolve)
     elseif is_imex_CTS_algo_type(alg_or_tableau)
         newtons_method = CTS.NewtonsMethod(; max_iters = max_newton_iters)
         return CTS.IMEXAlgorithm(alg_or_tableau(), newtons_method)
     else
-        return alg_or_tableau(; linsolve = linsolve!)
+        error("Uncaught case")
     end
 end

examples/ordinary_diff_eq_bug_fixes.jl

@@ -1,21 +1,5 @@
-using OrdinaryDiffEq:
-    AbstractNLSolver,
-    isnewton,
-    DIRK,
-    update_W!,
-    initialize!,
-    Divergence,
-    get_new_W!,
-    initial_η,
-    compute_step!,
-    apply_step!,
-    Convergence,
-    isJcurrent,
-    postamble!,
-    TryAgain
 using SciMLBase: RECOMPILE_BY_DEFAULT
 
-import OrdinaryDiffEq: nlsolve!
 import SciMLBase: SplitFunction
 
 #=
@@ -40,90 +24,3 @@ function SplitFunction{iip}(f1, f2; kwargs...) where {iip}
         tgrad = f1.tgrad,
     )
 end
-
-#=
-Issue: nlsolve! has a goto statement that restarts the solver when the
-       conditions for DiffEqBase.SlowConvergence (denoted by TryAgain) are
-       satisfied. These conditions include !isJcurrent, which is always true
-       when a custom Wfact or Wfact_t is used because of how calc_W! is written.
-       So, when a custom Wfact or Wfact_t is used, nlsolve! sometimes enters an
-       infinite loop (though there are also situations where it converges after
-       several iterations of the loop). This is probably not the intended
-       behavior, so this patch limits the number of loop iterations to 100,
-       after which divergence is declared.
-=#
-function nlsolve!(
-    nlsolver::AbstractNLSolver,
-    integrator,
-    cache = nothing,
-    repeat_step = false,
-)
-    num_redos = 0
-    @label REDO
-    if isnewton(nlsolver)
-        cache === nothing && throw(
-            ArgumentError(
-                "cache is not passed to `nlsolve!` when using NLNewton",
-            ),
-        )
-        if nlsolver.method === DIRK
-            γW = nlsolver.γ * integrator.dt
-        else
-            γW = nlsolver.γ * integrator.dt / nlsolver.α
-        end
-        update_W!(nlsolver, integrator, cache, γW, repeat_step)
-    end
-
-    (; maxiters, κ, fast_convergence_cutoff) = nlsolver
-
-    initialize!(nlsolver, integrator)
-    nlsolver.status = Divergence
-    η = get_new_W!(nlsolver) ? initial_η(nlsolver, integrator) : nlsolver.ηold
-
-    local ndz
-    for iter in 1:maxiters
-        nlsolver.iter = iter
-
-        # compute next step and calculate norm of residuals
-        iter > 1 && (ndzprev = ndz)
-        ndz = compute_step!(nlsolver, integrator)
-        if !isfinite(ndz)
-            nlsolver.status = Divergence
-            break
-        end
-
-        # check divergence (not in initial step)
-        if iter > 1
-            θ = ndz / ndzprev
-
-            # divergence
-            if θ > 2
-                nlsolver.status = Divergence
-                break
-            end
-        end
-
-        apply_step!(nlsolver, integrator)
-
-        # check for convergence
-        iter > 1 && (η = θ / (1 - θ))
-        if (iter == 1 && ndz < 1e-5) ||
-           (iter > 1 && (η >= zero(η) && η * ndz < κ))
-            nlsolver.status = Convergence
-            break
-        end
-    end
-
-    if isnewton(nlsolver) &&
-       nlsolver.status == Divergence &&
-       !isJcurrent(nlsolver, integrator)
-        num_redos += 1
-        if num_redos < 100
-            nlsolver.status = TryAgain
-            @goto REDO
-        end
-    end
-
-    nlsolver.ηold = η
-    postamble!(nlsolver, integrator)
-end