gistfile1.py

"""
Solve the linearised Boussinesq equations in a vertical slice using
an extruded mesh.
"""
from firedrake import *

m = UnitIntervalMesh(5)
mesh = ExtrudedMesh(m, layers=4, layer_height=0.25)

#1D spaces in horizontal
IUh0 = FiniteElement("CG", "interval", 1)
IUh1 = FiniteElement("DG", "interval", 0)

#1D spaces in vertical
IUv0 = FiniteElement("CG", "interval", 1)
IUv1 = FiniteElement("DG", "interval", 0)

##Setting up the finite element spaces

#Pressure space
V2_elt = OuterProductElement(IUh1, IUv1)
V2 = FunctionSpace(mesh, V2_elt)

#Velocity space
V1_v = OuterProductElement(IUh1, IUv0)
V1_h = OuterProductElement(IUh0, IUv1)
V1_elt = HDiv(V1_v) + HDiv(V1_h)
V1 = FunctionSpace(mesh,V1_elt)

#Temperature space
Vt = FunctionSpace(mesh,V1_v)

#Mixed function space

W = MixedFunctionSpace( (V1,V2,Vt) )

state = Function(W)
oldstate = Function(W)
u, p, theta = state.split()

u0 = project(Expression(("0.0","0.0")),W[0])
p0 = project(Expression("0.0"),W[1])
theta0 = project(Expression("sin(pi*x[0])*sin(pi*x[1])"),W[2])
oldu, oldp, oldtheta = oldstate.split()
oldu.assign(u0)
oldp.assign(p0)
oldtheta.assign(theta0)

w, phi, gamma = TestFunctions(W)

oldu, oldp, oldtheta = split(oldstate)
u, p, theta = split(state)
dt = 0.01
alpha = 0.5
Nsq = 1.
gammag = 1.0
Theta_z = 1.0

bbar = alpha*p + (1-alpha)*oldp
thetabar = alpha*theta + (1-alpha)*oldtheta
ubar = alpha*u + (1-alpha)*oldu

F = (inner(w,u-oldu) - dt*(div(w)*p + w[1]*gammag*thetabar) +
     phi*div(u) +
     gamma*(theta-oldtheta + dt*Theta_z*ubar[1])
    )*dx

t = 0.0
end = 0.5

file = File("temp.pvd")
while (t <= end):
    solve(F == 0, state)
    oldstate.assign(state)

    _, _, theta = state.split()
    file << theta
    t += dt