Method of Lines for simple 2D Boussinesq convection?

Specific Domains Modelling & Simulations
1

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:

@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):

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

2

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 (and corresponding Julia xPU code).

1 Like
3

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!

4

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, also in the intro book by Gerya, and maybe more refs regarding the iterative algorithm used in the code in this GMD paper. Hope this helps.

1 Like