Mixed FEM Spaces / LBB Stable Pairs / Product Spaces
Edit / Summary 08/14/23: Mixed FEM and the associated
combined
space work now as tested with different problems of (Navier-)Stokes type. A tutorial page for the Stokes problem can be found here.
Original Issue: For an important class of PDEs with saddle point structure, including (Navier-)Stokes equations, FEM Spaces with equal interpolation order violate the LBB condition and therefore do not yield stable results with the "canonical" weak forms. The current version of Dune-FEM doesn't seem to handle mixed / combined finite element spaces (at least not on the solver side?) and therefore excludes standard methods to solve flow problems on the C++ side, in particular steady flow problems where standard splitting techniques (to my knowledge) don't apply.
To start a discussion and to make it easier to check the existing code parts, here is a minimal code example using the product space to combine Q2/Q1 ("Taylor-Hood") elements for Velocity/Pressure in a stationary stokes problem with known forcing (Bercovier-Engelman Testcase, Section 3.1 here):
from dune.grid import structuredGrid
from dune.ufl import cell, DirichletBC, BoundaryId
from ufl import div, grad, dot, inner, dx, ds, as_vector, \
SpatialCoordinate, FacetNormal, TrialFunction, TestFunction, eq, conditional
from dune.fem.function import uflFunction
from dune.fem.space import lagrange, product
from dune.fem.scheme import galerkin as solutionScheme
gridView = structuredGrid([0, 0], [1, 1], [50, 50])
dim_space = gridView.dimGrid
spatial_coordinate = SpatialCoordinate(cell(dim_space))
x = spatial_coordinate[0]
y = spatial_coordinate[1]
n = FacetNormal(cell(dim_space))
def weak_form_stationary_stokes(u,v,p,q,f):
"""
Return weak form for the stationary stokes problem, given a set of trial and test basis
functions for velocity and pressure
:param u: (vector-valued) trial function for velocity (in R^d)
:param v: (vector-valued) test function for velocity (in R^d)
:param p: (scalar) trial function for pressure (in R)
:param q: (scalar) test function for pressure (in R)
:param f: (vector valued) external force (uflFunction in R^d)
"""
# Note: This assumes Dirichlet Velocity BCs everywhere:
viscous_stress = inner(grad(u), grad(v))*dx - inner(dot(grad(u),n), v)*ds
pressure_force = -p*div(v)*dx + p*dot(v,n)*ds
continuity = div(u)*q*dx
return viscous_stress + pressure_force + continuity == dot(f, v)*dx
### Analytical Solution for the stationary stokes testcase:
### Bercovier-Engelmann Testcase
# Reference case with analytical solution from https://numtourcfd.pages.math.cnrs.fr/doc/NumtourCFD/v1.0.1-alpha.2/applications/validation/bercovier-engelman.html
# original work: https://doi.org/10.1016/0021-9991(79)90098-6
# benchmark doc: https://doi.org/10.1007/978-3-319-57397-7_12
# (open access: https://www.researchgate.net/publication/318020390_FVCA8_Benchmark_for_the_Stokes_and_Navier-Stokes_Equations_with_the_TrioCFD_Code-Benchmark_Session)
# Compare equations to the setup on the unit square, Section 3.1 of the last reference.
u_exact = uflFunction(gridView,
name="analytical_solution_x",
order=5,
ufl=-256*y*(y-1)*(2*y-1)*x**2*(x-1)**2)
v_exact = uflFunction(gridView,
name="analytical_solution_y",
order=5,
ufl=256*x*(x-1)*(2*x-1)*y**2*(y-1)**2)
p_exact = uflFunction(gridView,
name="analytical_solution_pressure",
order=5,
ufl=(x-0.5)*(y-0.5))
# some random factor
def f_1(x,y):
return 256*(x**2*(x-1)**2*(12*y-6) + y*(y-1)*(2*y-1)*(12*x**2 -12*x + 2))
# Analytical reference, forcing
f_exact = uflFunction(gridView,
name="analytical_solution_force_term",
order=5,
ufl=as_vector([f_1(x,y) + (y-0.5), -f_1(y,x) +(x-0.5)]) )
velocity_order = 2
space_velocity = lagrange(gridView,
dimRange=dim_space,
order=velocity_order,
storage="numpy")
space_pressure = lagrange(gridView,
dimRange=1,
order=velocity_order-1,
storage="numpy")
# Q2-Q1 Polynomial space
taylor_hood_space = product([space_velocity, space_pressure], components=["velocity" , "pressure"])
is_left = eq(BoundaryId(taylor_hood_space),1)
is_right = eq(BoundaryId(taylor_hood_space),2)
is_bot = eq(BoundaryId(taylor_hood_space),3)
is_top = eq(BoundaryId(taylor_hood_space),4)
is_dirichlet = conditional(is_left + is_right + is_bot + is_top > 0.0, 1.0, 0.0)
U = TrialFunction(taylor_hood_space)
V = TestFunction(taylor_hood_space)
# Actually not sure how to do this?
u_h = as_vector([U[0], U[1]])
v_h = as_vector([V[0], V[1]])
p_h = U[2]
q_h = V[2]
dirichlet_bcs_exact = [DirichletBC(taylor_hood_space, [u_exact, v_exact, None], is_left),
DirichletBC(taylor_hood_space, [u_exact, v_exact, None], is_right),
DirichletBC(taylor_hood_space, [u_exact, v_exact, None], is_bot),
DirichletBC(taylor_hood_space, [u_exact, v_exact, None], is_top)]
print("velocity element order: " + str(taylor_hood_space.components[0].order))
print("pressure element order: " + str(taylor_hood_space.components[1].order))
solution_vector = taylor_hood_space.interpolate([0]*(dim_space+1), name="solution_vector")
print("how does this interpolation work?")
print("solution_vector.order = " + str(solution_vector.order))
print("solution_vector[2].order = " + str(solution_vector[2].order))
solver_parameters ={"newton.tolerance": 1e-4,
"newton.verbose": True,
"newton.linear.tolerance": 1e-6,
"newton.linear.preconditioning.method": "sor",
"newton.linear.verbose": False,
"newton.linear.maxiterations": 2000}
scheme_taylor_hood = solutionScheme( [weak_form_stationary_stokes(u_h, v_h, p_h, q_h, f_exact), *dirichlet_bcs_exact],
parameters=solver_parameters,
solver=("istl", "gmres"))
I am not exactly sure if this how things are supposed to look on the python side (design-wise) but just to get a discussion going solutionScheme
for this with an equal-order Lagrangian space with dimRange = dim_space + 1
(velocity components + scalar pressure) but this expectedly fails to converge to a solution. The above example crashes (possibly due to wrong syntax on my side).
Since many FEM spaces are already implemented and the product space
already exists it would be super nice to be able to combine spaces, at least for Lagrangian polynomials in CG / DG setting.