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

General
1

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

    # Read the mesh

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

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

    # 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):
    return sym(nabla_grad(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.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
    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)
    L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)

    # 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?
Thank you in advance!

2

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

3

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

    # Read the mesh

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):
    return sym(nabla_grad(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
#error = dolfinx.fem.form(ufl.inner(uh - u_ufl, uh - u_ufl) * ufl.dx(metadata={"quadrature_degree": 4}))
#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.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
    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)
    L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)

    # 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)
4

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

5

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…

6

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.

7

Thank you so much @dokken !!!