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
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).
3 Likes
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!
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 aMaxIterserror, indicating it doesn’t converge. - On the other hand, the system solves successfully when I switch to
NonlinearSolve.FastShortcutNonlinearPolyalg()with increasedmaxitersand 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:
- Why the system fails to converge with
NewtonRaphson(). - What adjustments I can make (e.g., in initial guesses, tolerances, scaling, or other settings) to ensure convergence with
NewtonRaphson().
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?
Hi, thank you for the input! I tried using TrustRegion, but I am still having issues with convergence.
I get the following error:
Do you have any tips on how to move forward and how to give initial guesses with TrustRegion or MethodOfLines?
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}
"
Can you open an issue on MethodOfLines.jl? I don’t want to lose this.