Add support for the initializealg argument in SciMLBase callbacks

[BenChung]

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

Additional context

Add any other context about the problem here.

[jmaedler]

Setting the leakage streams to zero will potentially lead to inifite solutions for the value between the valves, because you have "dirichlet" boundary values at the input of the first valve and the output of the second valve. The pressure difference between is distributed based the implicit relationship Vdot = Sf(\Delta p) + k_leakage\Delta p . Thus, there is no clear solution if all "S" and "k_leakage" are 0. I am not sure how the Julia solver will deal with this kind of equation system breakdown.

[jmaedler]

Let me suggest a proper set of test conditions until tomorrow ;-).

[jmaedler]

Ok, if I would test the behavior of the system, I would use the following checks. These numbers were calculated with NonlinearSolve assuming a k_leakage = 1e-08, furthermore I exchanged the valve equation to smooth its behavior. I'll make an extra post on that.

# Case 1: Both valves are open
@test ( (sol(45;idxs=sys.binary_valve_1.Δp) ≈ 0.5) &&
        (sol(45;idxs=sys.binary_valve_2.Δp) ≈ 0.5) &&
        (sol(45;idxs=sys.binary_valve_1.Vdot) ≈ 0.707107493292268) &&
        sol(45;idxs=sys.binary_valve_1.Vdot) == sol(45;idxs=sys.binary_valve_2.Vdot))
# Case 2: First valve is open, second valve is closed
@test ( isapprox(sol(45;idxs=sys.binary_valve_1.Δp),2.6913001348648953e-10;atol=1e-06) &&
        isapprox(sol(45;idxs=sys.binary_valve_2.Δp),0.99999999973087;atol=1e-06) &&
        isapprox(sol(45;idxs=sys.binary_valve_1.Vdot),9.9999999973087e-9;atol=1e-06) &&
        sol(45;idxs=sys.binary_valve_1.Vdot) == sol(45;idxs=sys.binary_valve_2.Vdot))
# Case 3: First valve is closed, second valve is open
@test ( isapprox(sol(45;idxs=sys.binary_valve_1.Δp),0.99999999973087;atol=1e-06) &&
        isapprox(sol(45;idxs=sys.binary_valve_2.Δp),2.6913001348648953e-10;atol=1e-06) &&
        isapprox(sol(45;idxs=sys.binary_valve_1.Vdot),9.9999999973087e-9;atol=1e-06) &&
        sol(45;idxs=sys.binary_valve_1.Vdot) == sol(45;idxs=sys.binary_valve_2.Vdot))
# Case 4: Both valves are closed
@test ( isapprox(sol(45;idxs=sys.binary_valve_1.Δp),0.5;atol=1e-06) &&
        isapprox(sol(45;idxs=sys.binary_valve_2.Δp),0.5;atol=1e-06) &&
        (sol(45;idxs=sys.binary_valve_1.Vdot) ≈ 5e-09) &&
        sol(45;idxs=sys.binary_valve_1.Vdot) == sol(45;idxs=sys.binary_valve_2.Vdot))

[jmaedler]

Smooth version of this equation

Vdot ~ S*k_V*tanh(1175*Δp/p_ref)*sqrt( sqrt(Δp^2 + 0.001^2)/p_ref * ρ_ref/ρ) + k_leakage*Δp

[jmaedler]

@BenChung We can also start with a simpler example with just one valve to avoid the numerical breakdown issue. I'll make an example - just let me know.

[jmaedler]

Hey @BenChung and @ChrisRackauckas,
as I wrote briefly in Slack I solved the issue with the discontinuous states during callbacks. So, if you want a running version with just switching but continuous state, you can use the following code. It even does not need any leakage anymore. A little compressibility between the valves does the trick :-). Makes the system more stiff but physical correct.
(...if reinitialization with state jumps would become a thing it is a trade-off between the computational cost of this stiffness the reinitialization.)

using ModelingToolkit
using ModelingToolkit: D_nounits as D, t_nounits as t
using DifferentialEquations
using Test

@connector LiquidPort begin
    p(t)::Float64, [    description = "Set pressure in bar", 
                        guess = 1.01325]
    Vdot(t)::Float64, [ description = "Volume flow rate in L/min", 
                        guess = 0.0, 
                        connect = Flow]
end

@mtkmodel PressureSource begin
    @components begin
        port = LiquidPort()
    end
    @parameters begin
        p_set::Float64 = 1.01325, [description = "Set pressure in bar"]
    end
    @equations begin
        port.p ~ p_set
    end
end

@mtkmodel BinaryValve begin
    @constants begin
        p_ref::Float64 = 1.0, [description = "Reference pressure drop in bar"]
        ρ_ref::Float64 = 1000.0, [description = "Reference density in kg/m^3"]
    end
    @components begin
        port_in = LiquidPort()
        port_out = LiquidPort()
    end
    @parameters begin 
        k_V::Float64 = 1.0, [description = "Valve coefficient in L/min/bar"]
        k_leakage::Float64 = 0.0, [description = "Leakage coefficient in L/min/bar"]
        ρ::Float64 = 1000.0, [description = "Density in kg/m^3"]
    end
    @variables begin
        S(t)::Float64, [description = "Valve state", guess = 0.0, irreducible = true]
        Δp(t)::Float64, [description = "Pressure difference in bar", guess = 1.0]
        Vdot(t)::Float64, [description = "Volume flow rate in L/min", guess = 1.0]
    end
    @equations begin
        # Port handling
        port_in.Vdot ~ -Vdot
        port_out.Vdot ~ Vdot
        Δp ~ port_in.p - port_out.p
        # System behavior
        D(S) ~ 0.0
        Vdot ~ S*k_V*tanh(1175*Δp/p_ref)*sqrt( sqrt(Δp^2 + 0.001^2)/p_ref * ρ_ref/ρ) + k_leakage*Δp
    end
end

@mtkmodel CompressibleHoldup begin
    @components begin
        port_in = LiquidPort()
        port_out = LiquidPort()
    end
    @parameters begin
        β::Float64 = 4.6e-05, [description = "Compressibility in 1/bar"]
        p_ref::Float64 = 1.01325, [description = "Reference pressure in bar"]
        V_holdup::Float64 = 1.0, [description = "Holdup volume in L"]
    end
    @variables begin
        V(t)::Float64, [description = "Volume in L", guess = 1e-03]
        p(t)::Float64, [description = "Pressure in bar", guess = 1.01325]
        Vdot_in(t)::Float64, [description = "Volume flow rate in L/min", guess = 0.0]
        Vdot_out(t)::Float64, [description = "Volume flow rate in L/min", guess = 0.0]
    end
    @equations begin
        # Port handling
        port_in.Vdot ~ -Vdot_in
        port_out.Vdot ~ Vdot_out
        port_in.p ~ p
        port_out.p ~ p
        # System behavior
        D(V) ~ Vdot_in - Vdot_out
        V_holdup ~ V*exp(-β*(p - p_ref))
    end
end

# Test System
@mtkmodel TestSystem begin
    @components begin
        pressure_source_1 = PressureSource(p_set = 2.0)
        binary_valve_1 = BinaryValve(S = 1.0, k_leakage = 0.0)
        compressible_holdup = CompressibleHoldup(V=1.0, p=1.0, β = 1.0e-01)
        binary_valve_2 = BinaryValve(S = 1.0, k_leakage = 0.0)
        pressure_source_2 = PressureSource(p_set = 1.0)
    end
    @equations begin
        connect(pressure_source_1.port, binary_valve_1.port_in)
        connect(binary_valve_1.port_out, compressible_holdup.port_in)
        connect(compressible_holdup.port_out, binary_valve_2.port_in)
        connect(binary_valve_2.port_out, pressure_source_2.port)
    end
    @discrete_events begin
        [30] => [binary_valve_1.S ~ 0.0]
        [60] => [binary_valve_1.S ~ 1.0, binary_valve_2.S ~ 0.0]
        [90] => [binary_valve_1.S ~ 0.0]
    end
end

# Test Simulation
@mtkbuild sys = TestSystem()

# Test Simulation
prob = ODEProblem(sys, [], (0.0, 120.0))
sol = solve(prob, QNDF())

# Case 1: Both valves are open
@test ( isapprox(sol(29;idxs=sys.binary_valve_1.Δp), 0.5; atol=0.01) &&
        isapprox(sol(29;idxs=sys.binary_valve_2.Δp), 0.5; atol=0.01) &&
        isapprox(sol(29;idxs=sys.binary_valve_1.Vdot), sqrt(0.5); atol=0.01) &&
        isapprox(sol(29;idxs=sys.binary_valve_1.Vdot), sol(29;idxs=sys.binary_valve_2.Vdot); atol=0.01))
# Case 2: First valve is closed, second valve is open
@test ( isapprox(sol(59;idxs=sys.binary_valve_1.Δp), 1.0; atol=0.01) &&
        isapprox(sol(59;idxs=sys.binary_valve_2.Δp), 0.0; atol=0.01) &&
        isapprox(sol(59;idxs=sys.binary_valve_1.Vdot), 0.0; atol=0.01) &&
        isapprox(sol(59;idxs=sys.binary_valve_2.Vdot), 0.0; atol=0.01))
# Case 3: First valve is open, second valve is closed
@test ( isapprox(sol(89;idxs=sys.binary_valve_1.Δp), 0.0; atol=0.01) &&
        isapprox(sol(89;idxs=sys.binary_valve_2.Δp), 1.0; atol=0.01) &&
        isapprox(sol(89;idxs=sys.binary_valve_1.Vdot), 0.0; atol=0.01) &&
        isapprox(sol(89;idxs=sys.binary_valve_2.Vdot), 0.0; atol=0.01))
# Case 4: Both valves are closed
@test ( isapprox(sol(119;idxs=sys.binary_valve_1.Δp), 0.0;atol=0.01) &&
        isapprox(sol(119;idxs=sys.binary_valve_2.Δp), 1.0;atol=0.01) &&
        isapprox(sol(119;idxs=sys.binary_valve_1.Vdot), 0.0; atol=0.01) &&
        isapprox(sol(119;idxs=sys.binary_valve_2.Vdot), 0.0; atol=0.01))

# using Plots
# p1 = plot(sol; idxs=[sys.binary_valve_1.S,
#                 sys.binary_valve_2.S]);
# p2 = plot(sol; idxs=[sys.binary_valve_1.Vdot,
#                 sys.binary_valve_2.Vdot]);     
# p3 = plot(sol; idxs=[sys.binary_valve_1.Δp, 
#                 sys.binary_valve_2.Δp]);
# plot(p1, p2, p3; layout=(3,1), size=(800, 800))