Why "Cyclic Guesses" error message when creating ODEProblem?

[B-LIE]

Question❓

A rather technical question, so perhaps better suited here?

I have a component of a pipe filled with fluid, where the fluid is compressible and the walls are elastic. I try to make a component using the "method of characteristic" for discretization of the PDEs, essentially because many in the industry use this method - and I'd like to compare results.

So the equations consist of 2 PDEs with spatial derivatives, one for mass flow rate md, and one for pressure. In addition, there is a "PDE" with variable "chi" which is independent of neighbor cells. I connect this pipe to a hydraulic port with mass flow rate as Flow, and pressure as potential.

In principle, for a discretization with N cells, there should be 3N unknowns, but 2 should be equal to port values. Possibly the pressures at each end?

---

The process of mtkcompile works. But when I try to create the ODEProblem, I get (a) a warning about "The initialization system is structurally singular. ", and (b) a crash with an error message:

Cyclic guesses detected in the system. Symbolic values were found for the following variables/parameters in the map: 
(pipe_moc₊mdˍt(t))[1]  => (pipe_moc₊mdˍt(t))[1]
(pipe_moc₊chiˍt(t))[5]  => (pipe_moc₊chiˍt(t))[5]
(pipe_moc₊chiˍt(t))[4]  => (pipe_moc₊chiˍt(t))[4]
(pipe_moc₊mdˍt(t))[7]  => (pipe_moc₊mdˍt(t))[7]
In order to resolve this, please provide additional numeric guesses so that the chain can be resolved to assign numeric values to each variable. Alternatively, the `missing_guess_value` keyword can be used to set a fallback guess for all variables. The keyword must be passed an instance of the `MissingGuessValue` sum-type.

It is not clear to me what causes the error message.

If it is any help, here is the component -- somewhat ugly coding, perhaps -- with lots of "help" from Google AI Studio...

@component function TaperedPipeMOC(; 
    name, media::HydraulicMedia, conduit::ConduitStructure, 
    L, Da, Db, z_a, z_b, 
    N = 10, eps_r = 1e-4, k_VB = 0.033, g_T = 15.0, friction_dynamic::Bool = true
)
    @named a = HydraulicPort()
    @named b = HydraulicPort()

    dx = L / (N - 1)
    g_acc = 9.80665
    mu    = media.mu_mix
    rho_0 = media.rho_mix_0
    p_0   = media.p_0

    pars = @parameters begin
        L=L; z_a=z_a; z_b=z_b; ϵ_r=eps_r; k=k_VB; g_T=g_T; f_u = friction_dynamic ? 1.0 : 0.0
    end

    # Profile Guesses for Pressure (Newton solver starting point)
    z_vec = [z_a + (z_b - z_a) * (i - 1) / (N - 1) for i in 1:N]
    z_high = max(z_a, z_b)
    p_guess_vec = [p_0 + rho_0 * g_acc * (z_high - z_vec[i]) for i in 1:N]

    # --- THE STATES ---
    # p is determined by the boundaries -> use [guess]
    @variables p(t)[1:N], [unit=u"Pa", guess=p_guess_vec]
    
    # md and chi are differential states -> use [=] to FIX initial conditions at 0.0
    # In MTK v9, providing a vector to the right-hand side is the correct way.
    @variables md(t)[1:N] = fill(0.0, N), [unit=u"kg/s"]
    @variables chi(t)[1:N] = fill(0.0, N), [unit=u"N/m"]

    # Geometric vectors
    D_vec = [Da + (Db - Da) * (i - 1) / (N - 1) for i in 1:N]
    A_vec = [pi * d^2 / 4 for d in D_vec]

    eqs = Equation[]
    
    for j in 1:N
        Dj = D_vec[j]; Aj = A_vec[j]
        
        # Physics Functors
        rho_j   = media(Val(:rho), p[j])
        beta_Tj = media(Val(:beta_T), p[j])
        beta_Cj = conduit(Dj)
        aj      = 1.0 / sqrt(rho_j * (beta_Tj + beta_Cj))
        vj      = md[j] / (rho_j * Aj)
        
        # MOC Coefficients
        Kp_j = Aj * (vj - aj) / aj^2
        Km_j = Aj * (vj + aj) / aj^2
        
        # Friction [N/m]
        N_Re = abs(vj) * Dj * rho_j / mu
        f_F  = 0.3086 / (log10(6.9/(N_Re + 1e-2) + (ϵ_r/(3.7*Dj))^1.11))^2
        F_fs_j = ifelse(N_Re < 2300, (32 * mu * md[j]) / (rho_j * Dj^2), (f_F * md[j] * abs(md[j])) / (2 * rho_j * Aj * Dj))
        
        # Relaxation
        T_z = ifelse(N_Re < 2300, (k * rho_j * Dj^2) / (8 * mu), (k * Dj) / (sqrt(vj^2 + 1e-4) * max(1e-4, f_F)))
        T_p = max(1e-3, T_z / g_T)
        
        push!(eqs, D(chi[j]) ~ ((1 - g_T) * F_fs_j - chi[j]) / T_p)
        
        F_f_j = ifelse(f_u > 0.5, g_T * F_fs_j + chi[j], F_fs_j)
        F_g_j = -rho_j * Aj * g_acc * (z_b - z_a) / L
        sigma_j = F_g_j - F_f_j

        if j == 1
            push!(eqs, p[1] ~ a.p)
            wm = -(vj - aj) * ((p[2] - p[1])/dx - (aj/Aj) * (md[2] - md[1])/dx)
            push!(eqs, D(p[1]) + (1/Km_j)*D(md[1]) ~ wm + sigma_j/Km_j)
        elseif j == N
            push!(eqs, p[N] ~ b.p)
            wp = -(vj + aj) * ((p[N] - p[N-1])/dx + (aj/Aj) * (md[N] - md[N-1])/dx)
            push!(eqs, D(p[N]) + (1/Kp_j)*D(md[N]) ~ wp + sigma_j/Kp_j)
        else
            wp = -(vj + aj) * ((p[j] - p[j-1])/dx + (aj/Aj) * (md[j] - md[j-1])/dx)
            wm = -(vj - aj) * ((p[j+1] - p[j])/dx - (aj/Aj) * (md[j+1] - md[j])/dx)
            push!(eqs, D(p[j]) ~ (wp - wm) / (Kp_j - Km_j))
            push!(eqs, D(md[j]) ~ -(Kp_j * wm - Km_j * wp) / (Kp_j - Km_j) + sigma_j)
        end
    end

    push!(eqs, a.md ~ md[1])
    push!(eqs, b.md ~ -md[N])

    # Ensure all scalar array elements are registered as states
    sts = [p; md; chi]
    
    compose(ODESystem(eqs, t, sts, pars; name=name), a, b)
end

Any ideas? Or too little information? I should add that I also have a simpler pipe version (inelastic). That one works fine with the connected components.

[oscardssmith]

That should be resolved by #4445 (if we go that direction). The TLDR is that currently for algebraic equations we require a guess for initial values as a place to start the newton rhapson iterations to find your initial conditions. With #4445, pseudorandom values will be chosen for equations where the user doesn't supply a guess which should be a fairly nice UX improvement.

[B-LIE]

Thanks for quick response. One question and one comment:

Question:
In the code I include, I have defined 3 vectors as "variables". And I have specified initial values for two of the vectors, and initial guesses for the third. That is the only variables that are passed on to ODESystem.

My problem is: which variables do I lack initial guesses for???

OK, it could be that since I also add other expressions of symbolic variables, such as:
f_F = 0.3086 / (log10(6.9/(N_Re + 1e-2) + (ϵ_r/(3.7*Dj))^1.11))^2,

that these also in reality become "variables" (should I have declared them as so) that require initial guesses???

Comment:
Automating insertion of guesses sounds like a good idea. However, in practice, using MTK to implement models of "mechanistic models" of physical systems, often there may be problems in choosing values with value zero -- often they should be positive. And for some variables (say, pressure), the values should probably rather be 1e5 than 1.

[B-LIE]

Update: I have converted every "label" of symbolic expressions in the code into @variables with guess values. So now there is no variable without guess value or initial value. Still the same error message.

[B-LIE]

New update: I reduced the problem to case N=2. Then I got a solution, but very wrong. With some rescaling to improve numerics, I get reasonable results for N=2. I need to do some more checking for N > 2.