# Solving a system of PDEs with MOL memory allocation problems

Hello, I´m new user of Julia and I´m trying to solve a system of PDES non lineal with method of lines. To check it I´m using a static solution not dependent on time. The idea is to recover to t_i+1 the same solutions. I am going to enumerate the different difficulties and questions that arose and that I did not find answers in the tutorials and pages about the package.

The system is of second order passed to first order. The unknowns are `Φ, Λ, Ψ` and `γ`. With

`∂ₜ Φ = Φₜ , ∂ₜΨ = Ψₜ , ∂ₜΛ = Λₜ`

The code is

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

#static analytical solutions
Λ_est = x -> 2/(2-(4/5)*x^2*(5-3*x^2))

Φ_est = x -> -(4/5 - 1)*(-2*4/5 + 2)^(1/4)*exp.(-5/2*sqrt(4/5)*atan.((6*x^2 - 5)*sqrt.(4/5)*sqrt.(-25*4/5 + 24)/(25*4/5 - 24))/sqrt.(-25*4/5 + 24) + 5/2*sqrt(4/5)*atan.(sqrt.(4/5)*sqrt.(-25*4/5 + 24)/(25*4/5 - 24))/sqrt.(-25*4/5 + 24))/(3*4/5*x^4 - 5*4/5*x^2 + 2)^(1/4)

Ψ_est = x -> 1

γ_est = x -> 1.0/x * (1.0  - (2/(2-(4/5)*x^2*(5-3*x^2)))^2)

``````
``````@parameters t x
@variables Φ(..) Φₜ(..) Λ(..) Λₜ(..) Ψ(..) Ψₜ(..) γ(..)
Dt = Differential(t)
Dx = Differential(x)

x_min = 0.001  #is not 0 because the system is irregular at the origin
t_min = 0.0
x_max = 1.0
t_max = 0.1

# q is know
q = x -> (2/5)*(5.0-3.0*x^2)

eq = [Dx(Φ(t,x)) ~ (q(x)*x*Λ(t,x)*Φ(t,x))/(Ψ(t,x)^2),
Dx(Φₜ(t,x)) ~ (q(x)*x)/(Ψ(t,x)^3) * (Φ(t,x)*Ψ(t,x)*Λₜ(t,x) + Φₜ(t,x)*Ψ(t,x)*Λ(t,x) - 2 * Φ(t,x)*Ψₜ(t,x)*Λ(t,x)),
Dt(Ψ(t,x)) ~ Ψₜ(t,x),
Dt(Ψₜ(t,x)) ~ (Φ(t,x)^2/Λ(t,x)^2) * (Dx(Ψ(t,x))/x + (1/2 * γ(t,x) * Ψ(t,x))/x) - ((Ψₜ(t,x)^2)/(2*Ψ(t,x))) + Ψₜ(t,x)*Φₜ(t,x)/Φ(t,x) + Φ(t,x)^2 * ((Dx(Ψ(t,x)))^2/(2*Λ(t,x)*Ψ(t,x)) + q(x)*x*Dx(Ψ(t,x))/(Λ(t,x)*Ψ(t,x)^2) + q(x)/(Λ(t,x)*Ψ(t,x)) + q(x)^2*x^2 /(2*Ψ(t,x)^3)),
Dt(Λ(t,x)) ~ Λₜ(t,x),
Dt(Λₜ(t,x)) ~ -(Φ(t,x)^2/Λ(t,x))*(2*Dx(Ψ(t,x))/(x*Ψ(t,x)) + γ(t,x)/x) + Λ(t,x) * Ψₜ(t,x)/(Ψ(t,x)^2) - Λₜ(t,x) * Φₜ(t,x)/Φ(t,x) - (Φ(t,x)^2/Ψ(t,x)^2)* ( Dx(Ψ(t,x))^2/Λ(t,x) + 2*q(x)*x*Dx(Ψ(t,x))/Ψ(t,x) + 2*q(x) + q(x)^2 * x^2 * Λ(t,x)/(Ψ(t,x)^2)),
Dt(γ(t,x)) ~ (2*Λ(t,x)/Ψ(t,x)^3) *(Ψ(t,x)^3*Λₜ(t,x)+ q(x)*x*Λ(t,x)^2*Ψₜ(t,x) - Dx(Ψₜ(t,x)))
]
bcs=[
Ψ(t_min,x) ~ 1.0,
Ψₜ(t_min,x) ~ 0.0,
Λ(t_min,x) ~ 2/(2-(4/5)*x^2*(5-3*x^2)),
Λₜ(t_min,x) ~ 0.0,
γ(t_min,x) ~ (1.0/x) * (1.0  - (2/(2-(4/5)*x^2*(5-3*x^2)))^2),
Φ(t,x_min) ~ Φ_est(x_min),
Φₜ(t,x_min) ~ 0.0,

Φₜ(t_min,x) ~ 0.0,
Φ(t_min,x) ~ -(4/5 - 1)*(-2*4/5 + 2)^(1/4)*exp.(-5/2*sqrt.(4/5)*atan.((6*x^2 - 5)*sqrt.(4/5)*sqrt.(-25*4/5 + 24)/(25*4/5 - 24))/sqrt.(-25*4/5 + 24) + 5/2*sqrt(4/5)*atan.(sqrt.(4/5)*sqrt.(-25*4/5 + 24)/(25*4/5 - 24))/sqrt.(-25*4/5 + 24))/(3*4/5*x^4 - 5*4/5*x^2 + 2)^(1/4)
]

domains = [t ∈Interval(t_min, t_max),
x ∈Interval(x_min, x_max)]
@named pdesys = PDESystem(eq, bcs, domains,[t, x],[Φ(t,x), Φₜ(t,x), Ψ(t,x), Ψₜ(t,x), Λ(t,x), Λₜ(t,x),γ(t,x)])

N = 32
dx = (x_max -x_min)/N # intervalo de discretización
discx = chebyspace(N, domains)

discretization = MOLFiniteDifference([x => N],t)
@time prob = discretize(pdesys, discretization)
alg = Rosenbrock23()
@time sol = solve(prob,alg)
``````
``````discrete_x = sol[x]
discrete_t = sol[t]

solΦ = sol[Φ(t, x)]
solΦₜ = sol[Φₜ(t, x)]

solΛ = sol[Λ(t, x)]
solΛₜ = sol[Λₜ(t, x)]

solΨ = sol[Ψ(t, x)]
solΨₜ = sol[Ψₜ(t, x)]
``````
``````Φ_real = [Φ_est(s) for s in discrete_x]
errΦ = scatter()
for i in 1:length(discrete_t)
scatter!(discrete_x, solΦ[i, :].-Φ_real,title= "Error en Φ", label="t=\$(discrete_t[i])")
end
display(errΦ)
``````

The code runs, one of the problems I have is that if I want to overcome the barrier of N greater than 100 I have memory allocation problems and it crashes. Do you know why this happens?. Is it because of the complexity of the system? Or is there something that as a newbie I don’t do (I must clarify that I do not have a great computer).
On the other hand, if I try to use a smaller grid but with a chebyshev-type discretization Non-Uniform Rectilinear Grids · MethodOfLines.jl , when I run the pdes discretization, I get the following error:
`BoundsError: attempt to access 32-element Vector{StaticArraysCore.SVector{2, Float64}} at index `

Lastly, as can be seen, I do not have an evolution equation for Φ, so errors accumulate in its calculation and it does not give a static solution. My question is why when I give only the bc Φ(x=0.001,t) = Φ_est(0.001) and Φ_t(x=0.001,t) = 0 it does not solve for Phi and Phi_t ? . But if I give it those conditions plus Φ(x,t=0) = Φ_est and Φ(x,t=0) = 0.0 if it gives me valid results

I would appreciate if at least one doubt can be resolved in advance
Sorry if my writing is bad, english is not my native language, thank you very much

You’ll need to try it with JuliaSimCompiler which should fix the scaling issue, once outlining lands at least.

https://help.juliahub.com/juliasimcompiler/dev/

This is very new stuff though so it’s a bit greenfield and no guarantees yet. We’re still benchmarking it on PDE examples like this.