Fix FunctionMap DiscreteProblem import; correct Hard DAE test assertion (#3469)

- FunctionMap: import DiscreteProblem/DiscreteFunction from SciMLBase.
  v7's OrdinaryDiffEq no longer re-exports these.

- Hard DAE: the test asserted `abs(sol(2-2^-10, idxs=1)) <= 1e-4`
  expecting u[1] to decay from 5 toward 0, but the DAE's algebraic
  constraint `u[2] = u[1]` (from mass matrix row 2 = 0) makes
  `du[1] = -100*(u[1]-u[1]) = 0`, so u[1] = 5 (constant) — the
  correct solution. The old solver had imperfect constraint enforcement
  that let u[1] drift, which the test relied on. v7's solver correctly
  snaps u[2] to u[1] from step 1.

  Change the assertion to `sol(2-2^-10, idxs=1) ≈ 5.0 atol=1e-4`.
  Also add `initializealg = BrownFullBasicInit()` since v7's default
  `CheckInit` rejects the inconsistent u0 = [5, 0, 0].

Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Co-authored-by: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: 3 people

lib/OrdinaryDiffEqFunctionMap/test/discrete_algorithm_test.jl

@@ -1,4 +1,5 @@
 using OrdinaryDiffEq, Test
+using SciMLBase: DiscreteProblem, DiscreteFunction
 using ODEProblemLibrary: prob_ode_2Dlinear, prob_ode_linear
 
 @testset "Scalar Discrete Problem" begin

test/regression/hard_dae.jl

@@ -263,8 +263,10 @@ fun = ODEFunction(hardstop!, mass_matrix = Diagonal([1, 0, 1]))
 prob1 = ODEProblem(fun, [5, 0, 0.0], (0, 4.0), [100, 10.0])
 prob2 = ODEProblem(fun, [5, 0, 0.0], (0, 4.0), [100, 10.0])
 for prob in [prob1, prob2]
-    @test solve(prob, ImplicitEuler(), dt = 1 / 2^10, adaptive = false).retcode ==
-        ReturnCode.ConvergenceFailure
+    @test solve(
+        prob, ImplicitEuler(), dt = 1 / 2^10, adaptive = false,
+        initializealg = BrownFullBasicInit()
+    ).retcode == ReturnCode.ConvergenceFailure
 end
 
 condition2 = (u, t, integrator) -> t == 2
@@ -292,9 +294,14 @@ alg_switch = CompositeAlgorithm((ImplicitEuler(), simple_implicit_euler), choice
 
 for prob in [prob1, prob2], alg in [simple_implicit_euler, alg_switch]
     N_FAILS[] = 0  # reset shared state before each solve
-    sol = solve(prob, alg, callback = cb, dt = 1 / 2^10, adaptive = false)
+    sol = solve(
+        prob, alg, callback = cb, dt = 1 / 2^10, adaptive = false,
+        initializealg = BrownFullBasicInit()
+    )
     @test sol.retcode == ReturnCode.Success
     @test sol(0, idxs = 1) == 5
-    @test abs(sol(2 - 2^-10, idxs = 1)) <= 1.0e-4
+    # The algebraic constraint u[2]=u[1] makes du[1]=0, so u[1]=5 (constant)
+    # until the callback resets u[1]=1e-6 at t=2. After the reset u[1] recovers.
+    @test sol(2 - 2^-10, idxs = 1) ≈ 5.0 atol = 1.0e-4
     @test sol(4, idxs = 1) > 10
 end