# "equation(s) resulted in a non-finite number.." using MethodOfLines

Hi all,
It’s my first time using MethodOfLines so I’m not sure how to debug it. The (nondimensionalized) equations I’m trying to solve are:

\begin{aligned} & \frac{\partial u}{\partial x}+\frac{\partial w}{\partial z}=0 \\ & \frac{\partial P}{\partial x}=\frac{\partial \tau_1}{\partial x}+\frac{\partial \tau_3}{\partial z} \\ & \frac{\partial P}{\partial z}=\frac{\partial \tau_3}{\partial x}-\frac{\partial \tau_1}{\partial z}-\operatorname{Ra}(1-T), \\ & \tau_1=2 \eta \frac{\partial u}{\partial x}, \\ & \tau_3=\eta\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right), \\ \frac{\partial T}{\partial t}+\boldsymbol{u} \cdot \nabla T+D T w &= \nabla^2 T+\frac{D}{\operatorname{Ra}} \frac{\tau^2}{2 \eta}, \\ \eta=\exp \left[\frac{1-T+\mu(1-z-T)}{\varepsilon T}\right] \end{aligned}

in a 2D Cartesian box, with the following boundary conditions:

My approach is:

@parameters t,x,z
@variables u(..), w(..), P(..), T(..), σ₁(..), σ₃(..), η(..);

Dt = Differential(t)
Dx = Differential(x)
Dz = Differential(z)
Dxx = Differential(x)^2
Dzz = Differential(z)^2

D = 0.6
Ra = 1e7
θ = 0.1
ε = 0.042
μ = 0.0

domain = [x ∈ Interval(0.0, 1.0),
z ∈ Interval(0.0, 1.0),
t ∈ Interval(0.0, 1.0)]

ic_bc = [u(t,0.0,z) ~ 0,
u(t,1.0,z) ~ 0,
w(t,x,0.0) ~ 0,
w(t,x,1.0) ~ 0,
σ₃(t,0.0,z) ~ 0,
σ₃(t,1.0,z) ~ 0,
σ₃(t,x,0.0) ~ 0,
σ₃(t,x,1.0) ~ 0,
T(t,x,0.0) ~ 1,
T(t,x,1.0) ~ θ,
Dx(T(t,0.0,z)) ~ 0,
Dx(T(t,1.0,z)) ~ 0,
u(0,x,z) ~ 0,
w(0,x,z) ~ 0,
P(0,x,z) ~ 0.1,
T(0,x,z) ~ 1 - 0.9*z,
σ₁(0,x,z) ~ 0.0,
σ₃(0,x,z) ~ 0.0,
η(0,x,z) ~ 0.2,
]

eqs = [Dx(u(t,x,z)) + Dz(w(t,x,z)) ~ 0,
Dx(P(t,x,z)) - Dx(σ₁(t,x,z)) - Dz(σ₃(t,x,z)) ~ 0,
Dz(P(t,x,z)) - Dx(σ₃(t,x,z)) + Dz(σ₁(t,x,z)) + Ra*(1-T(t,x,z)) ~ 0,
σ₁(t,x,z) - 2*η(t,x,z)*Dx(u(t,x,z)) ~ 0,
σ₃(t,x,z) - η(t,x,z)*( Dz(u(t,x,z)) + Dx(w(t,x,z)) ) ~ 0,
η(t,x,z) ~ exp( (1 - T(t,x,z) + μ*(1 - z - T(t,x,z))) / (ε*T(t,x,z)) ),
Dt(T(t,x,z)) + u(t,x,z)*Dx(T(t,x,z)) + w(t,x,z)*Dz(T(t,x,z)) +
D*T(t,x,z)*w(t,x,z) ~ Dxx(T(t,x,z)) + Dzz(T(t,x,z)) +
(D/Ra)*(σ₁(t,x,z) + σ₃(t,x,z))^2 / (2*η(t,x,z))
]

@named sys = PDESystem(eqs, ic_bc, domain, [t,x,z], [u(t,x,z), w(t,x,z), P(t,x,z),
T(t,x,z), σ₁(t,x,z), σ₃(t,x,z), η(t,x,z)])

dx = 0.1
dz = 0.1

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

prob = discretize(sys, discretization)


And the when I’m trying to solve it:

@time sol = solve(prob, FBDF(), progress = true, progress_steps = 1, saveat = 0.05)


I get the error:

During the resolution of the non-linear system, the evaluation of the following equation(s) resulted in a non-finite number: [1, 39, 41, 42, 51, 52, 53, 54]


My questions are:

1. What does it mean? How can I tell which equations are the problem according to the numbers in the error?

2. What steps I can take to try and solve this issue?

Thank you!!

Did you try a WENO discretization? This PDE looks a bit tricky so you may need more than just the finite difference discretization.