# Method of Lines for simple 2D Boussinesq convection?

**URL:** <https://discourse.julialang.org/t/method-of-lines-for-simple-2d-boussinesq-convection/101416>\
**Category:** Modelling & Simulations\
**Created:** [July 10, 2023, 5:14am UTC](https://discourse.julialang.org/t/method-of-lines-for-simple-2d-boussinesq-convection/101416 "2023-07-10T05:14:38Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![PeX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pex/32/49986_2.png) [@PeX](https://discourse.julialang.org/u/PeX)\
**Post date:** [July 10, 2023, 5:14am UTC](https://discourse.julialang.org/t/method-of-lines-for-simple-2d-boussinesq-convection/101416/1 "2023-07-10T05:14:38Z")

</div>

Hi everyone,  
I’m wondering if the method of lines can be used for a simple 2D convection:

\begin{aligned} \frac{\partial u}{\partial x}+\frac{\partial w}{\partial z} & =0 \\ \frac{\partial p}{\partial x} & =\frac{\partial \sigma\_{11}}{\partial x}+\frac{\partial \sigma\_{13}}{\partial z} \\ \frac{\partial p}{\partial z} & =\frac{\partial \sigma\_{13}}{\partial x}+ \frac{\partial \sigma\_{33}}{\partial z}+Ra(1-T) \\ \sigma\_{11} & =2 \eta \frac{\partial u}{\partial x} \\ \sigma\_{33} & =2 \eta \frac{\partial w}{\partial z} \\ \sigma\_{13} & =\eta\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right) \\ \frac{\partial T}{\partial t} + u\frac{\partial T}{\partial x} + w\frac{\partial T}{\partial z} & = \nabla^2 T \end{aligned}

Where Ra is the Rayleigh number and \eta is the viscosity (constants).

I’m trying to go in a straightforward manner to solve this:

```julia
@parameters t,x,z
@variables u(..), w(..), P(..), T(..), sig_11(..), sig_33(..), sig_13(..);

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

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

Ra = 1e2
η = 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,
		 sig_33(0,x,z) ~ 0,
         sig_11(0,x,z) ~ 0,
         sig_13(0,x,z) ~ 0,
		 T(t,x,0.0) ~ 1,
		 T(t,x,1.0) ~ 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,
		 T(0,x,z) ~ (1.0-z) - 0.01*cos(π*x/1.0)*sin(π*z)
		]

eqs = [Dx(u(t,x,z)) + Dz(w(t,x,z)) ~ 0,
       sig_11(t,x,z) ~ 2*η*Dx(u(t,x,z)),
       sig_33(t,x,z) ~ 2*η*Dz(w(t,x,z)),
       sig_13(t,x,z) ~ η*( Dz(u(t,x,z)) + Dx(w(t,x,z)) ),
       Dx(P(t,x,z)) ~ Dx(sig_11(t,x,z)) + Dz(sig_13(t,x,z)),
       Dz(P(t,x,z)) ~ Dx(sig_13(t,x,z)) + Dz(sig_33(t,x,z)) + Ra*(1-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)) ~ Dxx(T(t,x,z)) + Dzz(T(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), sig_11(t,x,z), sig_33(t,x,z), sig_13(t,x,z)])

dx = 0.1
dz = 0.1

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

prob = discretize(sys, discretization)

@time sol = solve(prob, Rodas5(), progress = true, saveat = 0.05)

```

However, I’m getting the following error (for all variables):

```julia
Warning: Solution has length 1 in dimension t. Interpolation will not be possible for variable T(t, x, z). Solution return code is InitialFailure.

```

I’m not sure how to debug this, one of my guesses was that the solver doesn’t iterate in time, so I changed the initial conditions but I keep getting the same error with every initial condition I’m trying.

Any insights? I’m quite frustrated

---

<div class="post-metadata">

**Author:** ![luraess](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luraess/32/16189_2.png) [@luraess](https://discourse.julialang.org/u/luraess)\
**Post date:** [July 10, 2023, 6:44am UTC](https://discourse.julialang.org/t/method-of-lines-for-simple-2d-boussinesq-convection/101416/2 "2023-07-10T06:44:41Z")

</div>

It may be that the proposed arrangement of the equations may not be appropriate to solve the problem iteratively. One way is to use the first equation to satisfy incompressibility and solve for p, while equations 2-3 are used to retrieve u and v.

Suggesting an alternative way to tackle your problem: check out the [Thermo-Mechanical miniapp in ParallelStencil](https://github.com/omlins/ParallelStencil.jl#thermo-mechanical-convection-2-d-app) (and [corresponding Julia xPU code](https://github.com/omlins/ParallelStencil.jl/blob/main/miniapps/ThermalConvection2D.jl)).

---

<div class="post-metadata">

**Author:** ![PeX](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pex/32/49986_2.png) [@PeX](https://discourse.julialang.org/u/PeX)\
**Post date:** [July 10, 2023, 7:02am UTC](https://discourse.julialang.org/t/method-of-lines-for-simple-2d-boussinesq-convection/101416/3 "2023-07-10T07:02:45Z")

</div>

Wow! The link you shared actually solves the same type of problem! Thank you so much for that! Great reference!  
Can you also expand a little bit more on how to solve for p using incompressibility, to make it possible for an iterative solution? This will help a lot!

---

<div class="post-metadata">

**Author:** ![luraess](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luraess/32/16189_2.png) [@luraess](https://discourse.julialang.org/u/luraess)\
**Post date:** [July 10, 2023, 8:41am UTC](https://discourse.julialang.org/t/method-of-lines-for-simple-2d-boussinesq-convection/101416/4 "2023-07-10T08:41:51Z")

</div>

Welcome. It’s a rather “common” approach in the field of computational geodynamics where one seeks at mostly incompressible solution of viscous Stokes flow. The “split” velocity pressure formulation gives somehow a natural way of solving the system. You may find some useful references in the first part of the following [M2Di paper](https://doi.org/10.1002/2016GC006727), also in the intro [book by Gerya](https://www.cambridge.org/ch/universitypress/subjects/earth-and-environmental-science/structural-geology-tectonics-and-geodynamics/introduction-numerical-geodynamic-modelling-2nd-edition?format=HB&isbn=9781107143142), and maybe more refs regarding the iterative algorithm used in the code in this [GMD paper](https://doi.org/10.5194/gmd-15-5757-2022). Hope this helps.