Hello everyone,
I am working on solving the incompressible Navier-Stokes equations for the benchmark problem of flow around a cylinder using monolithic mixed finite elements in FEniCSx. Previously, I successfully implemented both a nonlinear solver using Newton’s method and a linear solver (using an approximation of the convection term) in FEniCS, and everything worked perfectly.
After transitioning to FEniCSx, I first re-implemented the nonlinear solver, which works well. However, I am now trying to solve the problem as a linear one by approximating the nonlinear convection term. The issue I’m facing is that while the code converges for a specific mesh and function space degree, it starts to diverge (producing inf values) when I refine the mesh or change the degree of the function spaces.
I am unsure of what might be causing this instability. Has anyone encountered similar issues or could you guide how to stabilize the solver with mesh refinement or function space degree changes?
Here is the code I’m using for reference:
import gmsh
import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import time
import tqdm.autonotebook
from mpi4py import MPI
from petsc4py import PETSc
from basix.ufl import element
from dolfinx.cpp.mesh import to_type, cell_entity_type
from dolfinx.fem import (Constant, Function, functionspace,
assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,create_vector,
create_matrix, set_bc, NonlinearProblem)
from dolfinx.graph import adjacencylist
from dolfinx.geometry import bb_tree, compute_collisions_points, compute_colliding_cells
from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio, XDMFFile)
from dolfinx.mesh import (create_mesh, meshtags_from_entities, CellType, compute_midpoints,
create_unit_cube, create_unit_square, meshtags, locate_entities_boundary)
from dolfinx.nls.petsc import NewtonSolver
from ufl import (FacetNormal, Identity, Measure, TestFunction, TrialFunction, TestFunctions, TrialFunctions,
as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym, system, split)
import ufl
import pyvista
import basix
import dolfinx.plot as plot
from dolfinx import io, nls, log
from dolfinx.cpp.log import set_log_level, LogLevel
#set_log_level(LogLevel.INFO)
# ********* Gmsh part *********
gmsh.initialize()
L = 2.2
H = 0.41
c_x = c_y = 0.2
r = 0.05
gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0
if mesh_comm.rank == model_rank:
rectangle = gmsh.model.occ.addRectangle(0, 0, 0, L, H, tag=1)
obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)
"""
The next step is to subtract the obstacle from the channel, such that we do not mesh the interior of the circle.
"""
if mesh_comm.rank == model_rank:
fluid = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle)])
gmsh.model.occ.synchronize()
"""
To get GMSH to mesh the fluid, we add a physical volume marker
"""
fluid_marker = 1
if mesh_comm.rank == model_rank:
volumes = gmsh.model.getEntities(dim=gdim)
assert (len(volumes) == 1)
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")
"""
To tag the different surfaces of the mesh, we tag the inflow (left hand side) with marker 2, the outflow (right hand side) with marker 3 and the fluid walls with 4 and obstacle with 5. We will do this by computing the center of mass for each geometrical entity.
"""
inlet_marker, outlet_marker, wall_marker, obstacle_marker = 2, 3, 4, 5
inflow, outflow, walls, obstacle = [], [], [], []
if mesh_comm.rank == model_rank:
boundaries = gmsh.model.getBoundary(volumes, oriented=False)
for boundary in boundaries:
center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
if np.allclose(center_of_mass, [0, H / 2, 0]):
inflow.append(boundary[1])
elif np.allclose(center_of_mass, [L, H / 2, 0]):
outflow.append(boundary[1])
elif np.allclose(center_of_mass, [L / 2, H, 0]) or np.allclose(center_of_mass, [L / 2, 0, 0]):
walls.append(boundary[1])
else:
obstacle.append(boundary[1])
gmsh.model.addPhysicalGroup(1, walls, wall_marker)
gmsh.model.setPhysicalName(1, wall_marker, "Walls")
gmsh.model.addPhysicalGroup(1, inflow, inlet_marker)
gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
gmsh.model.addPhysicalGroup(1, outflow, outlet_marker)
gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
gmsh.model.addPhysicalGroup(1, obstacle, obstacle_marker)
gmsh.model.setPhysicalName(1, obstacle_marker, "Obstacle")
"""
In our previous meshes, we have used uniform mesh sizes. In this example, we will have variable mesh sizes to resolve the flow solution in the area of interest; close to the circular obstacle. To do this, we use GMSH Fields.
"""
# Create distance field from obstacle.
# Add threshold of mesh sizes based on the distance field
# LcMax - /--------
# /
# LcMin -o---------/
# | | |
# Point DistMin DistMax
res_min = r / 7
if mesh_comm.rank == model_rank:
distance_field = gmsh.model.mesh.field.add("Distance")
gmsh.model.mesh.field.setNumbers(distance_field, "EdgesList", obstacle)
threshold_field = gmsh.model.mesh.field.add("Threshold")
gmsh.model.mesh.field.setNumber(threshold_field, "IField", distance_field)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMin", res_min)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMax", 0.15 * H)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMin", r)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMax", 2.0 * H)
min_field = gmsh.model.mesh.field.add("Min")
gmsh.model.mesh.field.setNumbers(min_field, "FieldsList", [threshold_field])
gmsh.model.mesh.field.setAsBackgroundMesh(min_field)
"""
Generating the mesh
We are now ready to generate the mesh. However, we have to decide if our mesh should consist of triangles or quadrilaterals. In this demo, to match the DFG 2D-3 benchmark, we use second order quadrilateral elements.
"""
if mesh_comm.rank == model_rank:
gmsh.option.setNumber("Mesh.Algorithm", 5)
#gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2) #
#gmsh.option.setNumber("Mesh.RecombineAll", 1) # Uncomment if you want a quadrilateral mesh
#gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1) #
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize("Netgen")
#
# ********* end of the Gmsh part *********
#
mesh, _, ft = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
ft.name = "Facet markers"
t = 0.0
T = 0.1 # Final time
dt = 1 / 100 # Time step size
num_steps = int(T / dt)
k = Constant(mesh, PETSc.ScalarType(dt))
nu = Constant(mesh, PETSc.ScalarType(0.001)) # Dynamic viscosity
rho = Constant(mesh, PETSc.ScalarType(1)) # Density
n = FacetNormal(mesh) # Normal pointing out of obstacle
Ve = element("CG", mesh.topology.cell_name(), 2, shape=(mesh.topology.dim, ))
Pe = element("CG", mesh.topology.cell_name(), 1)
We = basix.ufl.mixed_element([Ve,Pe])
W = functionspace(mesh, We)
V_, _ = W.sub(0).collapse()
Q_, _ = W.sub(1).collapse()
w00 = Function(W)
nx = len(w00.x.array)
print('DoFs on this level:', nx)
#-------------------------------
# Define the inlet velocity
#------------------------------
#
def InletVelocity(t):
g = Function(V_)
g.interpolate(
lambda x : np.vstack((
4 * 1.5 * np.sin(t * np.pi / 8) * x[1] * (0.41 - x[1]) / (0.41**2),
np.full(x.shape[1], 0.0, dtype=np.float64))))
return g
#
# --------------------------------
# Define the boundary conditions
#---------------------------------
u_in = Function(V_)
u_zero = Function(V_)
p_zero = Function(Q_)
#
u_in.interpolate(InletVelocity(t))
u_zero.interpolate(lambda x : np.vstack((np.full(x.shape[1], 0.0, dtype=np.float64), np.full(x.shape[1], 0.0, dtype=np.float64))))
p_zero.interpolate(lambda x : np.full(x.shape[1], 0.0, dtype=np.float64))
#
fdim = mesh.topology.dim - 1
#
inlet = locate_entities_boundary(mesh, fdim, lambda x: np.isclose(x[0], 0.0))
outlet = locate_entities_boundary(mesh, fdim, lambda x: np.isclose(x[0], 2.2))
walls = locate_entities_boundary(mesh, fdim, lambda x: np.isclose(x[1], 0.0) | np.isclose(x[1], 0.41))
cylin = locate_entities_boundary(mesh, fdim, lambda x: np.isclose((x[0]-0.2)**2 + (x[1]-0.2)**2, 0.05**2))
#
inlet_bc = dirichletbc( u_in, locate_dofs_topological((W.sub(0), V_), fdim, inlet), W.sub(0))
Walls_bc = dirichletbc( u_zero, locate_dofs_topological((W.sub(0), V_), fdim, walls), W.sub(0))
cylin_bc = dirichletbc( u_zero, locate_dofs_topological((W.sub(0), V_), fdim, cylin), W.sub(0))
outlet_bc = dirichletbc( p_zero, locate_dofs_topological((W.sub(1), Q_), fdim, outlet), W.sub(1))
#
bcu = [inlet_bc, Walls_bc, cylin_bc, outlet_bc]
#
#----------------------------------
# Define trial and test functions
#----------------------------------
#
u, p = TrialFunctions(W)
v, q = TestFunctions(W)
#
w = Function(W)
#
w_old = Function(W)
u_old, p_old = split(w_old)
#
# Weak formulation
#
F = ufl.dot((u - u_old)/k, v)*ufl.dx
F += ufl.dot(ufl.dot(u_old, ufl.nabla_grad(u)), v) * ufl.dx
F += ufl.dot(ufl.dot(u, ufl.nabla_grad(u_old)), v) * ufl.dx
F -= ufl.dot( ufl.dot(u_old, ufl.nabla_grad(u_old)), v) * ufl.dx
F += ufl.inner(nu*ufl.nabla_grad(u), ufl.nabla_grad(v))*ufl.dx
F -= ufl.dot(p,ufl.div(v))*ufl.dx
F -= ufl.dot(ufl.div(u),q)*ufl.dx
a = form(lhs(F))
L = form(rhs(F))
A = create_matrix(a)
b = create_vector(L)
# set the solver
solver = PETSc.KSP().create(mesh.comm)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
cells, types, x = plot.vtk_mesh(V_)
gridu = pyvista.UnstructuredGrid(cells, types, x)
# Create plotter
plotter = pyvista.Plotter()
plotter.show(interactive_update=True)
from pathlib import Path
folder = Path("results")
folder.mkdir(exist_ok=True, parents=True)
#
progress = tqdm.autonotebook.tqdm(desc="Solving PDE", total=num_steps)
#
for i in range(num_steps):
progress.update(1)
# Update current time step
t += dt
#
# Update inlet velocity
u_in.interpolate(InletVelocity(t))
#
# solve the system
A.zeroEntries()
assemble_matrix(A, a, bcs=bcu)
A.assemble()
with b.localForm() as loc:
loc.set(0)
assemble_vector(b, L)
#
apply_lifting(b, [a], [bcu])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
#
set_bc(b, bcu)
solver.setOperators(A)
solver.solve(b, w.vector)
w.x.scatter_forward()
# Update variable with solution form this time step
with w.vector.localForm() as loc_w, w_old.vector.localForm() as loc_w_old:
loc_w.copy(loc_w_old)
#
# Plot the solution
uw, pw = w_old.split()
u_ = uw.collapse()
p_ = pw.collapse()
#
plotter.clear()
#
values = np.zeros((x.shape[0], 3), dtype=np.float64)
values[:, :len(u_)] = u_.x.array.real.reshape((x.shape[0], len(u_)))
gridu["u"] = values
glyphs = gridu.glyph(orient="u", factor=0.08)
#c
plotter.add_mesh(glyphs, cmap='jet')
plotter.add_bounding_box(line_width=0.5, color='black')
plotter.update()
plotter.view_xy()
plotter.show()
Any advice would be greatly appreciated! Thanks in advance.
Best regards,
Abdelouahed