# Solving differential equations for multi-component fluid flow

**URL:** <https://discourse.julialang.org/t/solving-differential-equations-for-multi-component-fluid-flow/116254>\
**Category:** Modelling & Simulations\
**Tags:** physics, differentialequation, fluid-flow\
**Created:** [June 26, 2024, 11:40am UTC](https://discourse.julialang.org/t/solving-differential-equations-for-multi-component-fluid-flow/116254 "2024-06-26T11:40:06Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![DanaO](https://avatars.discourse-cdn.com/v4/letter/d/b5ac83/32.png) [@DanaO](https://discourse.julialang.org/u/DanaO)\
**Post date:** [June 26, 2024, 11:40am UTC](https://discourse.julialang.org/t/solving-differential-equations-for-multi-component-fluid-flow/116254/1 "2024-06-26T11:40:06Z")

</div>

Hello everyone,

I would like to model a **multi-component** , **isothermal gas flow** in a **1D** pipeline. To do this, I would like to solve the following simplified **differential equations** :

\quad \frac{\partial p}{\partial t} = - a \frac{\partial q\_{m}}{\partial x} \\ \quad \frac{\partial p}{\partial x} = -b \frac{q\_{m}|q\_{m}|}{p} \\ \quad \frac{\partial c}{\partial t} = - a \frac{q\_m}{p}\frac{\partial c}{\partial x}

The variables are the pressure p(x,t), the mass flow q\_m(x,t) and the methane mass fraction c(x,t). a and b are constants. The **initial state** is the steady-state flow, i.e.

\quad q\_m(x,0) = q\_0 \\ \quad p(x,0) = \sqrt{p\_0^2 - 2b\, q\_0|q\_0|x} \\ \quad c(x,0) = 1 \,.

For the **boundary conditions** I choose the outlet mass flow, the inlet pressure and the methane fraction,

\quad p(0,t) = p\_{in}(t) \ \quad c(0,t) = c\_{in}(t) \ \quad q\_m(L,t) = q\_{out}(t)\, .  
In a first step, these boundary conditions can be chosen as constant. Later I would like to use time-dependent input.

How do I solve this problem with the Julia (DifferentialEquations)? In particular, what do I need to avoid **numerical diffusion** and get a numerically **stable system**?

> **Input**
>
> Here are some realistic input values:
> 
> ```julia
> L = 10e3
> x = [0,L]
> t = [0,43200]  
> p_0 = 50e5
> q_0 = 18 
> a = 1.3e6
> b = 1.9e5
> 
> ```

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [June 27, 2024, 6:47am UTC](https://discourse.julialang.org/t/solving-differential-equations-for-multi-component-fluid-flow/116254/2 "2024-06-27T06:47:53Z")

</div>

Not sure if this is your complete problem, or just a small starting point. I would suggest to consider using [Home · ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/) , because it allows you just to write down your differential equations in symbolic form and it simplifies your system of differential equations symbolically before solving it numerically.

If this is is your complete problem, using ModelingToolkit might be overkill, though.

---

<div class="post-metadata">

**Author:** ![Stephen](https://avatars.discourse-cdn.com/v4/letter/s/96bed5/32.png) [@Stephen](https://discourse.julialang.org/u/Stephen)\
**Post date:** [June 27, 2024, 10:43am UTC](https://discourse.julialang.org/t/solving-differential-equations-for-multi-component-fluid-flow/116254/3 "2024-06-27T10:43:16Z")

</div>

```julia

using MethodOfLines, ModelingToolkit, DomainSets

@parameters t, x, a, b
@variables p(..), q_m(..), c(..)  

Dt = Differential(t); Dx = Differential(x)

a, b = 1.3e6, 1.9e5 
p0, q0 = 5e6, 18

eq = [
    Dt(p(t, x)) ~ -a * Dx(q_m(t, x)), 
    Dx(p(t, x)) ~ -b * q_m(t, x) * abs(q_m(t, x)) / p(t, x),
    Dt(c(t, x)) ~ -a * q_m(t, x) / p(t, x) * Dx(c(t, x)), 
]

L = 1000.0

domain = [x ∈ Interval(0.0, L),
          t ∈ Interval(0.0, 43200.0)
]

ic_bc = [q_m(0.0, x) ~ q0,
         p(0, x) ~ sqrt(p0^2 - 2b * q0 * abs(q0) * x),
         c(0, x) ~ 1,
         q_m(t, L) ~ 1, # boundary condition be chosen as constant
         p(t, 0) ~ p0,
         c(t, 0) ~ 1 # boundary condition be chosen as constant
]

@named sys = PDESystem(eq, ic_bc, domain, [t, x], [p(t, x), q_m(t, x), c(t, x)])

dx = 2.0

discretization = MOLFiniteDifference([x => dx], t, approx_order = 2)

prob = discretize(sys, discretization)

using DifferentialEquations

sol = solve(prob, saveat = 1.0)

```

which throws error:

```julia
ERROR: MethodError: no method matching SciMLBase.PDETimeSeriesSolution(::SciMLBase.PDETimeSeriesSolution{Float64, 3, Dict{Num, Matrix{Float64}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, ModelingToolkit.MTKParameters{Tuple{}, Tuple{}, Tuple{}, Tuple{}, Tuple{}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out,), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x856c221d, 0x2ab738c3, 0x90c8aeff, 0x564502a8, 0x9bbfff35), Expr}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xc9c2ab17, 0x3395fe14, 0x86c47b71, 0xce7f7866, 0xb9d3bbeb), Expr}}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#740"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x244bcc7c, 0x4ad8c85d, 0xd984c452, 0xe91bc568, 0x17bbcbce), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x29d5008a, 0x0fae669a, 0xc0886d32, 0xb3a1ef51, 0x6483659f), Nothing}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}}, CompositeAlgorithm{1, Tuple{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{1, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, Rodas5P{1, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, FBDF{5, 1, false, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, typeof(OrdinaryDiffEq.trivial_limiter!)}, FBDF{5, 1, false, KrylovJL{typeof(Krylov.gmres!), Int64, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, typeof(OrdinaryDiffEq.trivial_limiter!)}}, OrdinaryDiffEq.AutoSwitchCache{Tuple{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, Tuple{Rosenbrock23{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, Rodas5P{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, FBDF{5, 0, false, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, typeof(OrdinaryDiffEq.trivial_limiter!)}, FBDF{5, 0, false, KrylovJL{typeof(Krylov.gmres!), Int64, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, typeof(OrdinaryDiffEq.trivial_limiter!)}}, Rational{Int64}, Int64}}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#740"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x244bcc7c, 0x4ad8c85d, 0xd984c452, 0xe91bc568, 0x17bbcbce), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x29d5008a, 0x0fae669a, 0xc0886d32, 0xb3a1ef51, 0x6483659f), Nothing}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, Vector{Int64}, OrdinaryDiffEq.DefaultCache{OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, uNoUnitsType, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False} where uNoUnitsType<:(Vector), OrdinaryDiffEq.Vern7Cache{Vector{Float64}, Vector{Float64}, uNoUnitsType, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False} where uNoUnitsType<:(Vector), OrdinaryDiffEq.Rosenbrock23Cache{Vector{Float64}, Vector{Float64}, uNoUnitsType, Matrix{Float64}, Matrix{Float64}, TabType, SciMLBase.TimeGradientWrapper{true, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#740"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x244bcc7c, 0x4ad8c85d, 0xd984c452, 0xe91bc568, 0x17bbcbce), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x29d5008a, 0x0fae669a, 0xc0886d32, 0xb3a1ef51, 0x6483659f), Nothing}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Vector{Float64}, ModelingToolkit.MTKParameters{Tuple{}, Tuple{}, Tuple{}, Tuple{}, Tuple{}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out,), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x856c221d, 0x2ab738c3, 0x90c8aeff, 0x564502a8, 0x9bbfff35), Expr}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xc9c2ab17, 0x3395fe14, 0x86c47b71, 0xce7f7866, 0xb9d3bbeb), Expr}}}, SciMLBase.UJacobianWrapper{true, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#740"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x244bcc7c, 0x4ad8c85d, 0xd984c452, 0xe91bc568, 0x17bbcbce), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x29d5008a, 0x0fae669a, 0xc0886d32, 0xb3a1ef51, 0x6483659f), Nothing}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Float64, ModelingToolkit.MTKParameters{Tuple{}, Tuple{}, Tuple{}, Tuple{}, Tuple{}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out,), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x856c221d, 0x2ab738c3, 0x90c8aeff, 0x564502a8, 0x9bbfff35), Expr}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xc9c2ab17, 0x3395fe14, 0x86c47b71, 0xce7f7866, 0xb9d3bbeb), Expr}}}, F, FiniteDiff.JacobianCache{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, UnitRange{Int64}, Nothing, Val{:forward}(), Float64}, FiniteDiff.GradientCache{Nothing, Vector{Float64}, Vector{Float64}, Float64, Val{:forward}(), Float64, Val{true}()}, Float64, Rosenbrock23{1, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, Vector{Bool}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)} where {uNoUnitsType<:(Vector), TabType<:OrdinaryDiffEq.Rosenbrock23Tableau, F<:(LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, LinearSolve.DefaultLinearSolver})}, OrdinaryDiffEq.Rosenbrock5Cache{Vector{Float64}, Vector{Float64}, uNoUnitsType, Matrix{Float64}, Matrix{Float64}, TabType, SciMLBase.TimeGradientWrapper{true, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#740"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x244bcc7c, 0x4ad8c85d, 0xd984c452, 0xe91bc568, 0x17bbcbce), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x29d5008a, 0x0fae669a, 0xc0886d32, 0xb3a1ef51, 0x6483659f), Nothing}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Vector{Float64}, ModelingTool…
Stacktrace:

wrap_sol(sol::SciMLBase.PDETimeSeriesSolution{…}, metadata::MethodOfLines.MOLMetadata{…}) at pde_solutions.jl

wrap_sol(sol::SciMLBase.PDETimeSeriesSolution{…}) at basic_solutions.jl

solve(::ODEProblem{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::Base.Pairs{…}) at solve.jl

kwcall(::@NamedTuple{…}, ::typeof(solve), ::ODEProblem{…}) at solve.jl

top-level scope at multi-component fluid flow.jl

```

---

<div class="post-metadata">

**Author:** ![DanaO](https://avatars.discourse-cdn.com/v4/letter/d/b5ac83/32.png) [@DanaO](https://discourse.julialang.org/u/DanaO)\
**Post date:** [June 27, 2024, 12:24pm UTC](https://discourse.julialang.org/t/solving-differential-equations-for-multi-component-fluid-flow/116254/4 "2024-06-27T12:24:57Z")

</div>

This is just a starting point for a larger problem. I’ve tried using the ModelingToolkit, but I haven’t managed to get a stable solution and the propagation stops early. The mass flow of the last time step shows some unphysical oscillations, it seems to me that the boundary conditions might be the problem, but I have not managed to fix it.

Code

```julia
using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets, Plots

function test2(L::Float64, T::Float64)
    a = 1.277e6
    b = 1.865e5
    q = 15.0
    p_in = 5e6
 
    @parameters x t
    @variables p(..) q_m(..) c(..)
    ∂t = Differential(t)
    ∂x = Differential(x)
 
    domains = [
        x ∈ Interval(0, L),
        t ∈ Interval(0, T),
    ]
 
    eq = [
        ∂t(p(x, t)) ~ (-a) * ∂x(q_m(x, t)),
        ∂x(p(x, t)) ~ - b * q_m(x, t) * abs(q_m(x, t)) / p(x, t),
        ∂t(c(x, t)) ~ - a * q_m(x, t) / p(x, t) * ∂x(c(x, t)),
    ]
 
    p_0(x) = sqrt((p_in)^2 - 2 * b * q * q * x)

    bcs = [
        # -------------------
        # Boundary conditions
        # -------------------
        p(0, t) ~ p_in, 
        q_m(L, t) ~ q, 
        c(0, t) ~ 1.0,
        # -------------------
        # Initial conditions
        # -------------------
        p(x, 0) ~ p_0(x),
        q_m(x, 0) ~ q_m0(x),
        c(x, 0) ~ 1.0,
    ]
 
    @named pde_system = PDESystem(
        eq, bcs, domains,
        [x, t],
        [p(x, t), q_m(x, t), c(x, t)],
    )
 
    return pde_system
end
 
pde_system = test2(500.0, 3600.0 * 1)
x, t = pde_system.ivs
p, q_m, c = pde_system.dvs
 
discretization = MOLFiniteDifference([x => 50], t; approx_order=2)
 
prob = discretize(pde_system, discretization)
sol = solve(prob, Rosenbrock23())

```

Output/Error

```julia
PDESystem
Equations: Equation[Differential(t)(p(x, t)) ~ -1.277e6Differential(x)(q_m(x, t)), Differential(x)(p(x, t)) ~ (-186500.0abs(q_m(x, t))*q_m(x, t)) / p(x, t), Differential(t)(c(x, t)) ~ (-1.277e6q_m(x, t)*Differential(x)(c(x, t))) / p(x, t)]
Boundary Conditions: Equation[p(0, t) ~ 5.0e6, q_m(500.0, t) ~ 15.0, c(0, t) ~ 1.0, p(x, 0) ~ sqrt(2.5e13 - 8.3925e7x), q_m(x, 0) ~ 15.0, c(x, 0) ~ 1.0]
Domain: Symbolics.VarDomainPairing[Symbolics.VarDomainPairing(x, 0.0 .. 500.0), Symbolics.VarDomainPairing(t, 0.0 .. 3600.0)]
Dependent Variables: Num[p(x, t), q_m(x, t), c(x, t)]
Independent Variables: Num[x, t]
Parameters: SciMLBase.NullParameters()
Default Parameter ValuesDict{Any, Any}()

2-element Vector{Num}:
 x
 t

3-element Vector{Num}:
   p(x, t)
 q_m(x, t)
   c(x, t)

MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}(Dict{Num, Int64}(x => 50), t, 2, UpwindScheme(1), MethodOfLines.CenterAlignedGrid(), true, false, MethodOfLines.ScalarizedDiscretization(), true, Any[], Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}())

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 3600.0)
u0: 147-element Vector{Float64}:
  4.999914361511503e6
  4.999828721556165e6
  4.999743080133912e6
  4.999657437244668e6
  4.999571792888356e6
  4.999486147064904e6
  4.999400499774233e6
  4.999314851016271e6
  4.99922920079094e6
  4.999143549098165e6
  4.999057895937871e6
  4.998972241309983e6
  4.998886585214425e6
  ⋮
 15.0
 15.0
 15.0
 15.0
 15.0
 15.0
 15.0
 15.0
 15.0
 15.0
 15.0
 15.0
 15.0

┌ Warning: At t=0.00031538038626707164, dt was forced below floating point epsilon 5.421010862427522e-20, and step error estimate = 3.677287959149432. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64).
└ @ SciMLBase C:\Users\OrsolitsD\.julia\packages\SciMLBase\sakPO\src\integrator_interface.jl:633
PDESolution:
  Return Code:
    Unstable
  Dependent variables:
    p(x, t): (50, 533) sized solution
    q_m(x, t): (50, 533) sized solution
    c(x, t): (50, 533) sized solution
  Domain:
    t ∈ (0.0, 0.00031538038626707164) with 533 non-uniform points. average step size 5.928202749381046e-7
    x ∈ (0.0, 500.0) with 50 points, step size 10.204081632653061
  From system:
    Equations:
L"\begin{align}
1.277 \cdot 10^{6} \frac{\mathrm{d}}{\mathrm{d}x} q_{m}\left( x, t \right) + \frac{\mathrm{d}}{\mathrm{d}t} p\left( x, t \right) =& 0 \\
\frac{1.865 \cdot 10^{5} q_{m}\left( x, t \right) \left|q_{m}\left( x, t \right)\right|}{p\left( x, t \right)} + \frac{\mathrm{d}}{\mathrm{d}x} p\left( x, t \right) =& 0 \\
\frac{1.277 \cdot 10^{6} q_{m}\left( x, t \right) \frac{\mathrm{d}}{\mathrm{d}x} c\left( x, t \right)}{p\left( x, t \right)} + \frac{\mathrm{d}}{\mathrm{d}t} c\left( x, t \right) =& 0
\end{align}
"
    Boundary/Initial Conditions:
L"\begin{align}
p\left( 0, t \right) =& 5 \cdot 10^{6} \\
q_{m}\left( 500, t \right) =& 15 \\
c\left( 0, t \right) =& 1 \\
p\left( x, 0 \right) =& \sqrt{2.5 \cdot 10^{13} - 8.3925 \cdot 10^{7} x} \\
q_{m}\left( x, 0 \right) =& 15 \\
c\left( x, 0 \right) =& 1
\end{align}

```

Plot Output

```julia
plot(sol[x][:], sol[q_m][:,end], label="mass flow")

```

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [June 27, 2024, 12:32pm UTC](https://discourse.julialang.org/t/solving-differential-equations-for-multi-component-fluid-flow/116254/5 "2024-06-27T12:32:06Z")

</div>

Some remarks:  
DifferentialEquations.jl supports supports hundreds of solvers, you might need one or two days to find the one that is best suited for your problem. In addition it might be important to set the absolute and relative accuracy of the solvers to a lower value, e.g. 1e-7.

I use (for a very different problem) code like this:

```julia
    prob = ODEProblem(m.sys_simple, u0, tspan, p)
    tol=1e-7
    sol = solve(prob, Rodas5(), dt=dt, tstops=tstops, abstol=tol, dtmax=1, reltol = tol, saveat=ts)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 29, 2024, 2:19am UTC](https://discourse.julialang.org/t/solving-differential-equations-for-multi-component-fluid-flow/116254/6 "2024-06-29T02:19:34Z")

</div>

Rodas5P or FBDF would be a better choice. Rosenbrock23 isn’t great at DAEs.

---

<div class="post-metadata">

**Author:** ![gstein3m](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gstein3m/32/45441_2.png) [@gstein3m](https://discourse.julialang.org/u/gstein3m)\
**Post date:** [June 30, 2024, 11:25am UTC](https://discourse.julialang.org/t/solving-differential-equations-for-multi-component-fluid-flow/116254/7 "2024-06-30T11:25:37Z")

</div>

It is possible that swapping the boundary conditions will lead to more stable results: The flow should be specified upstream and the pressure downstream.  
A weno discretization of the space derivatives could also be helpful.

And  
0 ~ ∂x(p(x, t)) + b \* q\_m(x, t) \* abs(q\_m(x, t)) / p(x, t),  
instead of  
∂x(p(x, t)) ~ - b \* q\_m(x, t) \* abs(q\_m(x, t)) / p(x, t),  
seems preferable.

---

<div class="post-metadata">

**Author:** ![DanaO](https://avatars.discourse-cdn.com/v4/letter/d/b5ac83/32.png) [@DanaO](https://discourse.julialang.org/u/DanaO)\
**Post date:** [July 4, 2024, 8:28am UTC](https://discourse.julialang.org/t/solving-differential-equations-for-multi-component-fluid-flow/116254/8 "2024-07-04T08:28:12Z")

</div>

I didn’t change the BC’s as they wouldn’t fit into the structure of my problem. However, rewriting the second equation and using FBDF()/Rodas5P() leads to a stable solution. I could not observe any difference for different advection schemes.