# Problem to solve a Navier-Stokes problem with a Custom Newton solver in FEniCSx

Hi everone!
I’m trying to solve a flow around a cylinder with FEniCSx.
The boundaty conditions:
u=(1,0) at inflow and walls
u=(0,0) around the cylinder
p=0 at outflow

This is my code:

import dolfinx
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import math
import ufl
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement,
as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym, MixedElement)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
create_vector, create_matrix, set_bc)

from dolfinx.fem import (Constant, Function, FunctionSpace,
assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc)

##########################################################################################################3

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "Mesh/mesh.xdmf", "r") as xdmf:

########################################################################################################

# Definition of Spaces. (dim and polinomial type)

v_cg1 = VectorElement("CG", domain.ufl_cell(), 2)
s_cg1 = FiniteElement("CG", domain.ufl_cell(), 1)
V = dolfinx.fem.FunctionSpace(domain, v_cg1)
Q = dolfinx.fem.FunctionSpace(domain, s_cg1)
Z = MixedElement([v_cg1, s_cg1])
W = dolfinx.fem.FunctionSpace(domain, Z)

# Boundary subdomains and coditions

########################################################################################################

# Definition of subdomains

x_1 = -40
x_2 = 50
y_1 = -40
y_2 = 40

def upstream(x):
return np.isclose(x[0],x_1)

def downstream(x):
return np.isclose(x[0],x_2)

def walls(x):
return np.logical_or(np.isclose(x[1],y_1),np.isclose(x[1],y_2))

def cylinder(x):
return np.logical_and(np.greater(x[0]**2+x[1]**2, 0.2),np.less(x[0]**2+x[1]**2, 1))

fdim = domain.topology.dim -1

up_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, upstream)
down_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, downstream)
wall_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, walls)
cylinder_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, cylinder)

def u_exact(x):
values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = 1
return values

u_icondition = Function(V)
p_icondition = Function(Q)
u_cylindercondition = Function(V)
u_icondition.interpolate(u_exact)

up_f = dolfinx.fem.locate_dofs_topological(W.sub(0), fdim, up_facets)
wall_f = dolfinx.fem.locate_dofs_topological(W.sub(0), fdim, wall_facets)
cylinder_f = dolfinx.fem.locate_dofs_topological(W.sub(0), fdim, cylinder_facets)
down_f = dolfinx.fem.locate_dofs_topological(W.sub(1), fdim, down_facets)########################################################################################################

# Definition of Spaces. (dim and polinomial type)

v_cg1 = VectorElement("CG", domain.ufl_cell(), 2)
s_cg1 = FiniteElement("CG", domain.ufl_cell(), 1)
V = dolfinx.fem.FunctionSpace(domain, v_cg1)
Q = dolfinx.fem.FunctionSpace(domain, s_cg1)
Z = MixedElement([v_cg1, s_cg1])
W = dolfinx.fem.FunctionSpace(domain, Z)

# Boundary subdomains and coditions

########################################################################################################

# Definition of subdomains

x_1 = -40
x_2 = 50
y_1 = -40
y_2 = 40

def upstream(x):
return np.isclose(x[0],x_1)

def downstream(x):
return np.isclose(x[0],x_2)

def walls(x):
return np.logical_or(np.isclose(x[1],y_1),np.isclose(x[1],y_2))

def cylinder(x):
return np.logical_and(np.greater(x[0]**2+x[1]**2, 0.2),np.less(x[0]**2+x[1]**2, 1))

fdim = domain.topology.dim -1

up_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, upstream)
down_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, downstream)
wall_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, walls)
cylinder_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, cylinder)

def u_exact(x):
values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = 1
return values

u_icondition = Function(V)
p_icondition = Function(Q)
u_cylindercondition = Function(V)
u_icondition.interpolate(u_exact)

up_f = dolfinx.fem.locate_dofs_topological(W.sub(0), fdim, up_facets)
wall_f = dolfinx.fem.locate_dofs_topological(W.sub(0), fdim, wall_facets)
cylinder_f = dolfinx.fem.locate_dofs_topological(W.sub(0), fdim, cylinder_facets)
down_f = dolfinx.fem.locate_dofs_topological(W.sub(1), fdim, down_facets)

#######################################################################################################################

# Bounday conditions

bc_u = dirichletbc(u_icondition, up_f)
bc_w = dirichletbc(u_icondition, wall_f)
bc_cylinder = dirichletbc(u_cylindercondition, cylinder_f)
bc_p = dirichletbc(p_icondition, down_f)

bc = [bc_u, bc_w, bc_cylinder, bc_p]

mu = Constant(domain, PETSc.ScalarType(0.02))
n = FacetNormal(domain)
def epsilon(u):

# Define stress tensor
def sigma(u, p):
return 2*mu*epsilon(u) - p*Identity(len(u))

w = Function(W)
u, p = ufl.split(w)

m = TestFunction(W)

v, q = ufl.split(m)

F = dot(dot(u,nabla_grad(u)),v)*dx + div(u)*q*dx + inner(sigma(u,p),epsilon(v))*dx + dot(p*n, v)*ds - dot(mu*nabla_grad(u)*n, v)*ds

J = ufl.derivative(F, w)
residual = dolfinx.fem.form(F)
jacobian = dolfinx.fem.form(J)

deltax = dolfinx.fem.Function(Wsssssssssssssssssssssssssssssssssssssssssssssssssssss)
A = dolfinx.fem.petsc.create_matrix(jacobian)
L = dolfinx.fem.petsc.create_vector(residual)
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)

i = 0
dx_norm = []
max_iterations = 25
while i < max_iterations:
# Assemble Jacobian and residual
with L.localForm() as loc_L:
loc_L.set(0)
A.zeroEntries()
dolfinx.fem.petsc.assemble_matrix(A, jacobian, bcs=bc)
A.assemble()
dolfinx.fem.petsc.assemble_vector(L, residual)
L.scale(-1)

# Compute b - J(u_D-u_(i-1))
dolfinx.fem.petsc.apply_lifting(L, [jacobian], [bc], x0=[w.vector], scale=1)
# Set dx|_bc = u_{i-1}-u_D
dolfinx.fem.petsc.set_bc(L, bc, w.vector, 1.0)

# Solve linear problem
solver.solve(L, deltax.vector)
deltax.x.scatter_forward()

# Update u_{i+1} = u_i + delta x_i
w.x.array[:] += deltax.x.array
i += 1

# Compute norm of update
correction_norm = deltax.vector.norm(0)

# Compute L2 error comparing to the analytical solution
#L2_error.append(np.sqrt(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(error), op=MPI.SUM)))
dx_norm.append(correction_norm)

print(f"Iteration {i}: Correction norm {correction_norm}")
if correction_norm < 1e-10:
break



and the mesh is here: 446.5 KB folder on MEGA

It finishes at the first iteration and return 0.0

I dont find the problem. does anyone know why?

Could you simplify your example to trying to solve a problem, with a built in mesh, say: Test problem 1: Channel flow (Poiseuille flow) — FEniCSx tutorial
as this would make it easier to eliminate programming errors, and easier for people to reproduce

Here is:

import dolfinx
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import math
import ufl
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement,
as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym, MixedElement)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
create_vector, create_matrix, set_bc)

from dolfinx.fem import (Constant, Function, FunctionSpace,
assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc)
from dolfinx.mesh import create_unit_square

##########################################################################################################3

domain = create_unit_square(MPI.COMM_WORLD, 10, 10)

v_cg1 = VectorElement("CG", domain.ufl_cell(), 2)
s_cg1 = FiniteElement("CG", domain.ufl_cell(), 1)
V = dolfinx.fem.FunctionSpace(domain, v_cg1)
Q = dolfinx.fem.FunctionSpace(domain, s_cg1)
Z = MixedElement([v_cg1, s_cg1])
W = dolfinx.fem.FunctionSpace(domain, Z)

# Boundary subdomains and coditions

########################################################################################################

# Definition of subdomains

def upstream(x):
return np.isclose(x[0],0)

def downstream(x):
return np.isclose(x[0],1)

def walls(x):
return np.logical_or(np.isclose(x[1],0),np.isclose(x[1],1))

fdim = domain.topology.dim -1

up_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, upstream)
down_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, downstream)
wall_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, walls)

def u_exact(x):
values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = 1
return values

u_icondition = Function(V)
p_icondition = Function(Q)
u_icondition.interpolate(u_exact)

up_f = dolfinx.fem.locate_dofs_topological(W.sub(0), fdim, up_facets)
wall_f = dolfinx.fem.locate_dofs_topological(W.sub(0), fdim, wall_facets)
down_f = dolfinx.fem.locate_dofs_topological(W.sub(1), fdim, down_facets)

#######################################################################################################################

# Bounday conditions

bc_u = dirichletbc(u_icondition, up_f)
bc_w = dirichletbc(u_icondition, wall_f)
bc_p = dirichletbc(p_icondition, down_f)

bc = [bc_u, bc_w, bc_p]
mu = Constant(domain, PETSc.ScalarType(0.02))
n = FacetNormal(domain)
def epsilon(u):

# Define stress tensor
def sigma(u, p):
return 2*mu*epsilon(u) - p*Identity(len(u))

w = Function(W)
u, p = ufl.split(w)

m = TestFunction(W)

v, q = ufl.split(m)

F = dot(dot(u,nabla_grad(u)),v)*dx + div(u)*q*dx + inner(sigma(u,p),epsilon(v))*dx + dot(p*n, v)*ds - dot(mu*nabla_grad(u)*n, v)*ds

J = ufl.derivative(F, w)
residual = dolfinx.fem.form(F)
jacobian = dolfinx.fem.form(J)

deltax = dolfinx.fem.Function(W)
A = dolfinx.fem.petsc.create_matrix(jacobian)
L = dolfinx.fem.petsc.create_vector(residual)
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
i = 0
#L2_error = []
dx_norm = []
max_iterations = 25
while i < max_iterations:
# Assemble Jacobian and residual
with L.localForm() as loc_L:
loc_L.set(0)
A.zeroEntries()
dolfinx.fem.petsc.assemble_matrix(A, jacobian, bcs=bc)
A.assemble()

dolfinx.fem.petsc.assemble_vector(L, residual)
L.scale(-1)

# Compute b - J(u_D-u_(i-1))
dolfinx.fem.petsc.apply_lifting(L, [jacobian], [bc], x0=[w.vector], scale=1)
# Set dx|_bc = u_{i-1}-u_D
dolfinx.fem.petsc.set_bc(L, bc, w.vector, 1.0)

# Solve linear problem
solver.solve(L, deltax.vector)
deltax.x.scatter_forward()

# Update u_{i+1} = u_i + delta x_i
w.x.array[:] += deltax.x.array
i += 1

# Compute norm of update
correction_norm = deltax.vector.norm(0)

# Compute L2 error comparing to the analytical solution
#L2_error.append(np.sqrt(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(error), op=MPI.SUM)))
dx_norm.append(correction_norm)

print(f"Iteration {i}: Correction norm {correction_norm}")
if correction_norm < 1e-10:
break


in this case I do not know if the same problem is present because now the error is

carlos@carlos-PROX15-AMD:~/Descargas\$ python3 scrpit.py
malloc(): unsorted double linked list corrupted
[carlos-PROX15-AMD:09475] *** Process received signal ***
[carlos-PROX15-AMD:09475] Signal: Aborted (6)
[carlos-PROX15-AMD:09475] Signal code:  (-6)



You are not setting your boundary condition correctly, as you are using a mixed finite element. Ref: Incompressible Navier-Stokes with particular scheme & mixed elements FEniCSx - #6 by dokken

That solves the kernel issue, but now I have a problem with defining the boundary condition.

bc_p = dirichletbc(ScalarType(0), down_f,W.sub(1))
bc_u = dirichletbc(ScalarType((1,0)), up_f, W.sub(0))
bc_w = dirichletbc(ScalarType((1,0)), wall_f, W.sub(0))


returns

 Creating a DirichletBC using a Constant is not supported when the Constant size is not equal to the block size of the constrained (sub-)space. Use a fem::Function to create the fem::DirichletBC.


I don’t have a problem defining the condition to a Scalar function.

I tried

bc_p = dirichletbc(ScalarType(0), down_f,W.sub(1))
bc_u1 = dirichletbc(ScalarType(1), up_f, W.sub(0).sub(0))
bc_u1 = dirichletbc(ScalarType(0), up_f, W.sub(0).sub(1))
bc_w0 = dirichletbc(ScalarType(1), wall_f,W.sub(0).sub(0))
bc_w1 = dirichletbc(ScalarType(0), wall_f,W.sub(0).sub(1))

bc = [bc_u0,bc_u1,bc_w0, bc_w1, bc_p]


but return w = 0…

As stated you should then use a function to define the bc values, as for instance done in

you can interpolate any expression into the function.

Thank you so much @dokken !!!