# Issues with Navier-Stokes problem with monolithic solver in mixed elements space (dolfinx)

Dear all,

I’m having an issue solving Navier-Stokes problem in a inf-sup stable mixed elements space with a monolithic solver in dolfinx. More specifically, I’m trying to refactor the following code (applied to a toy problem) that works correctly in legacy fenics:

from dolfin import *
import numpy as np

msh = UnitSquareMesh(20,20)

# set fluid parameters
mu = 0.035
rho = 1.0

# spaces
VE = VectorElement('P', msh.ufl_cell(), 2)
PE = FiniteElement('P', msh.ufl_cell(), 1)

elem = MixedElement([VE, PE])
W = FunctionSpace(msh, elem)

# define gamma function
V = FunctionSpace(msh, 'P', 1)
gamma_function = Function(V)
vec = gamma_function.vector()
values = vec.get_local()

dofmap = V.dofmap()
my_first, my_last = dofmap.ownership_range() # global

visited = []

# 'Handle' API change of tabulate coordinates
X = V.tabulate_dof_coordinates()
X = X.reshape((-1, msh.geometry().dim()))

for cell in cells(msh):
dofs = dofmap.cell_dofs(cell.index()) # local
for dof in dofs:
if not dof in visited:
global_dof = dofmap.local_to_global_index(dof)  # global
if my_first <= global_dof < my_last:
if np.isclose(X[dof][1], 0.4, atol=1e-2) and X[dof][0] >= 0.3:
visited.append(dof)
values[dof] = 1.

# apply values to the vectorize function
vec.set_local(values)
vec.apply('insert')

R = 1e8
gamma_function.vector()[:] = R*gamma_function.vector()[:]

# time settings
T = 0.4
t_end = 0.2
ts = 20
Dt = t_end/ts
t = 0.

# boundary conditions
def inflow(x):
return np.isclose(x[1], 0.)

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

u_in = Expression('(16*x[0]-16*pow(x[0],2))*sin(M_PI/T*t)', t=t, T=T, degree=2)
bc_y = DirichletBC(W.sub(0).sub(1), u_in, inflow)
bc_x = DirichletBC(W.sub(0).sub(0), Constant(0.), inflow)
bc_walls = DirichletBC(W.sub(0), Constant((0., 0.)), walls)
bcs = [bc_x, bc_y, bc_walls]

# mark outflow to apply backflow stabilization
sub_domains = MeshFunction("size_t", msh, msh.topology().dim() - 1)
sub_domains.set_all(0)

class Outflow(SubDomain):
def inside(self, x, on_boundary):
return np.isclose(x[0], 1.) and on_boundary

outflow = Outflow()
outflow.mark(sub_domains, 3)

ds = Measure('ds', domain=msh, subdomain_data=sub_domains)

# variational form
w = Function(W, name='w')
u0, p0 = w.split()

u, p = TrialFunctions(W)
v, q = TestFunctions(W)
nu = FacetNormal(msh)

a = rho*inner(u,v)*dx \
- inner(p, div(v))*dx \
+ inner(div(u), q)*dx) \
+ Dt*inner(gamma_function*u,v)*dx
L = rho*inner(u0,v)*dx

a += 0.25*rho*Dt*( abs( inner( u0,nu ) - abs(inner( u0,nu )) )*inner( u,v )*ds(3) ) #backflow stabilization
a += 0.5*rho*Dt*( div(u0)*inner( u,v )*dx ) # Temam stabilization

while (t <= t_end):
t += Dt
u_in.t = t

solve(a == L, w, bcs)

u0, p0 = w.split(deepcopy=True)

xdmf = XDMFFile('./u_2D_fenics.xdmf')
xdmf.write(u0)


Such code returns this velocity field at t_end, with the expected convection induced behaviour (jet formation and recirculation):

I also have a dolfinx version of this code:

import numpy as np

from dolfinx.mesh import create_unit_square, locate_entities, meshtags
from dolfinx.fem import (Constant, Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,
locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block, LinearProblem
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
from ufl import (FacetNormal, FiniteElement, MixedElement, Identity, TestFunction, TrialFunction, VectorElement,
TrialFunctions, TestFunctions, Measure, split, div, grad, dot, ds, dx, inner, lhs, nabla_grad, rhs, sym)

msh = create_unit_square(MPI.COMM_WORLD, 20, 20)

# set fluid parameters
mu = 0.035
rho = 1.0

# spaces
VE = VectorElement('P', msh.ufl_cell(), 2)
PE = FiniteElement('P', msh.ufl_cell(), 1)

elem = MixedElement([VE, PE])
W = FunctionSpace(msh, elem)

# define gamma function
V = FunctionSpace(msh, ('P', 1))
gamma_function = Function(V, name='f', dtype=np.float64)

dofmap = V.dofmap
X = V.tabulate_dof_coordinates()
visited = []

Map = msh.topology.index_map(msh.topology.dim)
num_cells = Map.size_local + Map.num_ghosts # total number of cells

indices = []
values = []

for cell_idx in range(num_cells):
for dof in vertex_glob_idx:
if not dof in visited:
if np.isclose(X[dof][1], 0.4, atol=1e-2) and X[dof][0] >= 0.3:
visited.append(dof)
indices.append(dof)
values.append(1.)

indices = np.array(indices, dtype=np.int32)
values = np.array(values)
gamma_function.x.array[indices] = values
gamma_function.x.scatter_forward()

R = 1e8
gamma_function.x.array[:] = R*gamma_function.x.array[:]

# time settings
T = 0.4
t_end = 0.2
ts = 20
Dt = t_end/ts
t = 0.

# boundary conditions
def inflow(x):
return np.isclose(x[1], 0.)

class time_dependent_inflow:
def __init__(self):
self.t = 0.0

def eval(self, x):
return np.stack((np.zeros(x.shape[1]), np.full(x.shape[1], (16*x[0]-16*x[0]**2)*np.sin(np.pi/T*self.t)))) #-U/(rr**2)*(x[0]**2-rr**2)*np.sin(np.pi/T*self.t)

V0, _ = W.sub(0).collapse()
inflow_expr = time_dependent_inflow()
inflow_expr.t = 0.
bc_fun = Function(V0)
bc_fun.interpolate(inflow_expr.eval)
inflow_dofs = locate_dofs_geometrical((W.sub(0),V0), inflow)
bcu = dirichletbc(bc_fun, inflow_dofs, W.sub(0))

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

u_bc = Function(V0)
facets = df.mesh.locate_entities_boundary(msh, 1, walls)
bc_noslip = dirichletbc(u_bc, df.fem.locate_dofs_topological((W.sub(0), V0), 1, facets), W.sub(0))

bcs = [bcu, bc_noslip]

bnd = [(1, lambda x: np.isclose(x[0], 1.)),
(2, lambda x: np.isclose(x[0], 0.)),
(3, lambda x: np.isclose(x[1], 1.)),
(4, lambda x: np.isclose(x[1], 0.))]
facet_indices, facet_markers = [], []
fdim = msh.topology.dim - 1
for (marker, locator) in bnd:
facets = locate_entities(msh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(msh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = Measure("ds", domain=msh, subdomain_data=facet_tag) # this is the important part!
# facet_tag = tagged surface mesh

# variational form
w = Function(W)
u_n, p_n = split(w)

u, p = TrialFunctions(W)
v, q = TestFunctions(W)
n = FacetNormal(msh)

a = rho*inner(u,v)*dx \
- inner(p, div(v))*dx \
+ inner(div(u), q)*dx) \
+ Dt*inner(gamma_function*u,v)*dx
L = rho*inner(u_n,v)*dx

a += 0.25*rho*Dt*( abs( inner( u_n,n ) - abs(inner( u_n,n )) )*inner( u,v )*ds(3) ) # backflow stabilization
a += 0.5*rho*Dt*( div(u_n)*inner( u,v )*dx ) # temam stabilization

# time loop
while t < t_end:

t += Dt
inflow_expr.t = t
bc_fun.interpolate(inflow_expr.eval)

problem = LinearProblem(a, L, bcs, petsc_options={"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "superlu_dist"})
w = problem.solve()

u_n, p_n = w.split() # deepcopy=True not available in dolfinx

with df.io.XDMFFile(msh.comm, './u_new.xdmf', 'w') as xdmf:
u_n.x.scatter_forward()
xdmf.write_mesh(msh)
xdmf.write_function(u_n)


However, this code doesn’t seem to be working correctly: the velocity field at t_end looks pretty different, as if there was no convection involved. Actually, looking only at the velocity field I’d say it comes from a Stokes problem. Here you can take a look at it:

My best guess is that there’s something wrong in the way I use the functions coming from the mixed elements space, but I couldn’t find anything helpful on the forum.

Giorgia

Here you are making many assumptions with respect to the dofmap of the finite element space V and the vertices.

Are you simply trying to find those dofs in V where y is in the range (0.38-0.42) and x>=0.3?

Secondly. What happens if you remove the object? how does the flow field look like?

Thirdly, why is the linear problem defined inside the loop? this just adds overhead by recreating existing structures.

Finally, how does the solution look like after a single time step?

1 Like

First, what I’m trying to do is creating a function defined on the nodes of the mesh that will be used in the variational formulation to penalize the momentum equation. Such penalization will act in some sense as an internal homogeneous Dirichlet boundary condition, allowing us to work with very simple meshes and defining the obstacle(s) through the problem formulation. This is just a toy problem so yes, in this case I’m looking for those nodes. However, the actual problem we are solving is 3D and the obstacle has a much more complex shape than this one, but the way it is computed is exactly the same (using if conditions to determine which nodes are non-zero), as well as the weird results.

Second, removing the object produces the following results at t = t_end (top legacy fenics, bottom fenicsx):

Even in this case I still don’t see the expected jet flow, even if it reaches the same velocity values at the inlet, which makes me assume that the inflow boundary condition is working correctly.

Thirdly, thanks for pointing this out, I now took it out of the loop but the solution did not change.

Finally, here’s the solution after the first time step, t = 0.01 (top legacy fenics, bottom fenicsx):

The streamlines seem to be more similar, even if we have some differences in the velocity magnitude.

Another hypothesis that comes to my mind is that u_n is not correctly updated at each time iteration, hence the convective behaviour doesn’t show up throughout the loop, but I have no idea how to fix it.
Moreover, u_n comes from the w = ufl.split() that should permit direct access to the components of the mixed function, that’s why it should work.

The line

    w = problem.solve()


defines a new function that is not at all “connected” to the function w defined at the beginning.
The convecting velocity/previous time step function u_n is never updated.

This should work:

V0, subspace_map = W.sub(0).collapse()
u_n = fem.Function(V0)   # convecting/last time step velocity

# def. linear forms etc
# ...

while t < t_end:
# ...
w = problem.solve()
u_n.x.array[:] = w.x.array[subspace_map]

1 Like

Thank you very much for the explanation, this solved the problem! I was genuinely convinced the two w referred to the same object, that wasn’t the case.