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 \
+ Dt*rho*inner(grad(u)*u0,v)*dx +Dt*(mu * inner(grad(u), grad(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):
vertex_glob_idx = msh.topology.connectivity(msh.topology.dim,0).links(cell_idx) # np.array
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 \
+ Dt*rho*inner(grad(u)*u_n,v)*dx + Dt*(mu * inner(grad(u), grad(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.
Any advice would be great.
Thank you in advance.
Giorgia
