From 895c0c1b662b766b4df0c34cd01c7fac1d906aa9 Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Mon, 3 Aug 2026 15:43:35 -0400
Subject: [PATCH 1/2] logging for isoutofdomain and add interface for
 domain_checks

---
 lib/DiffEqBase/src/DiffEqBase.jl              |  10 +-
 lib/DiffEqBase/src/domain_checks.jl           | 253 ++++++++++++++++++
 lib/DiffEqBase/src/solve.jl                   |  28 ++
 .../src/OrdinaryDiffEqCore.jl                 |   7 +-
 .../src/integrators/integrator_utils.jl       |  22 +-
 lib/OrdinaryDiffEqCore/src/solve.jl           |  29 ++
 src/OrdinaryDiffEq.jl                         |   4 +-
 7 files changed, 349 insertions(+), 4 deletions(-)
 create mode 100644 lib/DiffEqBase/src/domain_checks.jl

diff --git a/lib/DiffEqBase/src/DiffEqBase.jl b/lib/DiffEqBase/src/DiffEqBase.jl
index d94b4d266a9..086cfe60f0e 100644
--- a/lib/DiffEqBase/src/DiffEqBase.jl
+++ b/lib/DiffEqBase/src/DiffEqBase.jl
@@ -153,6 +153,7 @@ include("dae_initialization.jl")
 
 include("callbacks.jl")
 include("common_defaults.jl")
+include("domain_checks.jl")
 include("solve.jl")
 include("internal_euler.jl")
 include("norecompile.jl")
@@ -194,6 +195,9 @@ export SensitivityADPassThrough
 
 export AutoDePSpecialize
 
+# written at `solve` call sites.
+export @isoutofdomain
+
 # Declare DiffEqBase-owned, documented API names `public` so downstream packages can
 # drop their `DiffEqBase.X` non-public ExplicitImports ignores. The `public` keyword is
 # only parseable on Julia >= 1.11.0-DEV.469, so it is gated to keep the 1.10 floor parsing.
@@ -220,7 +224,11 @@ export AutoDePSpecialize
             :prepare_alg, :prob2dtmin, :timedepentdtmin, :check_prob_alg_pairing,
             :default_factorize, :stripunits,
             # Solver-author wrapper/tag types and convergence-testing entry type
-            :EvalFunc, :OrdinaryDiffEqTag, :ConvergenceSetup
+            :EvalFunc, :OrdinaryDiffEqTag, :ConvergenceSetup,
+            # Domain-violation mechanism (`domain_checks` solver keyword)
+            :DomainCheckedFunction, :domain_checks_failing,
+            # Source-retaining `isoutofdomain` predicates, for failure diagnostics
+            :TracedPredicate, :TracedPredicateLeaf, :isoutofdomain_report
         )
     )
 end
diff --git a/lib/DiffEqBase/src/domain_checks.jl b/lib/DiffEqBase/src/domain_checks.jl
new file mode 100644
index 00000000000..61a5d1ac134
--- /dev/null
+++ b/lib/DiffEqBase/src/domain_checks.jl
@@ -0,0 +1,253 @@
+"""
+    domain_checks_failing(checks, u, p, t)
+
+Evaluate `checks` against `(u, p, t)` and return one message per failing check, naming each by
+its position and by its source expression when it has one, or `nothing` if all pass. Every
+check is evaluated rather than stopping at the first failure; the all-passing case allocates
+nothing. Each check is `(u, p, t) -> Bool` or
+`(u, p, t) -> Pair{Bool, <:AbstractString}`; a plain `Bool` uses a generic default failure
+message. Used both by [`DomainCheckedFunction`](@ref) (to decide whether to call the real
+right-hand-side) and by diagnostics (to explain which predicate failed).
+"""
+@inline function domain_checks_failing(checks, u, p, t)
+    failures = nothing
+    for (idx, check) in enumerate(checks)
+        result = check(u, p, t)
+        ok, msg = result isa Pair ? (result.first, result.second) : (result, "Domain check failed")
+        ok && continue
+        line = check isa TracedPredicate ? "[$idx] $msg: $(check.src)" : "[$idx] $msg"
+        failures === nothing ? (failures = [line]) : push!(failures, line)
+    end
+    return failures
+end
+
+"""
+    DomainCheckedFunction{iip, F, C}(f, pre_checks)
+
+Wraps a raw ODE right-hand-side callable `f` together with a tuple/vector of `pre_checks`
+predicates, each of the form `(u, p, t) -> Bool` or `(u, p, t) -> Pair{Bool, <:AbstractString}`.
+Before calling `f`, every predicate is evaluated against `u`. If any predicate fails, `f` is
+not called at all, and instead, a `NaN` is injected into the state. A failing check therefore 
+falls into the NaN-detection in the adaptive step controller
+(`integrator.isout = ... || isnan(EEst) || isinf(EEst)`) straight into a
+shrink-and-retry.
+
+Constructed automatically from the `domain_checks` solver keyword.
+"""
+struct DomainCheckedFunction{iip, F, C}
+    f::F
+    pre_checks::C
+end
+
+function DomainCheckedFunction{iip}(f::F, pre_checks::C) where {iip, F, C}
+    return DomainCheckedFunction{iip, F, C}(f, pre_checks)
+end
+
+function (dcf::DomainCheckedFunction{true})(du, u, p, t)
+    if domain_checks_failing(dcf.pre_checks, u, p, t) === nothing
+        dcf.f(du, u, p, t)
+    else
+        du .= NaN
+    end
+    return nothing
+end
+
+function (dcf::DomainCheckedFunction{false})(u, p, t)
+    domain_checks_failing(dcf.pre_checks, u, p, t) === nothing && return dcf.f(u, p, t)
+    #maintain units of u
+    return u .* (zero(eltype(u)) / zero(eltype(u)))
+end
+
+
+#part 2: deal with isoutofdomain for logging purposes
+"""
+    TracedPredicateLeaf(f, src, operand_srcs, operand_f)
+
+One atomic sub-expression of an `isoutofdomain` predicate built by [`@isoutofdomain`](@ref),
+retaining its source text so diagnostics can name it.
+
+# Fields
+- `f`: `(u, p, t) -> Bool` evaluating just this sub-expression.
+- `src`: the sub-expression's source text, e.g. `"u[1] < 0"`.
+- `operand_srcs`: source text of each operand, e.g. `("u[1]", "0")`. Empty when the leaf is
+  not a comparison.
+- `operand_f`: `(u, p, t) -> Tuple` of the operand *values*, or `nothing` when the leaf is
+  not a comparison and so has no operands to report.
+"""
+struct TracedPredicateLeaf{F, S, O}
+    f::F
+    src::String
+    operand_srcs::S
+    operand_f::O
+end
+
+"""
+    TracedPredicate(f, src, leaves)
+
+An `isoutofdomain` predicate that carries its own source text, built by
+[`@isoutofdomain`](@ref).
+
+Purely a logging device: `f` is the user's original closure, stored verbatim and called
+unchanged, so this costs nothing on the integrator's hot path (a single-field forward through
+an immutable struct is erased by the compiler — the emitted code is instruction-for-instruction
+identical to calling `f` directly). `src` and `leaves` are touched only by
+[`isoutofdomain_report`](@ref) when a solve fails, so the tree never participates in deciding
+whether a step is rejected.
+
+It has to be the object stored in `opts.isoutofdomain` because that is the only handle
+diagnostics have to reach the metadata from an integrator.
+"""
+struct TracedPredicate{F, L}
+    f::F
+    src::String
+    leaves::L
+end
+
+@inline (tp::TracedPredicate)(u, p, t) = tp.f(u, p, t)
+
+const COMPARISON_OPS = (:<, :>, :<=, :>=, :(==), :!=, :≤, :≥, :≠, :isless)
+
+function is_connective(ex)
+    ex isa Expr || return false
+    (ex.head === :|| || ex.head === :&&) && return true
+    ex.head === :call || return false
+    length(ex.args) == 3 && ex.args[1] in (:|, :&) && return true
+    return length(ex.args) == 2 && ex.args[1] === :!
+end
+
+# Strip LineNumberNodes and single-statement `begin ... end` wrappers, which the parser
+# inserts around lambda bodies.
+function strip_lines(ex)
+    ex isa Expr || return ex
+    if ex.head === :block
+        stmts = filter(a -> !(a isa LineNumberNode), ex.args)
+        length(stmts) == 1 && return strip_lines(only(stmts))
+    end
+    return Expr(ex.head, map(strip_lines, ex.args)...)
+end
+
+# Flatten a boolean expression into its atomic sub-expressions, in source order.
+function collect_leaves!(acc, ex)
+    if is_connective(ex)
+        operands = (ex.head === :call) ? ex.args[2:end] : ex.args
+        for operand in operands
+            collect_leaves!(acc, operand)
+        end
+    else
+        push!(acc, ex)
+    end
+    return acc
+end
+
+# Operands of a comparison leaf worth reporting the value of, or `nothing` if the leaf isn't a
+# comparison or has no such operands.
+function leaf_operands(ex)
+    ex isa Expr || return nothing
+    operands = if ex.head === :call && length(ex.args) == 3 && ex.args[1] in COMPARISON_OPS
+        (ex.args[2], ex.args[3])
+    elseif ex.head === :comparison
+        # `0 < u[1] < 1` => (0, :<, u[1], :<, 1); operands sit at the odd indices
+        Tuple(ex.args[1:2:end])
+    else
+        return nothing
+    end
+    reportable = filter(o -> !(o isa Union{Number, Bool, AbstractString, Char}), operands)
+    return isempty(reportable) ? nothing : Tuple(reportable)
+end
+
+"""
+    @isoutofdomain (u, p, t) -> expr
+
+Build a [`TracedPredicate`](@ref) from an `isoutofdomain` predicate, retaining the source
+text of the whole expression and of each atomic sub-expression so that a failed solve can
+report *which* part of the predicate rejected the step, and with what operand values:
+
+```julia
+sol = solve(prob, Tsit5();
+    isoutofdomain = @isoutofdomain (u, p, t) -> u[1] < 0 || u[2] > 10)
+```
+
+On a `dt <= dtmin` abort the diagnostics then read:
+
+```
+predicate: u[1] < 0 || u[2] > 10
+  [1] u[1] < 0  => false  (u[1] = 0.42)
+  [2] u[2] > 10 => true   (u[2] = 14.7)
+```
+
+The predicate is stored and called verbatim, so this is free at run time; the retained source
+is read only when a solve fails. 
+
+The argument must be an anonymous function of three parameters; their names are taken from
+the signature, so `(state, par, time) -> state[1] < 0` reports `state[1] < 0`. Collection
+recurses through `||`, `&&`, `!`, `|` and `&`; anything else (e.g. `any(x -> x < 0, u)`)
+becomes a single leaf, still labelled with its source text but without operand decomposition.
+"""
+macro isoutofdomain(ex)
+    if !(ex isa Expr && ex.head === :->)
+        throw(ArgumentError("@isoutofdomain expects an anonymous function, e.g. `@isoutofdomain (u, p, t) -> u[1] < 0`; got `$ex`"))
+    end
+    signature = ex.args[1]
+    argnames = signature isa Expr && signature.head === :tuple ? signature.args : [signature]
+    if length(argnames) != 3 || !all(a -> a isa Symbol, argnames)
+        throw(ArgumentError("@isoutofdomain expects exactly three plain parameters `(u, p, t)`; got `$signature`"))
+    end
+
+    body = strip_lines(ex.args[2])
+    argtuple = Expr(:tuple, argnames...)
+
+    leaf_exprs = map(collect_leaves!(Any[], body)) do leaf
+        operands = leaf_operands(leaf)
+        operand_srcs = operands === nothing ? () : Tuple(string.(operands))
+        operand_f = operands === nothing ? nothing :
+            esc(Expr(:->, argtuple, Expr(:tuple, operands...)))
+        # `esc` so the sub-expression's free variables resolve in the caller's scope: a
+        # predicate like `u[1] < lo` must find the caller's `lo`, not DiffEqBase's.
+        leaf_f = esc(Expr(:->, argtuple, leaf))
+        return :($TracedPredicateLeaf($leaf_f, $(string(leaf)), $operand_srcs, $operand_f))
+    end
+
+    return :($TracedPredicate($(esc(ex)), $(string(body)), ($(leaf_exprs...),)))
+end
+
+format_operand(v::Real) = @sprintf("%.4g", v)
+format_operand(v) = repr(v)
+
+function explain_leaf(idx, leaf, u, p, t, width)
+    label = "  [$idx] " * rpad(leaf.src, width) * " => "
+    value = try
+        leaf.f(u, p, t)
+    catch err
+        return label * "errored: " * string(typeof(err))
+    end
+    operands = ""
+    if leaf.operand_f !== nothing
+        operands = try
+            values = leaf.operand_f(u, p, t)
+            pairs = ("$src = $(format_operand(val))" for (src, val) in zip(leaf.operand_srcs, values))
+            "  (" * join(pairs, ", ") * ")"
+        catch err
+            "  (operands unavailable: " * string(typeof(err)) * ")"
+        end
+    end
+    return label * string(value) * operands
+end
+
+"""
+    isoutofdomain_report(pred, u, p, t) -> Vector{String}
+
+Explain what an `isoutofdomain` predicate does at the state `(u, p, t)`,
+for a diagnostic message. Returns the predicate's source text followed by one line per
+atomic sub-expression giving its value and operand values.
+"""
+function isoutofdomain_report(tp::TracedPredicate, u, p, t)
+    lines = ["predicate: " * tp.src]
+    isempty(tp.leaves) && return lines
+    width = maximum(length(leaf.src) for leaf in tp.leaves)
+    for (idx, leaf) in enumerate(tp.leaves)
+        push!(lines, explain_leaf(idx, leaf, u, p, t, width))
+    end
+    return lines
+end
+
+isoutofdomain_report(::Any, ::Any, ::Any, ::Any) = String[]
diff --git a/lib/DiffEqBase/src/solve.jl b/lib/DiffEqBase/src/solve.jl
index 471ffc8322d..a9c680ae856 100644
--- a/lib/DiffEqBase/src/solve.jl
+++ b/lib/DiffEqBase/src/solve.jl
@@ -39,6 +39,10 @@ DiffEq's standard merging rules:
 - Passed kwargs take precedence (i.e., they override problem kwargs)
 - If `merge_callbacks=true` and both prob and kwargs have callbacks, they are
   merged into a `CallbackSet` rather than one overriding the other
+- If both prob and kwargs supply an `isoutofdomain` predicate, they are combined
+  with a logical OR. This lets a problem-carried predicate
+  (e.g. one compiled from ModelingToolkit variable bounds) coexist with a
+  user-supplied `isoutofdomain` passed to `solve`.
 
 Returns the merged kwargs as a Base.pairs.
 
@@ -67,6 +71,22 @@ function merge_problem_kwargs(prob; merge_callbacks = true, kwargs...)
             )
             kwargs = merge(kwargs_temp, callbacks)
         end
+        # Compose isoutofdomain predicates with OR instead of overwriting
+        #used when MTK hooks var bounds into isoutofdomain
+        if haskey(prob.kwargs, :isoutofdomain) && haskey(kwargs, :isoutofdomain)
+            prob_iood = prob.kwargs[:isoutofdomain]
+            user_iood = values(kwargs).isoutofdomain
+            kwargs_temp = NamedTuple{
+                Base.diff_names(
+                    Base._nt_names(values(kwargs)),
+                    (:isoutofdomain,)
+                ),
+            }(values(kwargs))
+            composed = NamedTuple{(:isoutofdomain,)}(
+                ((u, p, t) -> prob_iood(u, p, t) || user_iood(u, p, t),)
+            )
+            kwargs = merge(kwargs_temp, composed)
+        end
         kwargs = isempty(prob.kwargs) ? kwargs : merge(values(prob.kwargs), kwargs)
     end
 
@@ -532,6 +552,14 @@ explanations of the timestepping algorithms, see the
   for more details.
 * `isoutofdomain`: Specifies a function `isoutofdomain(u,p,t)` where, when it
   returns true, it will reject the timestep. Disabled by default.
+* `domain_checks`: Specifies a vector of predicates, each `(u,p,t) -> Bool` or
+  `(u,p,t) -> Pair{Bool,<:AbstractString}` (to attach a custom message), checked
+  against the state immediately before every right-hand-side evaluation (including
+  intermediate stages and Newton iterations). Unlike `isoutofdomain`, which only
+  inspects the state once per tentative step, a failing `domain_checks` predicate
+  throws before the right-hand-side is ever called on the invalid state, and is
+  caught and converted into a normal step retry with a smaller `dt` (the same
+  recovery path as a failed nonlinear solve). Disabled by default.
 * `unstable_check`: Specifies a function `unstable_check(dt,u,p,t)` where, when
   it returns true, it will cause the solver to exit and throw a warning. Defaults
   to `any(isnan,u)`, i.e. checking if any value is a NaN.
diff --git a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
index 777f6c3a678..891a221ab04 100644
--- a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
+++ b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
@@ -43,7 +43,12 @@ export AutoDePSpecialize
 import DiffEqBase: ODE_DEFAULT_NORM,
     ODE_DEFAULT_ISOUTOFDOMAIN, ODE_DEFAULT_PROG_MESSAGE,
     ODE_DEFAULT_UNSTABLE_CHECK,
-    DEVerbosity, _process_verbose_param
+    DEVerbosity, _process_verbose_param,
+    DomainCheckedFunction, domain_checks_failing,
+    isoutofdomain_report
+
+import DiffEqBase: @isoutofdomain
+export @isoutofdomain
 
 import SciMLOperators: MatrixOperator, FunctionOperator,
     update_coefficients, update_coefficients!,
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
index bd276348be0..defd0d91e04 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
@@ -600,7 +600,7 @@ function _loopfooter!(integrator)
     elseif integrator.opts.adaptive
         q = stepsize_controller!(integrator, integrator.alg)
         apply_solve_step_limiter!(integrator, ttmp)
-        integrator.isout = integrator.opts.isoutofdomain(integrator.u, integrator.p, ttmp)
+        integrator.isout = integrator.opts.isoutofdomain(integrator.u, integrator.p, ttmp) || isnan(get_EEst(integrator)) || isinf(get_EEst(integrator))
         integrator.accept_step = (
             !integrator.isout &&
                 accept_step_controller(
@@ -821,6 +821,8 @@ function SciMLBase.log_numerical_instability(integrator::ODEIntegrator; jacobian
     # each analysis gets its own section
     state_analysis = String[]
     jacobian_analysis = String[]
+    domain_analysis = String[]
+    iood_analysis = String[]
     error_analysis = String[]
 
     # state diagnostics message
@@ -890,6 +892,22 @@ function SciMLBase.log_numerical_instability(integrator::ODEIntegrator; jacobian
         end
     end
 
+    # domain checks diagnostics, find out which checks failed and report them
+    if hasproperty(integrator.f, :f) && integrator.f.f isa DomainCheckedFunction
+        failing = domain_checks_failing(integrator.f.f.pre_checks, u, integrator.p, integrator.t)
+        if failing !== nothing
+            push!(domain_analysis, "$(length(failing)) domain_checks predicate(s) failed:")
+            append!(domain_analysis, failing)
+        end
+    end
+
+    # isoutofdomain diagnostics, the other way a step is rejected for leaving a valid domain.
+    if integrator.opts.adaptive && integrator.opts.isoutofdomain(u, integrator.p, integrator.t)
+        push!(iood_analysis, "isoutofdomain predicate returned true for the proposed state, causing failure of the step")
+        # empty unless the predicate was built with `@isoutofdomain`
+        append!(iood_analysis, isoutofdomain_report(integrator.opts.isoutofdomain, u, integrator.p, integrator.t))
+    end
+
     # error estimate analysis
     if integrator.opts.adaptive
         push!(error_analysis, "step error estimate EEst = $(@sprintf("%.4g", get_EEst(integrator))) (a step is accepted when EEst <= 1)")
@@ -920,6 +938,8 @@ function SciMLBase.log_numerical_instability(integrator::ODEIntegrator; jacobian
     sections = (
         ("State Analysis", state_analysis),
         ("Jacobian Analysis", jacobian_analysis),
+        ("Domain Check Analysis: Gates inputs to function evaluations", domain_analysis),
+        ("Is Out of Domain: Gates results of integrator steps", iood_analysis),
         ("Error Analysis", error_analysis)
     )
     all(isempty(msgs) for (_, msgs) in sections) && return ""
diff --git a/lib/OrdinaryDiffEqCore/src/solve.jl b/lib/OrdinaryDiffEqCore/src/solve.jl
index 0de89bc5adc..8c42d1cb8e1 100644
--- a/lib/OrdinaryDiffEqCore/src/solve.jl
+++ b/lib/OrdinaryDiffEqCore/src/solve.jl
@@ -125,6 +125,26 @@ Internal implementation of `__init` for ODE/DAE/SDE/RODE problems. This is
 separated from `__init` so that SDE packages can call it directly, bypassing
 method dispatch (which would otherwise re-enter SDE's more specific `__init`).
 """
+# domain_checks only wraps a plain callable `f`; algorithms/functions that call
+# multiple sub-functions or an operator directly (rather than through a single
+# wrapped `f`) would silently miss the checks, so we reject them explicitly.
+function _validate_domain_checks_compatible(prob, alg)
+    f = prob.f
+    if f isa SplitFunction || f isa DynamicalODEFunction
+        msg = "domain_checks is not supported for $(nameof(typeof(f))) right-hand-sides " *
+            "(e.g. IMEX/symplectic partitioned problems), since these call multiple " *
+            "sub-functions directly rather than through a single wrapped callable."
+        throw(ArgumentError(msg))
+    end
+    if alg isa ExponentialAlgorithm
+        msg = "domain_checks is not supported for $(nameof(typeof(alg))), which accesses " *
+            "the operator/right-hand-side directly rather than calling through the " *
+            "wrapped function."
+        throw(ArgumentError(msg))
+    end
+    return nothing
+end
+
 Base.@constprop :aggressive function _ode_init(
         prob,
         alg,
@@ -164,6 +184,7 @@ Base.@constprop :aggressive function _ode_init(
         stage_limiter = trivial_limiter!,
         step_limiter = trivial_limiter!,
         isoutofdomain = ODE_DEFAULT_ISOUTOFDOMAIN,
+        domain_checks = nothing,
         unstable_check = ODE_DEFAULT_UNSTABLE_CHECK,
         verbose = Standard(),
         timeseries_errors = true,
@@ -688,6 +709,14 @@ Base.@constprop :aggressive function _ode_init(
     end
     sol = SciMLBase.build_solution(prob, _alg, ts, timeseries; _sol_kwargs...)
 
+    # wrap `f` for domain_checks after `id`/`sol` are built above 
+    # so that the solution can be used in the domain check function (e.g. for interpolation)
+    if domain_checks !== nothing && !isempty(domain_checks)
+        _validate_domain_checks_compatible(prob, alg)
+        inner_f = DomainCheckedFunction{isinplace(prob)}(f.f, domain_checks)
+        @reset f.f = inner_f
+    end
+
     FType = typeof(f)
     SolType = typeof(sol)
     cacheType = typeof(cache)
diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl
index dd2832ec595..fbdea0c1f73 100644
--- a/src/OrdinaryDiffEq.jl
+++ b/src/OrdinaryDiffEq.jl
@@ -32,7 +32,7 @@ using SciMLLogging: SciMLLogging
 using ADTypes: ADTypes, AutoForwardDiff, AutoFiniteDiff, AutoSparse
 
 # Import from OrdinaryDiffEqCore
-using OrdinaryDiffEqCore: OrdinaryDiffEqCore
+using OrdinaryDiffEqCore: OrdinaryDiffEqCore, @isoutofdomain
 
 # Import from OrdinaryDiffEqDefault
 using OrdinaryDiffEqDefault: DefaultODEAlgorithm
@@ -60,6 +60,8 @@ export SciMLBase, SciMLLogging, remake, successful_retcode, reinit!, set_propose
 # Specialization levels
 export AutoDePSpecialize
 
+export @isoutofdomain
+
 # ADTypes
 export AutoForwardDiff, AutoFiniteDiff, AutoSparse
 

From 496db5b6a5dbe04b3bf474aebc2e31aad1d4fdde Mon Sep 17 00:00:00 2001
From: Shreyas-Ekanathan <shreyase39@gmail.com>
Date: Wed, 5 Aug 2026 14:25:40 -0400
Subject: [PATCH 2/2] refactor domain checks

---
 lib/DiffEqBase/src/DiffEqBase.jl              |   3 +-
 lib/DiffEqBase/src/domain_checks.jl           | 126 +++++++++++++++++-
 lib/DiffEqBase/src/solve.jl                   |  28 ++--
 .../src/OrdinaryDiffEqCore.jl                 |   2 +-
 lib/OrdinaryDiffEqCore/src/alg_utils.jl       |  21 +++
 lib/OrdinaryDiffEqCore/src/initdt.jl          |   4 +-
 .../src/integrators/integrator_utils.jl       |   7 +-
 lib/OrdinaryDiffEqCore/src/solve.jl           |  34 +----
 .../src/OrdinaryDiffEqDifferentiation.jl      |   2 +-
 .../src/derivative_utils.jl                   |   9 +-
 .../src/OrdinaryDiffEqNonlinearSolve.jl       |   2 +-
 .../src/newton.jl                             |   2 +-
 12 files changed, 187 insertions(+), 53 deletions(-)

diff --git a/lib/DiffEqBase/src/DiffEqBase.jl b/lib/DiffEqBase/src/DiffEqBase.jl
index 086cfe60f0e..554bbce2eaf 100644
--- a/lib/DiffEqBase/src/DiffEqBase.jl
+++ b/lib/DiffEqBase/src/DiffEqBase.jl
@@ -226,7 +226,8 @@ export @isoutofdomain
             # Solver-author wrapper/tag types and convergence-testing entry type
             :EvalFunc, :OrdinaryDiffEqTag, :ConvergenceSetup,
             # Domain-violation mechanism (`domain_checks` solver keyword)
-            :DomainCheckedFunction, :domain_checks_failing,
+            :DomainCheckedFunction, :domain_checks_failing, :apply_domain_checks,
+            :strip_domain_checks, :find_domain_checks, :supports_domain_checks,
             # Source-retaining `isoutofdomain` predicates, for failure diagnostics
             :TracedPredicate, :TracedPredicateLeaf, :isoutofdomain_report
         )
diff --git a/lib/DiffEqBase/src/domain_checks.jl b/lib/DiffEqBase/src/domain_checks.jl
index 61a5d1ac134..ab73caa74ba 100644
--- a/lib/DiffEqBase/src/domain_checks.jl
+++ b/lib/DiffEqBase/src/domain_checks.jl
@@ -32,7 +32,10 @@ falls into the NaN-detection in the adaptive step controller
 (`integrator.isout = ... || isnan(EEst) || isinf(EEst)`) straight into a
 shrink-and-retry.
 
-Constructed automatically from the `domain_checks` solver keyword.
+Constructed automatically from the `domain_checks` solver keyword, by
+[`apply_domain_checks`](@ref), as the outermost layer of `f.f` — outside any
+`FunctionWrappersWrapper` `promote_f` may have installed. [`strip_domain_checks`](@ref)
+removes it again for the code paths that must not see the checks.
 """
 struct DomainCheckedFunction{iip, F, C}
     f::F
@@ -43,6 +46,8 @@ function DomainCheckedFunction{iip}(f::F, pre_checks::C) where {iip, F, C}
     return DomainCheckedFunction{iip, F, C}(f, pre_checks)
 end
 
+unwrapped_f(dcf::DomainCheckedFunction) = unwrapped_f(dcf.f)
+
 function (dcf::DomainCheckedFunction{true})(du, u, p, t)
     if domain_checks_failing(dcf.pre_checks, u, p, t) === nothing
         dcf.f(du, u, p, t)
@@ -58,6 +63,125 @@ function (dcf::DomainCheckedFunction{false})(u, p, t)
     return u .* (zero(eltype(u)) / zero(eltype(u)))
 end
 
+"""
+    strip_domain_checks(f)
+
+Remove a [`DomainCheckedFunction`](@ref) layer, returning `f` unchanged when there is none.
+Accepts either the inner callable or the `AbstractSciMLFunction` wrapping it.
+
+Checks gate *states* the integrator might accept — RK stages, Newton residuals — which all
+go through `integrator.f` and keep them. They must not gate *derivative probes*: the point
+`uprev ± ε eᵢ` a finite-difference Jacobian evaluates is not a state, `uprev` itself is
+in-domain, and only the stencil crossed the boundary. Poisoning that Jacobian column would
+reject a valid step unrecoverably, since `ε ≈ sqrt(eps) * |u|` does not scale with `dt`, so
+every retry probes the same point and the solve marches to `dtmin`, and only under
+`AutoFiniteDiff`, since `AutoForwardDiff` perturbs the dual part and leaves every predicate
+comparison in-domain.
+
+So the differentiation wrappers (`UJacobianWrapper`/`UDerivativeWrapper`/`TimeGradientWrapper`
+and their `jac_config`s) are built from a stripped `f`, and `islinearfunction` strips before
+asking `islinear`, since it hands `f.f` to the solver *as* the Jacobian.
+"""
+strip_domain_checks(dcf::DomainCheckedFunction) = dcf.f
+strip_domain_checks(f) = f
+function strip_domain_checks(f::SciMLBase.AbstractSciMLFunction)
+    hasfield(typeof(f), :f) || return f
+    inner = strip_domain_checks(f.f)
+    inner === f.f && return f
+    return @set f.f = inner
+end
+
+"""
+    find_domain_checks(f) -> checks or nothing
+
+Recover the `domain_checks` predicates carried by `f`, or `nothing` if it carries none.
+Used by failure diagnostics, which have only the integrator to work from. Kept separate
+from `f.f isa DomainCheckedFunction` so that the diagnostics do not silently go quiet if
+the layering around `f.f` ever changes.
+"""
+find_domain_checks(dcf::DomainCheckedFunction) = dcf.pre_checks
+find_domain_checks(::Any) = nothing
+function find_domain_checks(f::SciMLBase.AbstractSciMLFunction)
+    # See `strip_domain_checks`: not every `AbstractSciMLFunction` has a single `f`.
+    hasfield(typeof(f), :f) || return nothing
+    return find_domain_checks(f.f)
+end
+
+"""
+    supports_domain_checks(alg) -> Bool
+
+Whether `alg` implements the step-rejection recovery that the `domain_checks` keyword
+relies on. `false` by default: DiffEqBase installs the checks during problem
+concretization, which every solver package routes through, but only the solvers that turn
+a NaN right-hand-side into a smaller-`dt` retry can honor them. Passing `domain_checks` to
+any other solver is an error rather than a silently different meaning.
+"""
+supports_domain_checks(alg) = false
+
+"""
+    apply_domain_checks(f, domain_checks, alg; opaque_params = false) -> f
+
+Install the `domain_checks` predicates on `f` as the outermost layer of `f.f`, or return
+`f` unchanged when no checks were requested. Called from `get_concrete_problem`, i.e.
+before the cache, the nonlinear solver and the solution object are built, so that every
+consumer of `f` sees one consistent function.
+
+Rejects the combinations that cannot honor the checks instead of silently dropping them:
+right-hand-sides whose call signature is not `(du, u, p, t)`/`(u, p, t)` (split, DAE and
+second-order/partitioned forms, which would `MethodError`, or evaluate only one of several
+sub-functions), operator right-hand-sides (where the solver uses `f.f` itself as the
+Jacobian), and algorithms without the retry path.
+"""
+function apply_domain_checks(f, domain_checks, alg; opaque_params::Bool = false)
+    (domain_checks === nothing || isempty(domain_checks)) && return f
+
+    if !supports_domain_checks(alg)
+        throw(
+            ArgumentError(
+                "domain_checks is not supported for $(alg === nothing ? "this solver" : nameof(typeof(alg))). " *
+                    "The keyword relies on a failing predicate being converted into a step retry with a " *
+                    "smaller dt, which only the OrdinaryDiffEq.jl solvers implement."
+            )
+        )
+    end
+    if !(f isa ODEFunction)
+        throw(
+            ArgumentError(
+                "domain_checks is only supported for ODEFunction right-hand-sides, got " *
+                    "$(nameof(typeof(f))). Split (IMEX), DAE, and second-order/partitioned " *
+                    "right-hand-sides either take a different call signature or dispatch to several " *
+                    "sub-functions, so a single wrapped callable would miss evaluations rather than " *
+                    "gate them."
+            )
+        )
+    end
+
+    if f.f isa SciMLBase.AbstractSciMLOperator || islinear(f.f)
+        throw(
+            ArgumentError(
+                "domain_checks is not supported for operator or linear right-hand-sides. For a " *
+                    "linear `f.f` the solvers use it directly as the Jacobian (see " *
+                    "`islinearfunction`), which wrapping would silently replace with a " *
+                    "numerically-formed dense Jacobian."
+            )
+        )
+    end
+    if opaque_params
+        throw(
+            ArgumentError(
+                "domain_checks is not supported together with AutoDePSpecialize's opaque parameter " *
+                    "packing: the predicates sit outside the function wrapper that unpacks `p`, so they " *
+                    "would receive the opaque container rather than the problem's parameters."
+            )
+        )
+    end
+    # Idempotent: `get_concrete_problem` can run more than once on the same problem (e.g. a
+    # sensitivity adjoint re-concretizing), and the checks must not stack up.
+    find_domain_checks(f) === nothing || return f
+
+    return @set f.f = DomainCheckedFunction{isinplace(f)}(f.f, domain_checks)
+end
+
 
 #part 2: deal with isoutofdomain for logging purposes
 """
diff --git a/lib/DiffEqBase/src/solve.jl b/lib/DiffEqBase/src/solve.jl
index a9c680ae856..8d128be596a 100644
--- a/lib/DiffEqBase/src/solve.jl
+++ b/lib/DiffEqBase/src/solve.jl
@@ -552,14 +552,16 @@ explanations of the timestepping algorithms, see the
   for more details.
 * `isoutofdomain`: Specifies a function `isoutofdomain(u,p,t)` where, when it
   returns true, it will reject the timestep. Disabled by default.
-* `domain_checks`: Specifies a vector of predicates, each `(u,p,t) -> Bool` or
-  `(u,p,t) -> Pair{Bool,<:AbstractString}` (to attach a custom message), checked
-  against the state immediately before every right-hand-side evaluation (including
-  intermediate stages and Newton iterations). Unlike `isoutofdomain`, which only
-  inspects the state once per tentative step, a failing `domain_checks` predicate
-  throws before the right-hand-side is ever called on the invalid state, and is
-  caught and converted into a normal step retry with a smaller `dt` (the same
-  recovery path as a failed nonlinear solve). Disabled by default.
+* `domain_checks`: Specifies a tuple of predicates, each `(u,p,t) -> Bool` or
+  `(u,p,t) -> Pair{Bool,<:AbstractString}` (to attach a custom message). Each is
+  checked against the state before every right-hand-side evaluation, including
+  intermediate stages and Newton iterations, where `isoutofdomain` only inspects the
+  state once per tentative step. A failing predicate rejects the step and retries with
+  a smaller `dt`, so it requires `adaptive = true`; under fixed steps the solve stops
+  with `ReturnCode.Unstable` instead. Requires an `ODEFunction` right-hand-side and an
+  OrdinaryDiffEq.jl solver — other combinations throw rather than silently ignore the
+  checks. A `Vector` of mixed predicate types works but dispatches dynamically on
+  every evaluation. Disabled by default.
 * `unstable_check`: Specifies a function `unstable_check(dt,u,p,t)` where, when
   it returns true, it will cause the solver to exit and throw a warning. Defaults
   to `any(isnan,u)`, i.e. checking if any value is a NaN.
@@ -737,6 +739,11 @@ function get_concrete_problem(prob, isadapt; alg = nothing, kwargs...)
         tspan_promote[1], Val(_uses_forwarddiff(alg)),
         _forwarddiff_chunksize(alg)
     )
+    #build domain checks early for coherent function wrapped state
+    f_promote = apply_domain_checks(
+        f_promote, get(kwargs, :domain_checks, nothing), alg;
+        opaque_params = p_promote !== p
+    )
     if isconcreteu0(prob, tspan[1], kwargs) && prob.u0 === u0 &&
             typeof(u0_promote) === typeof(prob.u0) &&
             prob.tspan == tspan && typeof(prob.tspan) === typeof(tspan_promote) &&
@@ -770,6 +777,9 @@ function get_concrete_problem(prob::DAEProblem, isadapt; alg = nothing, kwargs..
         tspan_promote[1], Val(_uses_forwarddiff(alg)),
         _forwarddiff_chunksize(alg)
     )
+
+    f_promote = apply_domain_checks(f_promote, get(kwargs, :domain_checks, nothing), alg;opaque_params = p_promote !== p)
+
     if isconcreteu0(prob, tspan[1], kwargs) && typeof(u0_promote) === typeof(prob.u0) &&
             isconcretedu0(prob, tspan[1], kwargs) && typeof(du0_promote) === typeof(prob.du0) &&
             prob.tspan == tspan && typeof(prob.tspan) === typeof(tspan_promote) &&
@@ -802,6 +812,8 @@ function get_concrete_problem(prob::DDEProblem, isadapt; kwargs...)
     u0 = promote_u0(u0, p, tspan[1])
     tspan = promote_tspan(u0, p, tspan, prob, kwargs)
 
+    apply_domain_checks(prob.f, get(kwargs, :domain_checks, nothing), get(kwargs, :alg, nothing))
+
     return remake(prob; u0 = u0, tspan = tspan, p = p, constant_lags = constant_lags)
 end
 
diff --git a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
index 891a221ab04..40e53821f56 100644
--- a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
+++ b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
@@ -44,7 +44,7 @@ import DiffEqBase: ODE_DEFAULT_NORM,
     ODE_DEFAULT_ISOUTOFDOMAIN, ODE_DEFAULT_PROG_MESSAGE,
     ODE_DEFAULT_UNSTABLE_CHECK,
     DEVerbosity, _process_verbose_param,
-    DomainCheckedFunction, domain_checks_failing,
+    domain_checks_failing, find_domain_checks,
     isoutofdomain_report
 
 import DiffEqBase: @isoutofdomain
diff --git a/lib/OrdinaryDiffEqCore/src/alg_utils.jl b/lib/OrdinaryDiffEqCore/src/alg_utils.jl
index d0bb5f3b358..4c46ffbcd84 100644
--- a/lib/OrdinaryDiffEqCore/src/alg_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/alg_utils.jl
@@ -279,6 +279,27 @@ isimplicit(alg::OrdinaryDiffEqAdaptiveImplicitAlgorithm) = true
 isimplicit(alg::OrdinaryDiffEqImplicitAlgorithm) = true
 isimplicit(alg::CompositeAlgorithm) = any(isimplicit.(alg.algs))
 
+"""
+    supports_domain_checks(alg) -> Bool
+
+Opt in to the `domain_checks` solver keyword: these algorithms turn a `NaN` right-hand-side
+into a step rejection and a retry with a smaller `dt`, which is the recovery path a failing
+predicate relies on.
+
+Exponential and linear-exponential methods are excluded. They reach into the operator or
+the right-hand-side directly rather than calling through `integrator.f`, so a wrapped
+callable would be bypassed and the checks would silently do nothing.
+"""
+DiffEqBase.supports_domain_checks(::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = true
+function DiffEqBase.supports_domain_checks(
+        ::Union{ExponentialAlgorithm, OrdinaryDiffEqLinearExponentialAlgorithm}
+    )
+    return false
+end
+function DiffEqBase.supports_domain_checks(alg::CompositeAlgorithm)
+    return all(DiffEqBase.supports_domain_checks, alg.algs)
+end
+
 """
     isdtchangeable(alg) -> Bool
 
diff --git a/lib/OrdinaryDiffEqCore/src/initdt.jl b/lib/OrdinaryDiffEqCore/src/initdt.jl
index d4a0abdb4ec..9c88aa19b28 100644
--- a/lib/OrdinaryDiffEqCore/src/initdt.jl
+++ b/lib/OrdinaryDiffEqCore/src/initdt.jl
@@ -19,7 +19,7 @@ end
         prob, g, noise_prototype, order, integrator
     )
     _tType = eltype(t)
-    f = prob.f
+    f = DiffEqBase.strip_domain_checks(prob.f)
     p = integrator.p
     oneunit_tType = oneunit(t)
     dtmax_tdir = tdir * dtmax
@@ -348,7 +348,7 @@ end
         prob, g, order, integrator
     )
     _tType = eltype(t)
-    f = prob.f
+    f = DiffEqBase.strip_domain_checks(prob.f)
     p = prob.p
     oneunit_tType = oneunit(t)
     dtmax_tdir = tdir * dtmax
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
index defd0d91e04..d2df6590dcd 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
@@ -892,9 +892,10 @@ function SciMLBase.log_numerical_instability(integrator::ODEIntegrator; jacobian
         end
     end
 
-    # domain checks diagnostics, find out which checks failed and report them
-    if hasproperty(integrator.f, :f) && integrator.f.f isa DomainCheckedFunction
-        failing = domain_checks_failing(integrator.f.f.pre_checks, u, integrator.p, integrator.t)
+    #domain checks diagnostics
+    domain_checks = find_domain_checks(integrator.f)
+    if domain_checks !== nothing
+        failing = domain_checks_failing(domain_checks, u, integrator.p, integrator.t)
         if failing !== nothing
             push!(domain_analysis, "$(length(failing)) domain_checks predicate(s) failed:")
             append!(domain_analysis, failing)
diff --git a/lib/OrdinaryDiffEqCore/src/solve.jl b/lib/OrdinaryDiffEqCore/src/solve.jl
index 8c42d1cb8e1..7de3252570f 100644
--- a/lib/OrdinaryDiffEqCore/src/solve.jl
+++ b/lib/OrdinaryDiffEqCore/src/solve.jl
@@ -125,26 +125,6 @@ Internal implementation of `__init` for ODE/DAE/SDE/RODE problems. This is
 separated from `__init` so that SDE packages can call it directly, bypassing
 method dispatch (which would otherwise re-enter SDE's more specific `__init`).
 """
-# domain_checks only wraps a plain callable `f`; algorithms/functions that call
-# multiple sub-functions or an operator directly (rather than through a single
-# wrapped `f`) would silently miss the checks, so we reject them explicitly.
-function _validate_domain_checks_compatible(prob, alg)
-    f = prob.f
-    if f isa SplitFunction || f isa DynamicalODEFunction
-        msg = "domain_checks is not supported for $(nameof(typeof(f))) right-hand-sides " *
-            "(e.g. IMEX/symplectic partitioned problems), since these call multiple " *
-            "sub-functions directly rather than through a single wrapped callable."
-        throw(ArgumentError(msg))
-    end
-    if alg isa ExponentialAlgorithm
-        msg = "domain_checks is not supported for $(nameof(typeof(alg))), which accesses " *
-            "the operator/right-hand-side directly rather than calling through the " *
-            "wrapped function."
-        throw(ArgumentError(msg))
-    end
-    return nothing
-end
-
 Base.@constprop :aggressive function _ode_init(
         prob,
         alg,
@@ -584,18 +564,20 @@ Base.@constprop :aggressive function _ode_init(
     else
         dt
     end
+    #build cache state from stripped f 
+    cache_f = DiffEqBase.strip_domain_checks(f)
     if _cache !== nothing
         cache = _cache
     elseif prob isa DAEProblem
         cache = alg_cache(
             _alg, du, u, res_prototype, rate_prototype, uEltypeNoUnits,
-            uBottomEltypeNoUnits, tTypeNoUnits, uprev, uprev2, f, t, _dt,
+            uBottomEltypeNoUnits, tTypeNoUnits, uprev, uprev2, cache_f, t, _dt,
             reltol_internal, p, calck, Val(isinplace(prob)), verbose_spec
         )
     else
         cache = alg_cache(
             _alg, u, rate_prototype, uEltypeNoUnits, uBottomEltypeNoUnits,
-            tTypeNoUnits, uprev, uprev2, f, t, _dt, reltol_internal, p, calck,
+            tTypeNoUnits, uprev, uprev2, cache_f, t, _dt, reltol_internal, p, calck,
             Val(isinplace(prob)), verbose_spec
         )
     end
@@ -709,14 +691,6 @@ Base.@constprop :aggressive function _ode_init(
     end
     sol = SciMLBase.build_solution(prob, _alg, ts, timeseries; _sol_kwargs...)
 
-    # wrap `f` for domain_checks after `id`/`sol` are built above 
-    # so that the solution can be used in the domain check function (e.g. for interpolation)
-    if domain_checks !== nothing && !isempty(domain_checks)
-        _validate_domain_checks_compatible(prob, alg)
-        inner_f = DomainCheckedFunction{isinplace(prob)}(f.f, domain_checks)
-        @reset f.f = inner_f
-    end
-
     FType = typeof(f)
     SolType = typeof(sol)
     cacheType = typeof(cache)
diff --git a/lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl b/lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl
index 8f2155aafe9..3f627ac5e84 100644
--- a/lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl
+++ b/lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl
@@ -51,7 +51,7 @@ using FastBroadcast: @..
 
 using ConcreteStructs: @concrete
 
-import DiffEqBase: OrdinaryDiffEqTag
+import DiffEqBase: OrdinaryDiffEqTag, strip_domain_checks
 
 # Functions for sparse array handling - will be overloaded by extension
 # Default implementations return false/error for non-sparse types
diff --git a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl
index ccb849a5112..7961ebf46cd 100644
--- a/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl
+++ b/lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl
@@ -339,7 +339,7 @@ function calc_J(integrator, cache, next_step::Bool = false)
         else
             (; uf) = cache
 
-            uf.f = nlsolve_f(f, alg)
+            uf.f = strip_domain_checks(nlsolve_f(f, alg)) #use non-domain checks formulation
             uf.p = p
             uf.t = t
             J = jacobian(uf, uprev, integrator)
@@ -402,7 +402,7 @@ function calc_J!(J, integrator, cache, next_step::Bool = false)
             f.jac(J, uprev, p, t)
         else
             (; du1, uf, jac_config) = cache
-            uf.f = nlsolve_f(f, alg)
+            uf.f = strip_domain_checks(nlsolve_f(f, alg))
             uf.t = t
             if !(p isa SciMLBase.NullParameters)
                 uf.p = p
@@ -513,8 +513,9 @@ islinearfunction(integrator) = islinearfunction(integrator.f, integrator.alg)
 return the tuple `(is_linear_wrt_odealg, islinearodefunction)`.
 """
 function islinearfunction(f::F, alg)::Tuple{Bool, Bool} where {F}
-    isode = f isa ODEFunction && islinear(f.f)
-    islin = isode || (issplit(alg) && f isa SplitFunction && islinear(f.f1.f))
+    isode = f isa ODEFunction && islinear(strip_domain_checks(f.f)) #drop domain checks
+    islin = isode ||
+        (issplit(alg) && f isa SplitFunction && islinear(strip_domain_checks(f.f1.f)))
     return islin, isode
 end
 
diff --git a/lib/OrdinaryDiffEqNonlinearSolve/src/OrdinaryDiffEqNonlinearSolve.jl b/lib/OrdinaryDiffEqNonlinearSolve/src/OrdinaryDiffEqNonlinearSolve.jl
index 332c7492bb3..d1bdf296989 100644
--- a/lib/OrdinaryDiffEqNonlinearSolve/src/OrdinaryDiffEqNonlinearSolve.jl
+++ b/lib/OrdinaryDiffEqNonlinearSolve/src/OrdinaryDiffEqNonlinearSolve.jl
@@ -12,7 +12,7 @@ using SciMLBase: DAEFunction, DEIntegrator, NonlinearFunction, NonlinearProblem,
     _vec, _reshape, postamble!, alg_order, isadaptive
 import DiffEqBase
 import DiffEqBase: OrdinaryDiffEqTag, calculate_residuals, calculate_residuals!,
-    BrownFullBasicInit, ShampineCollocationInit
+    BrownFullBasicInit, ShampineCollocationInit, strip_domain_checks
 import ConstructionBase
 import PreallocationTools: DiffCache, get_tmp
 using SimpleNonlinearSolve: SimpleTrustRegion, SimpleGaussNewton
diff --git a/lib/OrdinaryDiffEqNonlinearSolve/src/newton.jl b/lib/OrdinaryDiffEqNonlinearSolve/src/newton.jl
index c1093f07ab9..ce951dba86e 100644
--- a/lib/OrdinaryDiffEqNonlinearSolve/src/newton.jl
+++ b/lib/OrdinaryDiffEqNonlinearSolve/src/newton.jl
@@ -186,7 +186,7 @@ function _update_nlsolvealg_W!(nlcache, integrator, dtgamma, tstep, new_jac = tr
             if SciMLBase.has_jac(f)
                 f.jac(J, uprev, p, tstep)
             elseif uf !== nothing
-                uf.f = nlsolve_f(f, alg)
+                uf.f = strip_domain_checks(nlsolve_f(f, alg)) #strip domain checks
                 uf.t = tstep
                 if !(p isa SciMLBase.NullParameters)
                     uf.p = p