I am having trouble solving incompressible Navier-Stokes with a scheme I will detail below. The mesh generation is following Dokken’s FEniCSx tutorial for the flow past a cylinder (Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark) — FEniCSx tutorial).
Here are the equations governing my problem, first one is the momentum equation, second is the incompressibility constraint:
I am employing a BDF (Backward Differentiation formula) of order 1 or 2 to approximate the time derivative (equation (3)), relaxing the incompressibility constraint (equation (4)) and using a Newton approximation for the convective term (equation (5)):
With that in mind, we can derive the weak formulation (after some tensor calculations) and obtain equation (6):
where
\mathbf\nabla,\;\mathbf\nabla\cdot,\;\overline\otimes,\;\overline{\overline\otimes},\;\Omega,\;\mathbf u,\;p,\;\mathbf v,\;q,\;\mathbf f,\;n,\;k,\;\epsilon,Re
are respectively gradient, divergence operators, simple and double dot products, the space we integrate on, the velocity vector (2D), pressure, test functions for velocity and pressure, volume forces, step, order of the BDF (Backward Differential Formula, defining the time derivative of velocity), a relaxation constant (used to relax the incompressibility constraint) and the Reynolds number (equations are not dimensional).
In my script, I use first order BDF for the first time step then the second order BDF, a Reynolds number of 500, a small relaxation constant and no volumic forces such that:
Re=500, k\in\lbrace1,2\rbrace, \epsilon=10^{-6}, \mathbf f=\mathbf0
Let’s now describe the boundary conditions, they are almost the same as in the tutorial, I just modified the inlet profile so that it does not depend on time (equation (7)):
- null velocity on walls and cylinder
- null pressure on the outlet
- inlet velocity profile following
where H,\;L,\;U_0,\;\partial S are respectivly the heigth and length of the rectangle, the lengthwise/longitudinal velocity at the middle of the inlet and the circle delimiting the cylinder (the obstacle).
Now I have to explain what I do in the code. The error occurs at first step, so we only consider first order BDF as it is the one used to calculate the solution at step 1. Let’s go:
First I import everything and create the mesh like in the tutorial (I modified the sizes of the rectangle and cylinder, and modified some other stuff, but I am very confident the issue does not lie here)
############################################
#--------------- Imports ------------------#
############################################
import gmsh
import os
import numpy as np
import matplotlib.pyplot as plt
import tqdm.autonotebook
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.cpp.mesh import to_type, cell_entity_type
from dolfinx.fem import (Constant, Function, FunctionSpace, VectorFunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_topological, set_bc, Expression)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector, create_vector, create_matrix, set_bc, LinearProblem, NonlinearProblem)
from dolfinx.graph import create_adjacencylist
from dolfinx.geometry import BoundingBoxTree, compute_collisions, compute_colliding_cells
from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio)
from dolfinx.mesh import create_mesh, meshtags_from_entities
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TestFunctions, TrialFunction, VectorElement, MixedElement, as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym, curl, derivative, split)
from dolfinx.io import XDMFFile
from dolfinx.nls.petsc import NewtonSolver
from os import path, mkdir
from dolfinx import log
mesh_comm = MPI.COMM_WORLD
############################################
#--------------- Mesh creation ------------#
############################################
model_rank = 0
gmsh.initialize()
# Dimensions
L = 1.0 # Length of the rectangle
H = 0.6 # Height of the rectangle
c_x = c_y = 0.3 # coordinates of the center of the cylinder (also refered a the obstacle)
r = 0.05 # radius of the cylinder (obstacle)
gdim = 2
model_rank = 0
## Creation of the geometry entities (a cylinder within a rectangle)
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 substract 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")
## Here 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.
res_factor = 2
res_min = r / res_factor
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.5 * H)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMin", r)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMax", 2 * 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. Here we use second order quadrilateral elements.
if mesh_comm.rank == model_rank:
gmsh.option.setNumber("Mesh.Algorithm", 8)
gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
gmsh.option.setNumber("Mesh.RecombineAll", 1)
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize("Netgen")
## As we have generated the mesh, we now need to load the mesh and corresponding
## facet markers into DOLFINx. See https://jsdokken.com/src/tutorial_gmsh.html for explanations
mesh, c, ft = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
mesh.name = "Cylinder"
c.name = "Cell markers"
ft.name = "Facet markers"
Then I define some problem specific parameters and the boundary conditions. Those are exactly the same as in the tutorial. I use a mixed function space.
############################################
#--- Parameters & boudnary conditions -----#
############################################
## We define our problem specific parameters
dt = Constant(mesh, PETSc.ScalarType(1/1600))
nu = Constant(mesh, PETSc.ScalarType(1/500)) # nu = 1/Re
relax = Constant(mesh, PETSc.ScalarType(1e-6)) # relaxation constant
U_0 = 2 # longitunal/lenghtwise velocity at coordinate (x,y)=(0,H/2)
# Inlet velocity profile
class Inlet_Velocity():
def __init__(self, U_0, H):
self.U_0 = U_0
self.H = H
def __call__(self, x):
values = np.zeros((gdim, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = 4 * self.U_0 * x[1] * (self.H - x[1])/(self.H**2)
return values
## As we have created the mesh and relevant mesh tags, we can now specify the function
## spaces V aand Q along with the boundary conditions. As the ft contains markers for facets,
## we use this class to find the facets for the inlet and walls.
Velocity_element = VectorElement("Lagrange", mesh.ufl_cell(), 2)
Pressure_element = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement([Velocity_element, Pressure_element]))
V, V_to_W = W.sub(0).collapse()
Q, Q_to_W = W.sub(1).collapse()
fdim = mesh.topology.dim - 1
# Inlet
u_inlet = Function(V)
inlet_velocity = Inlet_Velocity(U_0, H)
u_inlet.interpolate(inlet_velocity)
bcu_inflow = dirichletbc(u_inlet, locate_dofs_topological(V, fdim, ft.find(inlet_marker)))
# Walls
u_nonslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bcu_walls = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(wall_marker)), V)
# Obstacle
bcu_obstacle = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(obstacle_marker)), V)
# grouping all boundary conditions involving the velocity u
bcu = [bcu_inflow, bcu_walls, bcu_obstacle]
# Outlet
p_outlet = PETSc.ScalarType(0)
bcp_outlet = dirichletbc(p_outlet, locate_dofs_topological(Q, fdim, ft.find(outlet_marker)), Q)
# grouping all boundary conditions involving the pressure p
bcp = [bcp_outlet]
bcs = bcu + bcp
Now I define the scheme as in equation (6) with BDF of order 1 (BDF 1). w_prev1 is the solution at step 0 (known) and w is the unknown of the problem (solution at step 1).
############################################
#--------- Discretization scheme ----------#
############################################
# defining the temporal scheme for the first time step (BDF1 is used, so BDF is set to 1)
v, q = TestFunctions(W)
w = Function(W)
u, p = split(w)
w_prev1 = Function(W)
u_prev1, _ = split(w_prev1)
# discretization scheme
F = dot(u - u_prev1, v) * dx\
+ dot(dot(grad(u), u_prev1) + dot(grad(u_prev1), u) - dot(grad(u_prev1), u_prev1), v) * dt * dx\
+ nu * inner(grad(u), grad(v)) * dt * dx\
- p * div(v) * dt * dx\
- q * div(u) * dt * dx\
- relax * p * q * dt * dx
Finally I solve the system for only one time step.
############################################
#----------------- Solving ----------------#
############################################
# creating the problem and solving it
J = derivative(F, w)
problem = NonlinearProblem(F, w, bcs=bcu+bcp, J=J)
solver = NewtonSolver(mesh_comm, problem)
solver.rtol = 1e-6
solver.report = True
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(w)
assert(converged)
print(f"Number of interations: {n:d}")
I would like to have a result like in Implementation — FEniCSx tutorial. Instead I have a core dump error, or this kind of result (no iteration at all):
2023-01-25 08:59:10.937 ( 0.165s) [main ] NewtonSolver.cpp:36 INFO| Newton iteration 0: r (abs) = 0 (tol = 1e-10) r (rel) = -nan(tol = 1e-06)
2023-01-25 08:59:10.937 ( 0.165s) [main ] NewtonSolver.cpp:255 INFO| Newton solver finished in 0 iterations and 0 linear solver iterations.
Number of interations: 0
2023-01-25 08:59:10.992 ( 0.220s) [main ] loguru.cpp:526 INFO| atexit
I am to some extent a beginner in FEniCSx and solvers in general. Maybe I am not using the right solver because, to me, the weak formulation I use is linear as the convective term was linearized (am I right ?). In fact this way of solving incompressible Navier-Stokes equations has already been used in old FEniCS using NonLinearVariationalProblem and NonLinearVariationalSolver as follow:
J = derivative(F_BDF, w)
problem = NonlinearVariationalProblem(F_BDF, w, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.parameters["newton_solver"]["linear_solver"] = "mumps"
solver.solve()
And it worked, but now with FEniCSx it throws error. Can anyone help me findout why is the solver not converging in my case ? I hope I am clear in my explanations.
Thanks and regards.
