Hi,
I am trying to visualise some divergence-free Navier-Stokes results, and pretty much all I find in this topic is to calculate any kind of velocity (not at all necessarily div-free), and then as a second step, project it to a div-free space. I find this strange, is there really not an option to create something like
V = VectorFunctionSpace(mesh, “Lagrange”, 2, constraint=divergence-free)?
This would be important because not only the solution function u should be divergence-free (which is achieved e.g. by a projection later), but also the velocity test function v.
How to make sure that both the u Function(V) and the v TestFunction(V) is divergence-free?
Thank you so much!
Usually the divergence-free constraint is enforced in a discrete sense, i.e., with respect to some finite dimensional set of pressure test functions:
(\nabla\cdot\mathbf{u}^h,q^h)_{L^2} = 0~~~\forall q^h\in\mathcal{Q}^h\text{ ,}
where \mathcal{Q}^h is the discrete pressure test space. A velocity satisfying such a condition can be solved for alongside the corresponding pressure by using a mixed element, as in the Stokes demo here:
which can be modified in a straightforward way to include the nonlinear convection term for the Navier–Stokes problem. Some care must be taken to choose the velocity and pressure spaces in a stable way, though. Alternatively, one can include stabilization in the formulation, with the most popular choice being the so-called pressure-stabilizing Petrov–Galerkin (PSPG) method (typically in conjunction with streamline-upwind Petrov–Galerkin (SUPG) when there is convection). With PSPG, you can use whatever pressure and velocity spaces you want, and the resulting equation system tends to be friendlier to iterative linear solvers, but this modifies the approximate div-free condition satisfied by \mathbf{u}^h (in a way that still converges toward a div-free exact solution under refinement).
It is possible to get an exactly (pointwise) div-free \mathbf{u}^h, if \mathcal{V}^h and \mathcal{Q}^h are chosen such that the divergence of every function in \mathcal{V}^h is in \mathcal{Q}^h. This can be done using standard finite elements in FEniCS as outlined here
but is subject to some special requirements on the mesh. Exactly div-free solutions can also be obtained using isogeometric spline spaces, which are accessible in FEniCS using a library that I wrote. In isogeometric analysis, you can also solve directly for a vector potential whose curl is the velocity, such that the velocity is div-free by construction. However, exactly div-free approaches are not currently widely-used.
2 Likes
Thank you so much for your detailed answer, I really appreciate it.
The Stokes demo you mentioned seems to have a 2-dimensional velocity field. I can’t seem to find what is the 3-dimensional version of
“P2 = VectorElement(“Lagrange”, mesh.ufl_cell(), 2)”
I use
mesh = UnitCubeMesh(20, 20, 20)
V = VectorElement(“Lagrange”, ???, degree = 2, dim = 3)
In place of “???” is there something like mesh.ufl3D that I could use? Thank you so much!
If the mesh is 3D, then its corresponding UFL cell, returned by its ufl_cell() method, will be 3D. So, there is no need to change this argument. If the dim argument is not given explicitly, it is automatically chosen to equal the dimension of the cell, so, on a 2D mesh, it will be 2 and on a 3D mesh it will be 3. You only need to set dim=... if you want it to be different than the mesh dimension.
Yes, I am using a 3D UnitCubeMesh. I thought maybe the element dimensions are the problem, but then maybe something else — I’ve tried to update the Stokes demo you mentioned. The solution seems to be “empty” inside. When I visualise it in ParaView, it is defined on the boundary, but blank in the inside. Would you have an idea what I am doing wrong? (the splits are necessary because I need to have the specific u_i components as later I want to add weights - that’s why the very detailed form)
mesh = UnitCubeMesh(20, 20, 20)
V = VectorElement(“Lagrange”, mesh.ufl_cell(), degree = 2, dim = 3)
P = FiniteElement(“Lagrange”, mesh.ufl_cell(), degree = 1)
TH = V * P
VP = FunctionSpace(mesh, TH)
noslipbasin = DirichletBC(VP.sub(0), (0, 0, 0), “on_boundary && x[2] < 1.0 - DOLFIN_EPS”)
zerotop_u = DirichletBC(VP.sub(0).sub(2), 0, “on_boundary && x[2] > 1.0 - DOLFIN_EPS”)
zerotop_p = DirichletBC(VP.sub(1), 0, “on_boundary && x[2] > 1.0 - DOLFIN_EPS”)
bcu = [noslipbasin, zerotop_u, zerotop_p]
# Define variational problem
up_ = Function(VP)
(v, q) = TestFunctions(VP)
(u_, p_) = up_.split(True)
(u1_, u2_, u3_) = u_.split(True)
(v1, v2, v3) = split(v)
F = inner(u_, grad(u1_)) * v1 * dx + inner(u_, grad(u2_)) * v2 * dx + inner(grad(u1_),grad(v1)) * dx + inner(grad(u2_),grad(v2)) * dx + inner(u_, grad(u3_)) * v3 * dx + inner(grad(u3_),grad(v3)) * dx - inner(p_, div(v))*dx + inner(q, div(u_))*dx
solve(F == 0, up_, bcu, solver_parameters = {“newton_solver”:{ “linear_solver” : “gmres”}})
(u,p) = up_.split(True)
print (“Norm of velocity coefficient vector: %.15g” % u.vector().norm(“l2”))
print (“Norm of pressure coefficient vector: %.15g” % p.vector().norm(“l2”))
(u,p) = up_.split()
ufile_pvd = File(“velocity.pvd”)
ufile_pvd << u
pfile_pvd = File(“pressure.pvd”)
pfile_pvd << p
Thank you so much if you can give me a hint.
I just get that the solution is zero everywhere. New Functions are initialized to zero, and this already satisfies the equations and BCs, so Newton’s method exits after zero iterations. Sometimes I’ve seen ParaView do strange things when trying to plot an exactly-uniform function. (I think it may depend on which version you’re using.) Some other notes:
-
I avoid setting Dirichlet BCs for the pressure (although I see the demo does this); it disrupts global mass conservation, since the constant function is no longer in your pressure test space. The one exception I can think of is if you are setting it at just one node, to get rid of an un-determined hydrostatic mode. (Mathematicians prefer to apply a constraint that the average presure is zero, which you can technically do in FEniCS with a
"Real"-type element, but then you have one DoF communicating with the entire domain.) To set the pressure on a boundary, it is better to integrate a weakly-enforced traction BC. -
Applying GMRES with a default preconditioner does not perform well for saddle-point problems (like Galerkin’s method for incompressible Navier–Stokes). It is possible to do something smarter with the solver, e.g., https://fenicsproject.org/docs/dolfin/latest/python/demos/stokes-iterative/demo_stokes-iterative.py.html . However, most of my own experience with Navier–Stokes has involved stabilized methods, which do not result in discrete saddle-point problems, so I could always get away with simple-minded iterative solvers.
1.) Okay, thank you, I have taken it out from the BCs.
2.) I am not sure what to do. The iterative Stokes demo you shared is used for a linear problem with a Krylov solver, but I think in my case I would need to keep the problem nonlinear. This is because the problem I actually want to solve is to visualise the stationary Navier-Stokes solutions for the case when a small epsilon weight is applied to the equation in the 3rd dimension (i.e. to inner(u, grad(u3)) * v3 * dx, etc.) Because of this, it seems to me that linearised approaches such as Picard don’t work because the scheme loses some dependencies as epsilon becomes small, and I need to make sure that the u I solve for is still there. I don’t know if that makes sense. This is why I am trying with this nonlinear approach.
When I modify F to contain the additional boundary traction term, i.e. F is
F = inner(u, grad(u1)) * v1 * dx + inner(u, grad(u2)) * v2 * dx + inner(grad(u1),grad(v1)) * dx + inner(grad(u2),grad(v2)) * dx + inner(u, grad(u3)) * v3 * dx + inner(grad(u3),grad(v3)) * dx - inner(p, div(v))*dx + inner(q, div(u))*dx - inner(theta, v) * ds(1),
then I get no solution.
With gmres I simply get the classic runtime “DOLFIN encountered an error.”, with mumps the process received a “Killed”, and if I don’t specify, I just very soon (within approx. 1-2 mins) get an “out of memory” error.
Alltogether this is what I am trying to do, as seen below. Do you think there is a way that I could move forward with it? Thanks a lot.
mesh = UnitCubeMesh(20, 20, 20)
V = VectorElement(“Lagrange”, mesh.ufl_cell(), degree = 2, dim = 3)
P = FiniteElement(“Lagrange”, mesh.ufl_cell(), degree = 1)
TH = V * P
VP = FunctionSpace(mesh, TH)
# Define boundary conditions
class UpperBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[2], 1.0)
upperboundary = UpperBoundary()
boundaries = MeshFunction(“size_t”, mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
upperboundary.mark(boundaries, 1)
ds = Measure(‘ds’)[boundaries]
wind_shear_x = 10.0
wind_shear_y = 10.0
theta = Constant((wind_shear_x, wind_shear_y, 0))
noslipbasin = DirichletBC(VP.sub(0), (0, 0, 0), “on_boundary && x[2] < 1.0 - DOLFIN_EPS”)
zerotop_u = DirichletBC(VP.sub(0).sub(2), 0, “on_boundary && x[2] > 1.0 - DOLFIN_EPS”)
bcu = [noslipbasin, zerotop_u]
# Define variational problem
up = TrialFunction(VP)
u,p = split(up)
u1, u2, u3 = split(u)
(v, q) = TestFunctions(VP)
v1, v2, v3 = split(v)
up_ = Function(VP)
(u_, p_) = split(up_)
(u1_, u2_, u3_) = split(u_)
eps = 1.0
F = inner(u, grad(u1)) * v1 * dx + inner(u, grad(u2)) * v2 * dx + inner(grad(u1),grad(v1)) * dx + inner(grad(u2),grad(v2)) * dx + eps * inner(u, grad(u3)) * v3 * dx + eps * inner(grad(u3),grad(v3)) * dx - inner(p, div(v))*dx + inner(q, div(u))*dx - inner(theta, v) * ds(1)
F = action(F, up_)
J = derivative(F, up_, up)
problem = NonlinearVariationalProblem(F, up_, bcu, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
solver.solve()
#solve(F == 0, up_, bcu, J, solver_parameters = {“newton_solver”:{ “linear_solver” : “gmres”}})
(u,p) = up_.split(True)
print (“Norm of velocity coefficient vector: %.15g” % u.vector().norm(“l2”))
print (“Norm of pressure coefficient vector: %.15g” % p.vector().norm(“l2”))
(u,p) = up_.split(True)
ufile_pvd = File(“velocity.pvd”)
ufile_pvd << u
pfile_pvd = File(“pressure.pvd”)
pfile_pvd << p
You could try applying the preconditioning strategy to the linear solve within each Newton iteration. There is probably a smart way to do this with the built-in nonlinear solvers, but, in worst case, you can always just program the Newton iteration manually, since it’s not too complicated.
Alternatively, you can use a pressure-stabilizing approach where the PSPG term will fill in the pressure–pressure zero block, which tends to perform okay on mid-size 3D problems with GMRES and a simple preconditioner, like Jacobi. However, you’d have to be a bit careful about this if you’re modifying the PDE; the stabilization terms would need to be modified in a consistent way.
1 Like