# Solution not converging for advection PDE

**URL:** <https://discourse.julialang.org/t/solution-not-converging-for-advection-pde/123349>\
**Category:** Numerics\
**Tags:** pde, differentialequation\
**Created:** [December 2, 2024, 9:58am UTC](https://discourse.julialang.org/t/solution-not-converging-for-advection-pde/123349 "2024-12-02T09:58:08Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![Alma\_Thindberg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alma_thindberg/32/213715_2.png) [@Alma\_Thindberg](https://discourse.julialang.org/u/Alma_Thindberg)\
**Post date:** [December 2, 2024, 9:58am UTC](https://discourse.julialang.org/t/solution-not-converging-for-advection-pde/123349/1 "2024-12-02T09:58:08Z")

</div>

```julia
using Pkg
#Pkg.add("ModelingToolkit")
#Pkg.add("MethodOfLines")
#Pkg.add("DomainSets")
#Pkg.add("OrdinaryDiffEq") 

Pkg.add("NonlinearSolve")

using ModelingToolkit, MethodOfLines, DomainSets, OrdinaryDiffEq, LinearAlgebra, NonlinearSolve
include("Properties.jl")

import ModelingToolkit: scalarize

#substances = [CH4, H2, CO2, CO]

Nc = 4
R = 8.314
Nf = 60000
D0 = 250*10^-6
L = 0.6
k = 1*10^-4
Tref = 273.15
Pref = 1*10^5 
T = 313.15
Di = 200*10^-6
D = 0.1
Pf_0 = 20*10^5 
Pp_L = 1*10^5 
a_i = [5.33*10^-12, 1.6*10^-10, 6.956*10^-11, 6.66*10^-12] 
Mw_i = [16.04, 2.016, 44.01, 28.01] 
μ_i = [7.29956708200112e-5, 9.211635626733832e-6, 1.568221747830425e-5, 1.8359577261147434e-5]
x_i_0 = [0.04, 0.75, 0.20, 0.01]
Ff_0 = 0.03 

abstract type WilkeChangMixingRule end
function μ_mix(mixing_rule::Type{T}, y, Mw, μᵢ) where T <: WilkeChangMixingRule
    _1_Mw = scalarize(1.0./Mw) 
    num = scalarize(y.*μᵢ)
    sqrt_Mw_mat = scalarize(.√(collect(Mw)*collect(_1_Mw')))
    _den = scalarize(y.*sqrt_Mw_mat)
    den = scalarize(sum([_den[i, :] for i in 1:size(_den, 1)]))
    return scalarize(sum(num./den))
end

@parameters z

# Define variables with updated array-valued syntax
vars1 = @variables Fₚ(..) Ff(..) Pf(..) Pp(..) 
vars2 = @variables (Fp_i(..))[1:Nc] (Ff_i(..))[1:Nc] (J_i(..))[1:Nc] (x_i(..))[1:Nc] (y_i(..))[1:Nc]
vars3 = @variables (xm_i(z))[1:Nc] (A_i(z))[1:Nc] (N_i(z))[1:Nc] (ρᵢ(z))[1:Nc] Nt(z) X(z)

callable_scalar = [Fₚ, Ff, Pf, Pp]

var_scalars = [var(z) for var in callable_scalar]

var_vectors = [collect(var(z)) for var in vars2]

vars = [var_scalars...; var_vectors...; vars3...]

Dz = Differential(z)

# Define the system of equations
eqn = [
    scalarize(J_i(z) .~ a_i .* (Pf(z) .* xm_i .- Pp(z) .* y_i(z)))...
    scalarize(N_i .~ A_i .* (Pf(z) .* xm_i .- Pp(z) .* y_i(z)))...
    scalarize(N_i .~ k .* (x_i(z) .- xm_i) .* X .+ x_i(z) .* Nt)...
    scalarize(A_i .~ a_i .* Mw_i ./ ρᵢ)...
    scalarize(ρᵢ .~ (x_i(z) .* Pp(z)) ./ (R .* T .+ 1*10^-10))... 
    Nt ~ scalarize(sum(N_i))
    X ~ ((Pf(z))+1*10^-10 / Pref) * (Tref / T)
    scalarize(x_i(z) .~ Ff_i(z) ./ Ff(z))...
    scalarize(y_i(z) .~ Fp_i(z) ./ Fₚ(z))...
    scalarize(Dz(Ff_i(z)) .~ -J_i(z) .* π .* D0 .* Nf)...  
    scalarize(Dz(Fp_i(z)) .~ -J_i(z) .* π .* D0 .* Nf)...   
    Dz(Pf(z)) ~ (-8 * μ_mix(WilkeChangMixingRule, x_i(z), Mw_i, μ_i)* L)*0.03/((Di/2)^2)
    Dz(Pp(z)) ~ (8 * μ_mix(WilkeChangMixingRule, y_i(z), Mw_i, μ_i)* L)*0.03/((Di/2)^2)
    Fₚ(z) ~ scalarize(sum(Fp_i(z)))
    Ff(z) ~ scalarize(sum(Ff_i(z)))
]

# Integral limits are defined with DomainSets.ClosedInterval
Ix = Integral(z in DomainSets.ClosedInterval(0.0, L)) 

bcs=[
    Pf(0) ~ Pf_0
    Pp(0) ~ Pp_L
    [scalarize(Ff_i(0) .~ Ff_0 * x_i_0)]...
    [Fp_i(0)[i] .~ Ix(J_i(z)[i])* π * D0 * Nf for i in 1:Nc]...

]

domains = [z ∈ Interval(0.0, L)]

@named pdesys = PDESystem(eqn, bcs, domains, z, vars, [])

Nz = 20
dz = L/Nz

discretization = MOLFiniteDifference([z => dz], nothing, approx_order = 2)

prob = discretize(pdesys, discretization, should_transform = true)

sol = solve(prob, NewtonRaphson())

sol = solve(prob, NonlinearSolve.FastShortcutNonlinearPolyalg(), maxiters=50000, abstol=1e-6)

using Plots

# Generate z values (spatial grid points)
z_values = range(0, stop=0.6, length=length(sol.u[Pf(z)]))

# Substances for labeling
substances = ["CH4", "H2", "CO2", "CO"]

# 1. Plot fluxes J_i for all components
plt_fluxes = plot()
for i in 1:Nc
    plot!(z_values, sol.u[J_i(z)[i]], label=substances[i])
end
xlabel!("z [m]")
ylabel!("J_i [mol/s]")
title!("Fluxes J_i vs z")
legend=:topright

# 2. Plot x_i and y_i mole fractions for all components
plt_mole_fractions = plot()
for i in 1:Nc
    plot!(z_values, sol.u[x_i(z)[i]], label="x_" * substances[i] * " (Feed)")
    plot!(z_values, sol.u[y_i(z)[i]], linestyle=:dash, label="y_" * substances[i] * " (Permeate)")
end
xlabel!("z [m]")
ylabel!("x_i, y_i [mol/mol]")
title!("Mole Fractions x_i and y_i vs z")
legend=:topright

# 3. Plot pressure differences ΔP_f and ΔP_p
deltaP_f = Pf_0 .- sol.u[Pf(z)]
deltaP_p = (sol.u[Pp(z)] .- Pp_L)

plt_pressure_differences = plot()
plot!(z_values, deltaP_f, label="ΔP_f = P_{f,0} - P_f")
plot!(z_values, deltaP_p, label="ΔP_p = P_{p,0} - P_p")
xlabel!("z [m]")
ylabel!("Pressure Difference [Pa]")
title!("Pressure Differences vs z")
legend=:topright

# Combine all plots into one figure
combined_plot = plot(plt_fluxes, plt_mole_fractions, plt_pressure_differences, layout=(3,1), size=(800, 1200))

display(combined_plot)

println(sol.u[Pf(z)])
println(sol.u[Pp(z)])
# Generate z values (spatial grid points)
z_values = range(0, stop=0.6, length=length(sol.u[Pp(z)]))

# Extract Pp(z) values
Pp_values = sol.u[Pp(z)]

# Print z and Pp(z) together
println("z [m] \t Pp [Pa]")
for (z, Pp) in zip(z_values, Pp_values)
    println("$z \t $Pp")
end

```

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [December 2, 2024, 9:59am UTC](https://discourse.julialang.org/t/solution-not-converging-for-advection-pde/123349/2 "2024-12-02T09:59:29Z")

</div>

You will likely get better help if you describe the problem you have, what your expected results are and what help you would like, rather than just posting a code snippet (that does not even run due to missing file `Properties.jl`).

---

<div class="post-metadata">

**Author:** ![Alma\_Thindberg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alma_thindberg/32/213715_2.png) [@Alma\_Thindberg](https://discourse.julialang.org/u/Alma_Thindberg)\
**Post date:** [December 2, 2024, 10:58am UTC](https://discourse.julialang.org/t/solution-not-converging-for-advection-pde/123349/3 "2024-12-02T10:58:23Z")

</div>

It runs without the code Properties.jl. I added the code snippet from that file in the code shown. But you are correct in that the “include line” should have been commented out before posting, thank you!

---

<div class="post-metadata">

**Author:** ![Alma\_Thindberg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alma_thindberg/32/213715_2.png) [@Alma\_Thindberg](https://discourse.julialang.org/u/Alma_Thindberg)\
**Post date:** [December 2, 2024, 11:10am UTC](https://discourse.julialang.org/t/solution-not-converging-for-advection-pde/123349/4 "2024-12-02T11:10:50Z")

</div>

I’m working on simulating membrane separation using the method of lines (MOL) for solving a system of PDEs. I’ve encountered an issue with the solver:

- When I attempt to solve the system using `NewtonRaphson()` as the nonlinear solver, it throws a `MaxIters` error, indicating it doesn’t converge.
- On the other hand, the system solves successfully when I switch to `NonlinearSolve.FastShortcutNonlinearPolyalg()` with increased `maxiters` and a higher absolute tolerance (`abstol=1e-6`).

My goal is to make the code work with `NewtonRaphson()` at a lower tolerance. I’d like to understand:

1. Why the system fails to converge with `NewtonRaphson()`.
2. What adjustments I can make (e.g., in initial guesses, tolerances, scaling, or other settings) to ensure convergence with `NewtonRaphson()`.

![image](https://global.discourse-cdn.com/julialang/optimized/3X/0/2/02dc09563413128afd49047705456582cc33e1e1_2_690x66.png)

---

<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:** [December 3, 2024, 6:32pm UTC](https://discourse.julialang.org/t/solution-not-converging-for-advection-pde/123349/5 "2024-12-03T18:32:51Z")

</div>

> [@Alma\_Thindberg](#):
>
> When I attempt to solve the system using `NewtonRaphson()` as the nonlinear solver, it throws a `MaxIters` error, indicating it doesn’t converge.

Likely the initial condition for the nonlinear solver is starting at a singularity, and thus the Jacobian is having issues. NewtonRaphson is generally not a good idea, I’m curious why choose that instead of say TrustRegion?

---

<div class="post-metadata">

**Author:** ![Alma\_Thindberg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alma_thindberg/32/213715_2.png) [@Alma\_Thindberg](https://discourse.julialang.org/u/Alma_Thindberg)\
**Post date:** [December 4, 2024, 9:04am UTC](https://discourse.julialang.org/t/solution-not-converging-for-advection-pde/123349/6 "2024-12-04T09:04:21Z")

</div>

Hi, thank you for the input! I tried using TrustRegion, but I am still having issues with convergence.  
I get the following error:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/e/8/e8325dbc15ccd6ea9656f7637816a804d027f5ea.png)

Do you have any tips on how to move forward and how to give initial guesses with TrustRegion or MethodOfLines?

```julia
PDESolution:
  Return Code:
    MaxIters
  Dependent variables:
    (x_i(z))[2]: (21,) sized solution
    (Fp_i(z))[2]: (21,) sized solution
    var"A_i_Any[3]"(z): (21,) sized solution
    var"J_i_Any[4]"(z): (21,) sized solution
    (Ff_i(z))[3]: (21,) sized solution
    var"Ff_i_Any[3]"(z): (21,) sized solution
    var"xm_i_Any[3]"(z): (21,) sized solution
    (y_i(z))[2]: (21,) sized solution
    (Ff_i(z))[4]: (21,) sized solution
    var"N_i_Any[3]"(z): (21,) sized solution
    Fₚ(z): (21,) sized solution
    (J_i(z))[3]: (21,) sized solution
    (J_i(z))[4]: (21,) sized solution
    var"xm_i_Any[2]"(z): (21,) sized solution
    var"Fp_i_Any[2]"(z): (21,) sized solution
    var"Fp_i_Any[4]"(z): (21,) sized solution
    var"Ff_i_Any[4]"(z): (21,) sized solution
    Pp(z): (21,) sized solution
    var"y_i_Any[4]"(z): (21,) sized solution
    var"ρᵢ_Any[4]"(z): (21,) sized solution
    var"Fp_i_Any[3]"(z): (21,) sized solution
    var"J_i_Any[3]"(z): (21,) sized solution
    (xm_i(z))[1]: (21,) sized solution
    (A_i(z))[3]: (21,) sized solution
    var"Ff_i_Any[1]"(z): (21,) sized solution
    (ρᵢ(z))[1]: (21,) sized solution
    (A_i(z))[4]: (21,) sized solution
    (N_i(z))[3]: (21,) sized solution
    var"ρᵢ_Any[3]"(z): (21,) sized solution
    (N_i(z))[4]: (21,) sized solution
    (xm_i(z))[2]: (21,) sized solution
    var"J_i_Any[2]"(z): (21,) sized solution
    (ρᵢ(z))[2]: (21,) sized solution
    var"xm_i_Any[1]"(z): (21,) sized solution
    var"x_i_Any[4]"(z): (21,) sized solution
    (Ff_i(z))[1]: (21,) sized solution
    Ff(z): (21,) sized solution
    var"x_i_Any[2]"(z): (21,) sized solution
    (J_i(z))[1]: (21,) sized solution
    var"A_i_Any[4]"(z): (21,) sized solution
    var"Fp_i_Any[1]"(z): (21,) sized solution
    var"N_i_Any[1]"(z): (21,) sized solution
    (Ff_i(z))[2]: (21,) sized solution
    var"A_i_Any[2]"(z): (21,) sized solution
    (x_i(z))[3]: (21,) sized solution
    (Fp_i(z))[3]: (21,) sized solution
    X(z): (21,) sized solution
    (J_i(z))[2]: (21,) sized solution
    (Fp_i(z))[4]: (21,) sized solution
    (x_i(z))[4]: (21,) sized solution
    var"Ff_i_Any[2]"(z): (21,) sized solution
    (y_i(z))[3]: (21,) sized solution
    (A_i(z))[1]: (21,) sized solution
    var"ρᵢ_Any[1]"(z): (21,) sized solution
    (y_i(z))[4]: (21,) sized solution
    var"y_i_Any[2]"(z): (21,) sized solution
    (N_i(z))[1]: (21,) sized solution
    var"A_i_Any[1]"(z): (21,) sized solution
    Pf(z): (21,) sized solution
    var"xm_i_Any[4]"(z): (21,) sized solution
    var"y_i_Any[3]"(z): (21,) sized solution
    (A_i(z))[2]: (21,) sized solution
    Nt(z): (21,) sized solution
    var"N_i_Any[2]"(z): (21,) sized solution
    (N_i(z))[2]: (21,) sized solution
    var"J_i_Any[1]"(z): (21,) sized solution
    var"x_i_Any[1]"(z): (21,) sized solution
    (xm_i(z))[3]: (21,) sized solution
    (Fp_i(z))[1]: (21,) sized solution
    (x_i(z))[1]: (21,) sized solution
    (xm_i(z))[4]: (21,) sized solution
    (ρᵢ(z))[3]: (21,) sized solution
    var"ρᵢ_Any[2]"(z): (21,) sized solution
    var"N_i_Any[4]"(z): (21,) sized solution
    var"x_i_Any[3]"(z): (21,) sized solution
    (ρᵢ(z))[4]: (21,) sized solution
    var"y_i_Any[1]"(z): (21,) sized solution
    (y_i(z))[1]: (21,) sized solution
  Domain:
    z ∈ (0.0, 0.6) with 21 points, step size 0.03
  From system:
    Equations:
"\\begin{align}\n\\mathtt{J\\_i\\_Any[1]}\\left( z \\right) - 5.33 \\cdot 10^{-12} \\left( - \\mathtt{Pp}\\left( z \\right) \\mathtt{y\\_i\\_Any[1]}\\left( z \\right) + \\mathtt{Pf}\\left( z \\right) \\mathtt{xm\\_i\\_Any[1]}\\left( z \\right) \\right) &= 0 \\\\\n\\mathtt{J\\_i\\_Any[2]}\\left( z \\right) - 1.6 \\cdot 10^{-10} \\left( \\mathtt{xm\\_i\\_Any[2]}\\left( z \\right) \\mathtt{Pf}\\left( z \\right) - \\mathtt{Pp}\\left( z \\right) \\mathtt{y\\_i\\_Any[2]}\\left( z \\right) \\right) &= 0 \\\\\n\\mathtt{J\\_i\\_Any[3]}\\left( z \\right) - 6.956 \\cdot 10^{-11} \\left( - \\mathtt{Pp}\\left( z \\right) \\mathtt{y\\_i\\_Any[3]}\\left( z \\right) + \\mathtt{Pf}\\left( z \\right) \\mathtt{xm\\_i\\_Any[3]}\\left( z \\right) \\right) &= 0 \\\\\n\\mathtt{J\\_i\\_Any[4]}\\left( z \\right) - 6.66 \\cdot 10^{-12} \\left( \\mathtt{xm\\_i\\_Any[4]}\\left( z \\right) \\math" ⋯ 7114 bytes ⋯ " z \\right) + 4.6723 \\mathtt{y\\_i\\_Any[3]}\\left( z \\right) + 2.8207 \\mathtt{y\\_i\\_Any[1]}\\left( z \\right)} + \\frac{7.2996 \\cdot 10^{-5} \\mathtt{y\\_i\\_Any[1]}\\left( z \\right)}{1.3215 \\mathtt{y\\_i\\_Any[4]}\\left( z \\right) + 0.35452 \\mathtt{y\\_i\\_Any[2]}\\left( z \\right) + 1.6564 \\mathtt{y\\_i\\_Any[3]}\\left( z \\right) + \\mathtt{y\\_i\\_Any[1]}\\left( z \\right)} \\right) &= 0 \\\\\n\\mathtt{F_p}\\left( z \\right) - \\mathtt{Fp\\_i\\_Any[2]}\\left( z \\right) - \\mathtt{Fp\\_i\\_Any[3]}\\left( z \\right) - \\mathtt{Fp\\_i\\_Any[4]}\\left( z \\right) - \\mathtt{Fp\\_i\\_Any[1]}\\left( z \\right) &= 0 \\\\\n - \\mathtt{Ff\\_i\\_Any[3]}\\left( z \\right) - \\mathtt{Ff\\_i\\_Any[4]}\\left( z \\right) - \\mathtt{Ff\\_i\\_Any[2]}\\left( z \\right) + \\mathtt{Ff}\\left( z \\right) - \\mathtt{Ff\\_i\\_Any[1]}\\left( z \\right) &= 0\n\\end{align}\n"
    Boundary/Initial Conditions:
L"\begin{align}
\mathtt{Pf}\left( 0 \right) &= 2000000 \\
\mathtt{Pp}\left( 0 \right) &= 100000 \\
\mathtt{Ff\_i\_Any[1]}\left( 0 \right) &= 0.0012 \\
\mathtt{Ff\_i\_Any[2]}\left( 0 \right) &= 0.0225 \\
\mathtt{Ff\_i\_Any[3]}\left( 0 \right) &= 0.006 \\
\mathtt{Ff\_i\_Any[4]}\left( 0 \right) &= 0.0003 \\
\mathtt{Fp\_i\_Any[1]}\left( 0 \right) &= 47.124 \int_{0.0}^{0.6)} z \mathtt{J\_i\_Any[1]}\left( z \right) \\
\mathtt{Fp\_i\_Any[2]}\left( 0 \right) &= 47.124 \int_{0.0}^{0.6)} z \mathtt{J\_i\_Any[2]}\left( z \right) \\
\mathtt{Fp\_i\_Any[3]}\left( 0 \right) &= 47.124 \int_{0.0}^{0.6)} z \mathtt{J\_i\_Any[3]}\left( z \right) \\
\mathtt{Fp\_i\_Any[4]}\left( 0 \right) &= 47.124 \int_{0.0}^{0.6)} z \mathtt{J\_i\_Any[4]}\left( z \right)
\end{align}
"

```

---

<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:** [December 20, 2024, 12:02am UTC](https://discourse.julialang.org/t/solution-not-converging-for-advection-pde/123349/8 "2024-12-20T00:02:48Z")

</div>

Can you open an issue on MethodOfLines.jl? I don’t want to lose this.

---

<div class="post-metadata">

**Author:** ![Alma\_Thindberg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alma_thindberg/32/213715_2.png) [@Alma\_Thindberg](https://discourse.julialang.org/u/Alma_Thindberg)\
**Post date:** [March 17, 2025, 9:17am UTC](https://discourse.julialang.org/t/solution-not-converging-for-advection-pde/123349/9 "2025-03-17T09:17:28Z")

</div>

Hi, I’m sorry for the delayed response! I wanted to let you know that I managed to solve the issue. It turned out the problem wasn’t with the code after all… I had an incorrect parameter approximation that caused unexpected behavior.

---

<div class="post-metadata">

**Author:** ![kkakosim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kkakosim/32/21546_2.png) [@kkakosim](https://discourse.julialang.org/u/kkakosim)\
**Post date:** [March 18, 2025, 5:47am UTC](https://discourse.julialang.org/t/solution-not-converging-for-advection-pde/123349/10 "2025-03-18T05:47:34Z")

</div>

Hey @Alma_Thindberg, do you mind posting the final code that is running? we are facing similar challenges from time and it is good to have working examples for verification.

---

<div class="post-metadata">

**Author:** ![Alma\_Thindberg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alma_thindberg/32/213715_2.png) [@Alma\_Thindberg](https://discourse.julialang.org/u/Alma_Thindberg)\
**Post date:** [March 25, 2025, 1:05pm UTC](https://discourse.julialang.org/t/solution-not-converging-for-advection-pde/123349/11 "2025-03-25T13:05:01Z")

</div>

Hi, yes of course!  
The code is not “cleaned up”, but this version converges and yields results that look close to expected values.

```julia
using Pkg

Pkg.add("NonlinearSolve")
Pkg.add("NonlinearSolveFirstOrder")
Pkg.add("LaTeXStrings")

using ModelingToolkit, MethodOfLines, DomainSets, OrdinaryDiffEq, LinearAlgebra, NonlinearSolve, NonlinearSolveFirstOrder, Statistics  
include("Properties.jl")

import ModelingToolkit: scalarize

substances = [CH4, H2, CO2, CO]
Nc = 4
R = 8.3145
Nf = 60000
D0 = 250*10^-6
L = 0.6
k = 1*10^-4
Tref = 273.15
Pref = 1*10^5 
T = 313.15
Di = 200*10^-6
D = 0.1
Pf_0 = 20*10^5 
Pp_0 = 1*10^5 
a_i = [5.33*10^-12, 1.6*10^-10, 6.956*10^-11, 6.66*10^-12] 
Mw_i = [16.04, 2.016, 44.01, 28.01] 
μ_i = [] 
x_i_0 = [0.04, 0.75, 0.20, 0.01]
Ff_0 = 0.015 
μ_methane = 12.126 *10^-6

for substance in substances
    println("Processing substance: ", substance.infos.name)

    if substance.infos.name == "methane"
        μ_value = μ_methane

    else
        viscosity_data = get_viscosity_data(substance)  
        μ_value = μ(viscosity_data, T) 
    end

    println("Calculated viscosity: ", μ_value)
    push!(μ_i, μ_value)  
end

println(μ_i)

@parameters z
vars1 = @variables Fₚ(..) Ff(..) Pf(..) Pp(..) 
vars2 = @variables (Fp_i(..))[1:Nc] (Ff_i(..))[1:Nc] (J_i(..))[1:Nc] (x_i(..))[1:Nc] (y_i(..))[1:Nc]
vars3 = @variables (xm_i(z))[1:Nc] (A_i(z))[1:Nc] (N_i(z))[1:Nc] (ρᵢ(z))[1:Nc] Nt(z) X(z)

callable_scalar = [Fₚ, Ff, Pf, Pp]
var_scalars = [var(z) for var in callable_scalar]
var_vectors = [collect(var(z)) for var in vars2]
vars = [var_scalars...; var_vectors...; vars3...]

Dz = Differential(z)
eqn = [
    scalarize(J_i(z) .~ a_i .* (Pf(z) .* xm_i .- Pp(z) .* y_i(z)))...
    scalarize(N_i .~ A_i .* (Pf(z) .* xm_i .- Pp(z) .* y_i(z)))...
    scalarize(N_i .~ k .* (x_i(z) .- xm_i) .* X .+ x_i(z) .* Nt)...
    scalarize(A_i .~ a_i ./ ρᵢ)...
    scalarize(ρᵢ .~ (x_i(z) .* Pf(z)) ./ (R .* T))... #add molar weight
    Nt ~ scalarize(sum(N_i))
    X ~ ((Pf(z))/ Pref) * (Tref / T)
    scalarize(x_i(z) .~ Ff_i(z) ./ Ff(z))...
    scalarize(y_i(z) .~ Fp_i(z) ./ Fₚ(z))...
    scalarize(Dz(Ff_i(z)) .~ -J_i(z) .* π .* D0 .* Nf)...  
    scalarize(Dz(Fp_i(z)) .~ -J_i(z) .* π .* D0 .* Nf)...
    Dz(Pf(z)) ~ (-192.0 * Nf * D0 * (D + Nf * D0) * R * T * μ_mix(WilkeChangMixingRule, x_i(z), Mw_i, μ_i) * Ff(z) ) / (π * (D^2 - Nf * D0^2)^3 * Pf(z)) # Call μ_mix with arguments
    Dz(Pp(z)) ~ (128.0 * R * T * μ_mix(WilkeChangMixingRule, y_i(z), Mw_i, μ_i) * Fₚ(z)) / (π * Di^4 * Nf*Pp(z))     
    Ff(z) ~ scalarize(sum(Ff_i(z)))
    Fₚ(z) ~ scalarize(sum(Fp_i(z)))
]

Ix = Integral(z in DomainSets.ClosedInterval(0.0, L)) 

bcs=[
    Pf(0) ~ Pf_0
    Pp(0) ~ Pp_0
    [scalarize(Ff_i(0) .~ Ff_0*x_i_0)]...
    [Fp_i(0)[i] .~ Ix(J_i(z)[i])* π * D0 * Nf for i in 1:Nc]...
]

domains = [z ∈ Interval(0.0, L)]
@named pdesys = PDESystem(eqn, bcs, domains, z, vars, [])

Nz = 30
dz = L/Nz

discretization = MOLFiniteDifference([z => dz], nothing, approx_order = 2)

prob = discretize(pdesys, discretization, should_transform = true)

sol = solve(prob, NonlinearSolve.FastShortcutNonlinearPolyalg(), maxiters=80000,abstol=1e-6)

```

Properties.jl

```julia
import ModelingToolkit: scalarize

abstract type WilkeChangMixingRule end

mutable struct ViscosityPoly{T}
    A::T
    B::T
    C::T
    D::T
end

function ViscosityPoly(;A, B, C, D)
    return ViscosityPoly(A, B, C, D)
end

mutable struct FluidInfos
    name::String
    Mw::Number
end

function FluidInfos(;name, Mw)
    return FluidInfos(name, Mw)
end

mutable struct SimpleMediumData{T}
    viscosity_data::ViscosityPoly{T} 
end

function SimpleMediumData(;viscosity_data)
    return SimpleMediumData(viscosity_data)
end

mutable struct SimpleGasMedium
    infos::FluidInfos
    data::SimpleMediumData
end

function SimpleGasMedium(;infos, data)
    return SimpleGasMedium(infos, data)
end

function μ(viscosity_model::ViscosityPoly, T)
   # _1 = eltype(T)
   _1 = one(T)
    A, B, C, D = viscosity_model.A, viscosity_model.B, viscosity_model.C, viscosity_model.D
    return A * T^B / (_1 + C/T + D/T^2)
end

function get_viscosity_data(medium::SimpleGasMedium)
    return medium.data.viscosity_data
end

function μ_mix(mixing_rule::Type{T}, y, Mw, μᵢ) where T <: WilkeChangMixingRule
    _1_Mw = scalarize(1.0./Mw) 
    num = scalarize(y.*μᵢ)
    sqrt_Mw_mat = scalarize(.√(collect(Mw)*collect(_1_Mw')))
    _den = scalarize(y.*sqrt_Mw_mat)
    den = scalarize(sum([_den[i, :] for i in 1:size(_den, 1)]))
    return scalarize(sum(num./den))
end

const CH4 = SimpleGasMedium(infos = FluidInfos(name = "methane", Mw = 16.04),
 data = SimpleMediumData(viscosity_data = ViscosityPoly(A=2.148*10^-6, B=0.685, C=159.0, D=0.0)))

const H2 = SimpleGasMedium(
    infos = FluidInfos(name = "hydrogen", Mw = 2.016),
    data = SimpleMediumData(viscosity_data = ViscosityPoly(A=1.797*10^-7, B=0.685, C=-0.59, D=140.0))
)

const CO2 = SimpleGasMedium(
    infos = FluidInfos(name = "carbondioxide", Mw = 44.01),
    data = SimpleMediumData(viscosity_data = ViscosityPoly(A=2.148*10^-6, B=0.46, C=290.0, D=0.0))
)

const CO = SimpleGasMedium(
    infos = FluidInfos(name = "carbonmonoxide", Mw = 28.01),
    data = SimpleMediumData(viscosity_data = ViscosityPoly(A=1.1127*10^-6, B=0.5338, C=94.7, D=0.0))
)

```