BC application when using non linear solver (dolfinx_mpc)

[gianmon99]

I’m trying to implement a RANS simulation on a periodic hill (using the k-epsilon turbulence model)
However it seems like I’m unable to apply the BC to the speed ad the pressure when I’m using the NonLinear solver of dolfinx_mpc.
(I’m using the following versions of dolfinx and dolfinx_mpc:
dolfinx_mpc —> Version: 0.10.0
fenics-dolfinx —> Version: 0.10.0)
(Also for the non_linear solver I took inspiration from to code: dolfinx_mpc/python/demos/demo_stokes_nonlinear_nest.py at main · jorgensd/dolfinx_mpc · GitHub)
Can someone spot what am I doing wrong?

This is a short version of the important part of the code:

# Construct function spaces
cellname = mesh.basix_cell()
Ve = basix.ufl.element(basix.ElementFamily.P, cellname, 2, shape=(gdim,), dtype=default_real_type) # Lagrange 2 function space
Qe = basix.ufl.element(basix.ElementFamily.P, cellname, 1, dtype=default_real_type) # Legrange 1 function space
el_k = basix.ufl.element("Lagrange", mesh.basix_cell(), 1) 
V = functionspace(mesh, Ve)
Q = functionspace(mesh, Qe)
W = ufl.MixedFunctionSpace(V, Q)

# Costruct the boundary conditions: 
bcu = []
bcp = []
bcw = [] # bcu+bcp

# Reminder of the facet markers:
# 1 ---> inlet 
# 2 ---> outlet 
# 3 ---> top 
# 4 ---> bottom 

def periodic_relation(x):
        out_x = np.zeros(x.shape)
        out_x[0] = x[0] - Lx 
        out_x[1] = x[1]
        return out_x

#----------------------   DIREHLET BC   -----------------------#
# Velocity (ux, uy)
# Top or bottom walls
u_no_slip= Function(V)
u_no_slip.x.array[:]=0
u_no_slip.x.scatter_forward() # <-- updates ghost values
# Pressure:
p_inlet = Constant(mesh, dolfinx.default_scalar_type(2.0))
p_outlet = Constant(mesh, dolfinx.default_scalar_type(0.0))

# k-epsilon:
# Top and bottom walls:
k_walls = Constant(mesh, dolfinx.default_scalar_type(0.0))

# VELOCITY (ux,uy):
# Top and bottom walls 
dofs_top_u = locate_dofs_topological(V, fdim, ft.find(3))
dofs_bottom_u = locate_dofs_topological(V, fdim, ft.find(4))
bc_top_u = dirichletbc(u_no_slip, dofs_top_u)
bc_bottom_u = dirichletbc(u_no_slip, dofs_bottom_u)
bcu.append(bc_top_u)
bcu.append(bc_bottom_u)
# Pressure(p)
dofs_in_p = locate_dofs_topological(Q, fdim, ft.find(1))
dofs_out_p = locate_dofs_topological(Q, fdim, ft.find(2))
bc_in_p = dirichletbc(p_inlet, dofs_in_p, Q) 
bc_out_p = dirichletbc(p_outlet, dofs_out_p, Q) 
bcp.append(bc_in_p)
bcp.append(bc_out_p)

bcw = np.concatenate((bcu, bcp))

#-------

# Defining the mpc:
mpc_u = MultiPointConstraint(V)
mpc_u.create_periodic_constraint_topological(V, ft, 2, periodic_relation, bcu) # or bcw?
mpc_u.finalize()

mpc_p = MultiPointConstraint(Q) # or, should I apply bcp????
mpc_p.finalize()


# Initialise the functions 
Wh = ufl.MixedFunctionSpace(mpc_u.function_space, mpc_p.function_space)
u0 = dolfinx.fem.Function(mpc_u.function_space)
p0 = dolfinx.fem.Function(mpc_p.function_space)

u1 = dolfinx.fem.Function(mpc_u.function_space)
p1 = dolfinx.fem.Function(mpc_p.function_space)
(v, q) = ufl.TestFunctions(W)

# SOLVING THE NS:
#  (these are actually solved in a loop, since they are "linearised")
FW  = dot(dot(u0, nabla_grad(u1)), v)*dx\
    + (nu + nu_t) * inner(nabla_grad(u1), nabla_grad(v))*dx \
    - div(v)*p1*dx \
    - div(u1)*q*dx \
    - dot(force, v)*dx \
# Surface terms of the NS equation:
FW+= dot(p1*normal, v)*ds - dot((nu + nu_t)*nabla_grad(u1)*normal, v)*ds  
  
tol = 1e-7
problem_ns = dolfinx_mpc.NonlinearProblem(
ufl.extract_blocks(FW),
[u1, p1],
mpc=[mpc_u, mpc_p],
bcs=bcw,
kind="nest",
petsc_options={
    "snes_type": "newtonls",
    "snes_atol": 1e-7,
    "snes_rtol": 1e-7,
    #"snes_monitor": None,

    "ksp_type": "preonly", # the linear solver just solves the linear system (no iterations)
    "pc_type": "lu", # LU factorization for the pre-conditioning
}
    )

start_ns = time.time()
_, converged, num_iterations = problem_ns.solve()

These are the velocity and pressure field!
(first image the velocity magnitude, second image the pressure field)
They are both zero! Since no deltaP is being applied:

[dokken]

First of all, the indentation of this code is all over the place, making it hard to say if our problem works or not.
Secondly, try with the standard non-linear problem in DOLFINx (without MPC) first, to ensure that your setup is correct.
Thirdly, start with a simple channel (built in box), and make the example reproducible for others.
Finally, if you are linearizing the equation

@gianmon99

FW  = dot(dot(u0, nabla_grad(u1)), v)*dx\
    + (nu + nu_t) * inner(nabla_grad(u1), nabla_grad(v))*dx \
    - div(v)*p1*dx \
    - div(u1)*q*dx \
    - dot(force, v)*

you don’t need to use the Nonlinearproblem, and can solve this with linearproblem. if you rather use dot(dot(u1, nabla_grad(u1))) you would need to use a non-linear solver.

[gianmon99]

Sorry for the strange indentation in the question (in the code the indentation is correct)
However, as you suggested I’ve implemented a simpler version of the code on a simple channel (also I’m now using a linear solver).
But, I’m having the same problem, no bc is applied on the pressure at inlet and outlet (is this because I have imposed periodic Bc for the velocity)?

This is the channel code; so everyone can replicate:

import ufl
import dolfinx
from mpi4py import MPI
from dolfinx import mesh, default_real_type, default_scalar_type
from dolfinx.fem import functionspace, Function, Constant, dirichletbc, locate_dofs_topological, assemble_scalar, form 
from dolfinx import mesh
from pathlib import Path
import numpy as np
import basix.ufl
from ufl import (dot, inner, grad, div, dx, ds, nabla_grad, lhs, rhs) 
import dolfinx_mpc
from dolfinx_mpc import MultiPointConstraint




mesh_comm = MPI.COMM_WORLD
model_rank = 0


# create mesh 
height = 2.0; width = 2.0; nx = 40; ny = 600

def meshing(lx,ly, Nx,Ny):
    m = mesh.create_unit_square(mesh_comm, Nx, Ny)
    x = m.geometry.x

    #Refine on top and bottom sides
    x[:, 1] = 0.5 * ( np.arctan(2 * np.pi * (x[:,1] - 0.5)) / np.arctan(np.pi) + 1)

    #Scale
    x[:,0] = x[:,0]*lx
    x[:,1] = x[:,1]*ly
    return m

mesh = meshing(width, height, nx, ny)

nodes_coords = mesh.geometry.x
print(f'nodes of the mesh: {nodes_coords.shape[0]} \n')

gdim = mesh.geometry.dim
fdim = gdim -1 

# Defining the domain markers (that will be used later for the BC):
facet_indices, facet_markers = [], []
mesh.topology.create_connectivity(fdim, gdim)

def walls_marker(x):
    return np.isclose(x[1], 0.0) | np.isclose(x[1], height)
walls_facets= dolfinx.mesh.locate_entities_boundary(mesh, fdim, walls_marker)
facet_indices.append(walls_facets)
facet_markers.append(np.full_like(walls_facets,3)) # 3-->walls marker 

def inl_marker(x):
    return np.isclose(x[0],0.0)
inlet_facets= dolfinx.mesh.locate_entities_boundary(mesh, fdim, inl_marker)
facet_indices.append(inlet_facets)
facet_markers.append(np.full_like(inlet_facets,1)) # 1-->inlet marker 

def out_marker(x):
    return np.isclose(x[0],width)
out_facets= dolfinx.mesh.locate_entities_boundary(mesh, fdim, out_marker)
facet_indices.append(out_facets)
facet_markers.append(np.full_like(out_facets,2)) # 2-->outlet marker 

# Creating the meshargs object for easier BC imposition
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)

sorted_facets = np.argsort(facet_indices) 
ft = dolfinx.mesh.meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets]) # Create a MeshTags object that associates data with a subset of mesh entities.
# ft.find(1) <-- inlet facets 
# ft.find(2) <-- outlet facets
# ft.find(3) <-- walls facets 


normal = ufl.FacetNormal(mesh) 
# Construct function spaces
P2 = basix.ufl.element("Lagrange", mesh.topology.cell_name(), 2, shape=(mesh.geometry.dim,), dtype=default_real_type)
P1 = basix.ufl.element("Lagrange", mesh.topology.cell_name(), 1, dtype=default_real_type)

# 2. Crea l'elemento misto
Mel = basix.ufl.mixed_element([P2, P1])

# 3. Crea il FunctionSpace misto DOLFINx
W = functionspace(mesh, Mel)  # Ora W è un dolfinx.fem.FunctionSpace

# Nota: Non ti servono V e Q separati qui sopra per creare W.
# Se ti servono per definire le funzioni (es. u_no_slip), li ottieni collassando i sottospazi di W:
V, V_to_W_map = W.sub(0).collapse()
Q, Q_to_W_map = W.sub(1).collapse()

# V = dolfinx.fem.functionspace(mesh, P2)
# Q = dolfinx.fem.functionspace(mesh, P1)
# W = MixedFunctionSpace(V, Q)

# Costruct the boundary conditions: 
bcu = []
bcp = []
bcw = [] # bcu+bcp

def periodic_relation(x):
        out_x = np.zeros(x.shape)
        out_x[0] = x[0] - width 
        out_x[1] = x[1]
        return out_x

#----------------------   DIREHLET BC   -----------------------#
# Velocity (ux, uy)
# Top or bottom walls
u_no_slip= Function(V)
u_no_slip.x.array[:]=0
u_no_slip.x.scatter_forward() # <-- updates ghost values
# Pressure:
p_in = Constant(mesh, dolfinx.default_scalar_type(2.0))
p_out = Constant(mesh, dolfinx.default_scalar_type(0.0))

#----------------------   DIRECT BC   -----------------------#
# --- BC Velocità ---
u_no_slip = Function(V) # Definito su V (collassato)
u_no_slip.x.array[:] = 0.0
dofs_walls_u = locate_dofs_topological((W.sub(0), V), fdim, ft.find(3)) # Nota la tupla (W.sub, V)
bc_walls_u = dirichletbc(u_no_slip, dofs_walls_u, W.sub(0)) # Applicato su W.sub(0)
bcu.append(bc_walls_u)
# Pressure(p)
dofs_in_p = locate_dofs_topological(W.sub(1), fdim, ft.find(1))
dofs_out_p = locate_dofs_topological(W.sub(1), fdim, ft.find(2))
bc_in_p = dirichletbc(p_in, dofs_in_p, W.sub(1))
bc_out_p = dirichletbc(p_out, dofs_out_p, W.sub(1)) 
bcp.append(bc_in_p)
bcp.append(bc_out_p) 

bcw = np.concatenate((bcu, bcp))


mpc = MultiPointConstraint(W)
mpc.create_periodic_constraint_topological(W.sub(0), ft, 2, periodic_relation, bcu) 
mpc.finalize()


# Initialize constants and expressions
Re = 10
nu = Constant(mesh, default_scalar_type(1/Re))
# --- Set forcing ---
alpha = 0  # iniziale ( U_bulk=..... con alpha= -1)
force = Constant(mesh, default_scalar_type((alpha, 0)))

# Initialise the functions 
(u,p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)
u0 = Function(V)
u1 = Function(V)
p0 = Function(Q)
p1 = Function(Q)

FW  = dot(dot(u0, nabla_grad(u)), v)*dx\
    + (nu) * inner(nabla_grad(u), nabla_grad(v))*dx \
    - div(v)*p*dx \
    - div(u)*q*dx \
    - dot(force, v)*dx
# Surface terms non considered:
#FW+= dot(p1*normal, v)*ds - dot((nu)*nabla_grad(u1)*normal, v)*ds 

a = form(lhs(FW)); L = form(rhs(FW))

problem = dolfinx_mpc.LinearProblem(a, L, mpc, bcs=bcw)
wh = problem.solve()
u1 = wh.sub(0).collapse()
p1 = wh.sub(1).collapse()

This is a screen of the pressure an velocity magnitude obtained (everything is zero!):

[dokken]

You have not set any solver options:

@gianmon99

problem = dolfinx_mpc.LinearProblem(a, L, mpc, bcs=bcw)

replacing this with

problem = dolfinx_mpc.LinearProblem(
    a,
    L,
    mpc,
    bcs=bcw,
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
        "ksp_monitor": None,
        "ksp_error_if_not_converged": True,
    },
)

yields a solution, but not what you would expect.

This is due to a bug in:

@gianmon99

mpc.create_periodic_constraint_topological(W.sub(0), ft, 2, periodic_relation, bcu)

which you can replace with:

pc = MultiPointConstraint(W)
mpc.create_periodic_constraint_topological(
    W.sub(0).sub(0), ft, 2, periodic_relation, bcu
)
mpc.create_periodic_constraint_topological(
    W.sub(0).sub(1), ft, 2, periodic_relation, bcu
)
mpc.finalize()

i.e. constrain each velocity component individually, which as a complete code yields:

import ufl
import dolfinx
from mpi4py import MPI
from dolfinx import default_real_type, default_scalar_type
from dolfinx.fem import (
    functionspace,
    Function,
    Constant,
    dirichletbc,
    locate_dofs_topological,
    assemble_scalar,
    form,
)
import numpy as np
import basix.ufl
from ufl import dot, inner, grad, div, dx, ds, nabla_grad, lhs, rhs
import dolfinx_mpc
from dolfinx_mpc import MultiPointConstraint


mesh_comm = MPI.COMM_WORLD
model_rank = 0


# create mesh
height = 1.0
width = 1.0
nx = 50
ny = 100


def meshing(lx, ly, Nx, Ny):
    m = dolfinx.mesh.create_unit_square(mesh_comm, Nx, Ny)
    x = m.geometry.x

    # Refine on top and bottom sides
    x[:, 1] = 0.5 * (np.arctan(2 * np.pi * (x[:, 1] - 0.5)) / np.arctan(np.pi) + 1)

    # Scale
    x[:, 0] = x[:, 0] * lx
    x[:, 1] = x[:, 1] * ly
    return m


mesh = meshing(width, height, nx, ny)

tdim = mesh.geometry.dim
fdim = tdim - 1

# Defining the domain markers (that will be used later for the BC):
facet_indices, facet_markers = [], []
mesh.topology.create_connectivity(fdim, tdim)


def walls_marker(x):
    return np.isclose(x[1], 0.0) | np.isclose(x[1], height)


walls_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, walls_marker)
facet_indices.append(walls_facets)
facet_markers.append(np.full_like(walls_facets, 3))  # 3-->walls marker


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


inlet_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, inl_marker)
facet_indices.append(inlet_facets)
facet_markers.append(np.full_like(inlet_facets, 1))  # 1-->inlet marker


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


out_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, out_marker)
facet_indices.append(out_facets)
facet_markers.append(np.full_like(out_facets, 2))  # 2-->outlet marker

# Creating the meshargs object for easier BC imposition
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)

sorted_facets = np.argsort(facet_indices)
ft = dolfinx.mesh.meshtags(
    mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets]
)  # Create a MeshTags object that associates data with a subset of mesh entities.
# ft.find(1) <-- inlet facets
# ft.find(2) <-- outlet facets
# ft.find(3) <-- walls facets

normal = ufl.FacetNormal(mesh)
# Construct function spaces
P2 = basix.ufl.element(
    "Lagrange",
    mesh.topology.cell_name(),
    2,
    shape=(mesh.geometry.dim,),
    dtype=default_real_type,
)
P1 = basix.ufl.element(
    "Lagrange", mesh.topology.cell_name(), 1, dtype=default_real_type
)

# 2. Crea l'elemento misto
Mel = basix.ufl.mixed_element([P2, P1])

# 3. Crea il FunctionSpace misto DOLFINx
W = functionspace(mesh, Mel)  # Ora W è un dolfinx.fem.FunctionSpace

# Nota: Non ti servono V e Q separati qui sopra per creare W.
# Se ti servono per definire le funzioni (es. u_no_slip), li ottieni collassando i sottospazi di W:
V, V_to_W_map = W.sub(0).collapse()
Q, Q_to_W_map = W.sub(1).collapse()

# V = dolfinx.fem.functionspace(mesh, P2)
# Q = dolfinx.fem.functionspace(mesh, P1)
# W = MixedFunctionSpace(V, Q)

# Costruct the boundary conditions:
bcu = []
bcp = []
bcw = []  # bcu+bcp


def periodic_relation(x):
    out_x = np.zeros(x.shape)
    out_x[0] = x[0] - width
    out_x[1] = x[1]
    return out_x


# ----------------------   DIREHLET BC   -----------------------#
# Velocity (ux, uy)
# Top or bottom walls
u_no_slip = Function(V)
u_no_slip.x.array[:] = 0
u_no_slip.x.scatter_forward()  # <-- updates ghost values
# Pressure:
p_in = Constant(mesh, dolfinx.default_scalar_type(2.0))
p_out = Constant(mesh, dolfinx.default_scalar_type(0.0))

# ----------------------   DIRECT BC   -----------------------#
# --- BC Velocità ---
u_no_slip = Function(V)  # Definito su V (collassato)
u_no_slip.x.array[:] = 0.0
dofs_walls_u = locate_dofs_topological(
    (W.sub(0), V), fdim, ft.find(3)
)  # Nota la tupla (W.sub, V)
bc_walls_u = dirichletbc(u_no_slip, dofs_walls_u, W.sub(0))  # Applicato su W.sub(0)
bcu.append(bc_walls_u)
# Pressure(p)
dofs_in_p = locate_dofs_topological(W.sub(1), fdim, ft.find(1))
dofs_out_p = locate_dofs_topological(W.sub(1), fdim, ft.find(2))
bc_in_p = dirichletbc(p_in, dofs_in_p, W.sub(1))
bc_out_p = dirichletbc(p_out, dofs_out_p, W.sub(1))
bcp.append(bc_in_p)
bcp.append(bc_out_p)

bcw = np.concatenate((bcu, bcp))

with dolfinx.io.XDMFFile(mesh_comm, "mesh_mwe_10.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(ft, mesh.geometry)

mpc = MultiPointConstraint(W)
mpc.create_periodic_constraint_topological(
    W.sub(0).sub(0), ft, 2, periodic_relation, bcu
)
mpc.create_periodic_constraint_topological(
    W.sub(0).sub(1), ft, 2, periodic_relation, bcu
)
mpc.finalize()


# Initialize constants and expressions
Re = 10
nu = Constant(mesh, default_scalar_type(1 / Re))
# --- Set forcing ---
alpha = 0  # iniziale ( U_bulk=..... con alpha= -1)
force = Constant(mesh, default_scalar_type((alpha, 0)))

# Initialise the functions
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)
u0 = Function(V)
u1 = Function(V)
p0 = Function(Q)
p1 = Function(Q)

FW = (
    dot(dot(u0, nabla_grad(u)), v) * dx
    + (nu) * inner(nabla_grad(u), nabla_grad(v)) * dx
    - div(v) * p * dx
    - div(u) * q * dx
    - dot(force, v) * dx
)
# Surface terms non considered:
# FW+= dot(p1*normal, v)*ds - dot((nu)*nabla_grad(u1)*normal, v)*ds

a = form(lhs(FW))
L = form(rhs(FW))

problem = dolfinx_mpc.LinearProblem(
    a,
    L,
    mpc,
    bcs=bcw,
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
        "ksp_monitor": None,
        "ksp_error_if_not_converged": True,
    },
)
wh = problem.solve()
u1 = wh.sub(0).collapse()
p1 = wh.sub(1).collapse()
with dolfinx.io.VTXWriter(mesh.comm, "u.bp", [u1]) as vtx:
    vtx.write(0.0)
with dolfinx.io.VTXWriter(mesh.comm, "p.bp", [p1]) as vtx:
    vtx.write(0.0)

which yields