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

[Karlosjusus]

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!

[dokken]

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

[Karlosjusus]

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)

[dokken]

@Karlosjusus

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)

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

[Karlosjusus]

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…

[dokken]

@Karlosjusus

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.

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

Boundary conditions about Mixed space

Please make the example reproducible. Currently you have only posted snippets of code, which doesn’t make it easy for people to reproduce the problem. An example of how to make your code reproducible is: import dolfinx from mpi4py import MPI import ufl from petsc4py import PETSc import numpy as np domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10) Ve = ufl.VectorElement("CG", domain.ufl_cell(), 2) Pe = ufl.FiniteElement("CG", domain.ufl_cell(), 1) We = ufl.MixedElement([Ve, Pe]) …

How to impose a Dirichlet bc on a vector sub-space?

As the error message says, the first input has to be a function. In your code you invoke it with a Constant value: V1, _ = ME.sub(1).collapse() u_D = fem.Function(V1) u_D.x.array[:] =1 boundary_dofs_b = fem.locate_dofs_topological((ME.sub(1), V1), fdim, boundary_facets) bc_b = fem.dirichletbc(u_D, boundary_dofs_b, ME.sub(1))

you can interpolate any expression into the function.

[Karlosjusus]

Thank you so much @dokken !!!