fix: promote resid_prototype using tunables

src/systems/nonlinear/nonlinearsystem.jl

@@ -295,7 +295,7 @@ function SciMLBase.NonlinearFunction(sys::NonlinearSystem, args...; kwargs...)
 end
 
 function SciMLBase.NonlinearFunction{iip}(sys::NonlinearSystem, dvs = unknowns(sys),
- ps = parameters(sys), u0 = nothing;
+ ps = parameters(sys), u0 = nothing, p = nothing;
  version = nothing,
  jac = false,
  eval_expression = false,
@@ -327,11 +327,22 @@ function SciMLBase.NonlinearFunction{iip}(sys::NonlinearSystem, dvs = unknowns(s
 
  observedfun = ObservedFunctionCache(sys; eval_expression, eval_module)
 
+ if length(dvs) == length(equations(sys))
+ resid_prototype = nothing
+ elseif u0 === nothing
+ u0ElType = Float64
+ if SciMLStructures.isscimlstructure(p)
+ u0ElType = promote_type(
+ eltype(SciMLStructures.canonicalize(SciMLStructures.Tunable(), p)[1]),
+ u0ElType)
+ end
+ resid_prototype = zeros(u0ElType, length(equations(sys)))
+ end
+
  NonlinearFunction{iip}(f,
  sys = sys,
  jac = _jac === nothing ? nothing : _jac,
- resid_prototype = length(dvs) == length(equations(sys)) ? nothing :
- zeros(length(equations(sys))),
+ resid_prototype = resid_prototype,
  jac_prototype = sparse ?
  similar(calculate_jacobian(sys, sparse = sparse),
  Float64) : nothing,
@@ -355,7 +366,7 @@ variable and parameter vectors, respectively.
 struct NonlinearFunctionExpr{iip} end
 
 function NonlinearFunctionExpr{iip}(sys::NonlinearSystem, dvs = unknowns(sys),
- ps = parameters(sys), u0 = nothing;
+ ps = parameters(sys), u0 = nothing, p = nothing;
  version = nothing, tgrad = false,
  jac = false,
  linenumbers = false,
@@ -376,8 +387,18 @@ function NonlinearFunctionExpr{iip}(sys::NonlinearSystem, dvs = unknowns(sys),
  end
 
  jp_expr = sparse ? :(similar($(get_jac(sys)[]), Float64)) : :nothing
- resid_expr = length(dvs) == length(equations(sys)) ? :nothing :
- :(zeros($(length(equations(sys)))))
+ if length(dvs) == length(equations(sys))
+ resid_expr = :nothing
+ elseif u0 === nothing
+ u0ElType = Float64
+ if SciMLStructures.isscimlstructure(p)
+ u0ElType = promote_type(
+ eltype(SciMLStructures.canonicalize(SciMLStructures.Tunable(), p)[1]),
+ u0ElType)
+ end
+
+ resid_expr = :(zeros($u0ElType, $(length(equations(sys)))))
+ end
  ex = quote
  f = $f
  jac = $_jac
@@ -412,7 +433,7 @@ function process_NonlinearProblem(constructor, sys::NonlinearSystem, u0map, para
  check_eqs_u0(eqs, dvs, u0; kwargs...)
  end
 
- f = constructor(sys, dvs, ps, u0; jac = jac, checkbounds = checkbounds,
+ f = constructor(sys, dvs, ps, u0, p; jac = jac, checkbounds = checkbounds,
  linenumbers = linenumbers, parallel = parallel, simplify = simplify,
  sparse = sparse, eval_expression = eval_expression, eval_module = eval_module,
  kwargs...)

test/nonlinearsystem.jl

@@ -3,6 +3,7 @@ using ModelingToolkit: get_metadata
 using DiffEqBase, SparseArrays
 using Test
 using NonlinearSolve
+using ForwardDiff
 using ModelingToolkit: value
 using ModelingToolkit: get_default_or_guess, MTKParameters
 
@@ -325,3 +326,21 @@ end
  prob = @test_nowarn NonlinearProblem(sys, nothing)
  @test_nowarn solve(prob)
 end
+
+@testset "resid_prototype when system has no unknowns and an equation" begin
+ @variables x
+ @parameters p
+ @named sys = NonlinearSystem([x ~ 1, x^2 - p ~ 0])
+ for sys in [
+ structural_simplify(sys, fully_determined = false),
+ structural_simplify(sys, fully_determined = false, split = false)
+ ]
+ @test length(equations(sys)) == 1
+ @test length(unknowns(sys)) == 0
+ T = typeof(ForwardDiff.Dual(1.0))
+ prob = NonlinearProblem(sys, [], [p => ForwardDiff.Dual(1.0)]; check_length = false)
+ @test prob.f(Float64[], prob.p) isa Vector{T}
+ @test prob.f.resid_prototype isa Vector{T}
+ @test_nowarn solve(prob)
+ end
+end