Recurse on concrete initialization parameters

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

Author: ChrisRackauckas

lib/ModelingToolkitBase/src/systems/nonlinear/initializesystem.jl

@@ -877,8 +877,15 @@ function DiffEqBase.get_updated_symbolic_problem(
     t0 = is_time_dependent(prob) ? current_time(prob) : nothing
 
     unwrapped_p = _unwrap_mtk_parameters(p)
-    if unwrapped_p isa MTKParameters
-        buffer = unwrapped_p.initials
+    if unwrapped_p !== p
+        updated_prob = DiffEqBase.get_updated_symbolic_problem(
+            sys, prob; u0, p = unwrapped_p, kw...
+        )
+        return @set updated_prob.p = p
+    end
+
+    if p isa MTKParameters
+        buffer = p.initials
     else
         buffer = p
     end

lib/ModelingToolkitBase/test/despecialized_parameters.jl

@@ -1,12 +1,18 @@
 using ModelingToolkitBase
 using ModelingToolkitBase: t_nounits as t, D_nounits as D, MTKParameters
 using DiffEqBase
+using ForwardDiff
 using OrdinaryDiffEq
 using Random
+using SciMLBase
 using SciMLStructures
 using SymbolicIndexingInterface
 using Test
 
+struct _OpaqueParameterWrapper end
+const _opaque_parameters = Ref{Any}()
+SciMLBase.unwrap_parameters(::_OpaqueParameterWrapper) = _opaque_parameters[]
+
 @testset "AutoDespecialize parameters" begin
     @parameters a = 2.0
     @variables x(t) = 1.0 z(t)
@@ -63,6 +69,20 @@ using Test
     @test SciMLStructures.isscimlstructure(concrete_prob.p)
     @test SciMLStructures.ismutablescimlstructure(concrete_prob.p)
 
+    dual_initials = map(
+        value -> ForwardDiff.Dual{Nothing}(value, one(value)), concrete_parameters.initials
+    )
+    dual_parameters = SciMLStructures.replace(
+        SciMLStructures.Initials(), concrete_parameters, dual_initials
+    )
+    _opaque_parameters[] = SciMLBase.DespecializedParameters(dual_parameters)
+    opaque_parameters = _OpaqueParameterWrapper()
+    updated_prob = DiffEqBase.get_updated_symbolic_problem(
+        sys, concrete_prob; u0 = concrete_prob.u0, p = opaque_parameters
+    )
+    @test eltype(updated_prob.u0) <: ForwardDiff.Dual
+    @test updated_prob.p === opaque_parameters
+
     seen_parameter_type = Ref{DataType}()
     gfw = ModelingToolkitBase.GeneratedFunctionWrapper{Tuple{2, 3, false}}(
         (u, p, t) -> (seen_parameter_type[] = typeof(p); u),