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:
\begin{aligned}
& w=0, \quad \tau_3=0, \quad T=1 \quad \text { on } \quad z=0, \\
& w=0, \quad \tau_3=0, \quad T=0.1 \quad \text { on } \quad z=1, \\
& u=0, \quad \tau_3=0, \quad \frac{\partial T}{\partial x}=0 \quad \text { on } \quad x=0,1 .
\end{aligned}
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:
-
What does it mean? How can I tell which equations are the problem according to the numbers in the error?
-
What steps I can take to try and solve this issue?
Thank you!!