Issues with time-dependent Dirichlet boundary conditions in mixed element space (dolfinx)

General
1

Dear all,

I’m trying to solve Navier Stokes equations in a 2D toy problem with time-dependent Dirichlet inflow boundary conditions using a monolithic solver. The code I created does not produce any error but when it comes to the velocity solution (at t_end for instance) it seems like the inflow boundary condition is not applied correctly to all the nodes, as you can see in the following picture:

picture1
picture11654×908 263 KB

The same phenomenon happens also when refining the mesh on what seem to be random nodes on the boundary. I ran the code in serial, if this may be of any interest.
Here you can find the minimal working example that produces this result:

import dolfinx as df
from mpi4py import MPI
from petsc4py import PETSc
import ufl
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, div, grad, dot, ds, dx, inner, lhs, nabla_grad, rhs, sym)

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

# set fluid parameters
mu = 0.035
rho = 1.0
Re = 1200.

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

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

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

# 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], self.t*20))) 

V0, _ = W.sub(0).collapse()
inflow_expr = time_dependent_inflow()
inflow_expr.t = 0.
bc_fun, _ = Function(W).split()
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))

# Set noslip BC
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)

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

u, p = TrialFunctions(W)
v, q = TestFunctions(W)
n = FacetNormal(msh)
h = ufl.CellDiameter(msh)

problem_type = 'NS'

if problem_type == 'NS': 
        
    a = rho*inner( u,v )*dx - Dt*( - rho*inner( grad(u)*u_n, v )*dx - mu*inner( grad(u),grad(v) )*dx + p*div(v)*dx - q*div(u)*dx) 
    L = rho*inner( u_n,v )*dx

    # add backflow stabilization
    boundary_id = [3] 
    for el in boundary_id:
        a += 0.25*rho*Dt*( abs( inner( u_n,n ) - abs(inner( u_n,n )) )*inner( u,v )*ds(el) )

# 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": "mumps"})
    w_n = problem.solve()
    u_n, p_n = w_n.split()

u_n.name = 'u'

# write function for visualization
filepath = './u_n.xdmf'      
with df.io.XDMFFile(msh.comm, filepath, 'w') as xdmf:
    u_n.x.scatter_forward()
    xdmf.write_mesh(msh) 
    xdmf.write_function(u_n)

Any advice on how to fix it would be great.
Thanks in advance!

2

I think this should be bc_fun = Function(V0), as you have done for u_bc.

3

Thank you, you are right, that small modification solved the issue.