gistfile1.py

from firedrake import *

dimension = 1
mesh = UnitIntervalMesh(10)
U = FunctionSpace(mesh, "CG", 2)
H = FunctionSpace(mesh, "CG", 1)
W = U * H # Mixed function space

q = Function(W)
# Define trial and test functions
u, h = split(q)
w, v = TestFunctions(W)

q_old = Function(W).interpolate(Expression(("2.0", "8.0")))
u_old, h_old = split(q_old)

F = 0

dt = 0.5
M_momentum = (1.0/dt)*(inner(w, u) - inner(w, u_old))*dx
F += M_momentum
 
M_continuity = (1.0/dt)*(inner(v, h) - inner(v, h_old))*dx
F += M_continuity
          
# Non-linear term
A_momentum = 0
for dim_j in range(dimension):
   A_momentum += inner(dot(grad(u)[dim_j], u), w)*dx
F += A_momentum

C_momentum = 0
for dim in range(dimension):
   C_momentum += -inner(w, grad(h)[dim])*dx
F -= C_momentum

divergence = 0
for dim in range(dimension):
   divergence += grad(h*u)[dim]
Ct_continuity = inner(v, divergence)*dx
F += Ct_continuity

print solve(F == 0, q, solver_parameters={'ksp_monitor': True, 'ksp_view': False, 'pc_view': False, 'ksp_type':'gmres', 'pc_type':'jacobi'})
fields = q.split()
print len(fields)