gistfile1.txt

from firedrake import *
 
m = PeriodicUnitIntervalMesh(5)
mesh = ExtrudedMesh(m, layers=5, layer_height=1.0/5)
 
#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)
 
#Mixed function space
 
W = MixedFunctionSpace( (V1,V2) )
 
state = Function(W)
u, p = state.split()
 
f = project(Expression("sin(pi*x[0])"),W[1])
 
 
w, phi = TestFunctions(W)
 
u, p = split(state)
 
#Equations
F = (inner(w,u) - (div(w)*p) +
phi*div(u) +
phi*f
)*dx
 
#Boundary conditions
bc1 = [DirichletBC(W.sub(0), Expression(("0.", "0.")), x)
       for x in ["top", "bottom"]]
 
#Nullspace

null_vec = Function(W)
_, null_p = null_vec.split()
null_p.assign(1/sqrt(null_p.function_space().dof_count))

nullspace = VectorSpaceBasis(vecs=[null_vec])
 
solver_parameters={'snes_monitor': True,
                   'snes_view': True,
                   'ksp_monitor_true_residual': True,
                   'ksp_max_it': 10,
                   'ksp_view': True,
                   'ksp_type': 'gmres',
                   'snes_converged_reason': True,
                   'ksp_converged_reason': True}

a = derivative(F, state)
L = assemble(a, bcs=bc1)
dV0 = V1.dof_count
dV1 = V2.dof_count

A = np.zeros((dV0+dV1, dV0+dV1))
A[:dV0, :dV0] = L.M[0, 0].values
A[:dV0, dV0:dV0+dV1] = L.M[0, 1].values
A[dV0:dV0+dV1, :dV0] = L.M[1, 0].values
A[dV0:dV0+dV1, dV0:dV0+dV1] = L.M[1, 1].values

u, s, v = np.linalg.svd(A)

rhs = assemble(F)
for bc in bc1:
    bc.zero(rhs)

rhs2 = nullspace.orthogonalize(rhs)
solve(L, state, rhs2, bcs=bc1, solver_parameters=solver_parameters)

print state.dat.data

state.assign(0)
solve(L, state, rhs, bcs=bc1, nullspace=nullspace, solver_parameters=solver_parameters)
print state.dat.data