[BUG]: no slip boundary condition is not successfully applied

[ruf10]

How to reproduce the bug

Hi guys,

I have been working on the no slip boundary condition using dolfinx for a long time. I don't think the no slip boundary condition is successfully applied.

Here is the setting I have: I am running a concentric 3D cylinders, inner cylinder with radius r and outer cylinder with radius R. I am trying to spin the inner cylinder with angular velocity 4 but gradually turning it on, and the outer cylinder is stay still. Also when I creating the mesh I use periodic boundary condition in Z direction/ the height in gmsh by the translation matrix.

I followed https://jsdokken.com/dolfinx-tutorial/chapter2/ns_code2.html#boundary-conditions carefully but when I see the velocity, it looks like near the outer cylinder the speed is faster. What I would like to see is inner is spinning and so the near inner boundary layer will be exactly the speed of my cylinder and outer near outer cylinder is 0.

I tried so many ways to apply the no slip boundary condition with different tutorials provided by dolfinx/ fenicsx, non of them return me the plausible result.

Could anyone take a look at it sooner?

Thanks!

Minimal Example (Python)

import numpy as np
import warnings
import dolfinx
from dolfinx.io import XDMFFile, gmshio
from mpi4py import MPI
from petsc4py import PETSc
import gmsh
import os
import ufl

#Cylinder setting
height = 3
R = 1
r = 0.5

dt = .01 #timestep size
nu = 0.001 
mu = .55
eps = 1e-6
tau = 1

#mesh generation with gmsh
#Create "taylor_couette" mesh
gmsh.initialize()
gmsh.model.add("taylor_couette")

if mesh_comm.rank == model_rank:
    cyl_in = gmsh.model.occ.addCylinder(0, 0, 0, 0, 0, height, r)
    cyl_out = gmsh.model.occ.addCylinder(0, 0, 0, 0,0, height, R)
if mesh_comm.rank == model_rank:
    fluid = gmsh.model.occ.cut([(3, cyl_out)], [(3, cyl_in)])
    gmsh.model.occ.synchronize()

fluid_marker = 1
if mesh_comm.rank == model_rank:
    volumes = gmsh.model.getEntities(dim=3)
    assert(len(volumes) == 1)#assert(volumes == fluid[0])
    gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
    gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid")
top_marker, bottom_marker, in_marker, out_marker =  2, 3, 4, 5
top_surface, bottom_surface, inner_surface, outer_surface  = [], [], [], []
if mesh_comm.rank == model_rank:
    boundaries = gmsh.model.getBoundary(volumes, oriented = False)
    for boundary in boundaries:
        center_of_mass = gmsh.model.occ.getCenterOfMass(boundary[0],boundary[1])
        mass = gmsh.model.occ.getMass(boundary[0], boundary[1])
        if np.allclose(center_of_mass,[0,0,height]):
            top_surface.append(boundary[1])
        elif np.allclose(center_of_mass,[0,0,0]):
            bottom_surface.append(boundary[1])
        elif np.allclose(mass, 2*np.pi*r*height):
            inner_surface.append(boundary[1])
        elif np.allclose(mass, 2*np.pi*R*height):
            outer_surface.append(boundary[1])
    gmsh.model.addPhysicalGroup(2, top_surface, top_marker)
    gmsh.model.setPhysicalName(2, top_marker, 'top surface')    
    gmsh.model.addPhysicalGroup(2, bottom_surface, bottom_marker)
    gmsh.model.setPhysicalName(2, bottom_marker, 'bottom surface')
    gmsh.model.addPhysicalGroup(2, inner_surface, in_marker)
    gmsh.model.setPhysicalName(2, in_marker, 'inner surface')
    gmsh.model.addPhysicalGroup(2, outer_surface, out_marker)
    gmsh.model.setPhysicalName(2, out_marker, 'out surface')
        # Define the periodic boundary condition
    # Define the translation vector
    translation_vector = [0., 0., height]
    # Create the affine transformation matrix
    affine_matrix = [1., 0., 0., translation_vector[0],
                    0., 1., 0., translation_vector[1],
                    0., 0., 1., translation_vector[2],
                    0, 0., 0., 1.]
    gmsh.model.mesh.setPeriodic(2, [top_marker], [bottom_marker], affine_matrix)
    tagMaster=gmsh.model.mesh.getPeriodic(2,[top_marker])
    print('tag Master')
    print(tagMaster) #checked the tag master is the bottom surface
if mesh_comm.rank == model_rank:
    inner_dist = gmsh.model.mesh.field.add("Distance")
    gmsh.model.mesh.field.setNumbers(inner_dist, "FacesList",inner_surface)
    resolution = r/8
    threshold = gmsh.model.mesh.field.add("Threshold")
    gmsh.model.mesh.field.setNumber(threshold, "IField", inner_dist)
    gmsh.model.mesh.field.setNumber(threshold, "LcMin", resolution)
    gmsh.model.mesh.field.setNumber(threshold, "LcMax", 3*resolution)
    gmsh.model.mesh.field.setNumber(threshold, "DistMin", 0.05*r)
    gmsh.model.mesh.field.setNumber(threshold, "DistMax", 0.1*r)

    outer_dist = gmsh.model.mesh.field.add("Distance")
    gmsh.model.mesh.field.setNumbers(outer_dist, "FacesList", outer_surface)
    outer_thre = gmsh.model.mesh.field.add("Threshold")
    gmsh.model.mesh.field.setNumber(outer_thre, "IField", outer_dist)
    gmsh.model.mesh.field.setNumber(outer_thre, "LcMin", resolution)
    gmsh.model.mesh.field.setNumber(outer_thre, "LcMax", 3*resolution)
    gmsh.model.mesh.field.setNumber(outer_thre, "DistMin", 0.05*r)
    gmsh.model.mesh.field.setNumber(outer_thre, "DistMax", 0.1*r)

    minimum = gmsh.model.mesh.field.add("Min")
    gmsh.model.mesh.field.setNumbers(minimum, "FieldsList", [threshold, outer_thre])
    gmsh.model.mesh.field.setAsBackgroundMesh(minimum)

    gmsh.model.occ.synchronize()
    gmsh.model.mesh.generate(3)
    gmsh.write("mesh3D_july1.msh")

msh, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim =3)
msh.name = "taylor_couette"
cell_markers.name = f"{msh.name}_cells"
facet_markers.name = f"{msh.name}_facets"

#2703 
from dolfinx import fem
#Taylor Hood elemFunctionnt 
P2 = ufl.VectorElement("CG", msh.ufl_cell(),2)
P1 = ufl.FiniteElement("CG", msh.ufl_cell(),1)
V = fem.FunctionSpace(msh, P2)
Q = fem.FunctionSpace(msh, P1)
K = fem.FunctionSpace(msh,P1) 

mixed = ufl.MixedElement([P2,P1])
W = fem.FunctionSpace(msh, mixed)
u, p = ufl.TrialFunctions(W)
v, q = ufl.TestFunctions(W)
k = ufl.TrialFunction(K)
phi = ufl.TestFunction(K)

#Solution vectors
w_ = fem.Function(W)
wnMinus1 = fem.Function(W)
unMinus1, pnMinus1 = wnMinus1.split() #time n
wn = fem.Function(W)
un,pn = wn.split()  #time n
wnPlus1 = fem.Function(W)
unPlus1, pnPlus1 = wnPlus1.split() #time n+1

fdim = msh.topology.dim - 1
msh.topology.create_connectivity(fdim, msh.topology.dim)
from dolfinx import mesh
boundary_facets = mesh.exterior_facet_indices(msh.topology)

#Define boundary conditions
#spin inner cylinder
class No_slip_inner:
    def __init__(self, mint):
        self.mint = mint
    def __call__(self, x):
        # Added some spatial variation here. 
        values = np.zeros((3, x.shape[1]),dtype=PETSc.ScalarType)
        values[0] =  self.mint* omega_inner*x[1]
        values[1] =  self.mint* omega_inner*(-x[0])
        return values

def smooth_bridge(t): 
    s = t/t_fullspeed
    #Smoothly increase from 0 at t=0 to 1 at t=x
    if s==0:
        return 0.0
    if s>=1:
        return 1.0
    else:
        return ufl.exp(-ufl.exp(-1./(1-s)**2)/s**2)#boundary condition for velocity
u_inner_spin = fem.Function(V)
inner_spin_velocity = No_slip_inner(mint)
u_inner_spin.interpolate(inner_spin_velocity)
bcu_inner = fem.dirichletbc(u_inner_spin, fem.locate_dofs_topological(V, fdim, facet_markers.find(in_marker)))
#outer cylinders
u_outer_nonslip = np.array((0,) * msh.geometry.dim, dtype=PETSc.ScalarType)
bcu_outer = fem.dirichletbc(u_outer_nonslip, fem.locate_dofs_topological(V, fdim,facet_markers.find(out_marker)), V)
bcu = [bcu_inner, bcu_outer] #boundary condtion for velocity

def update_bcu(t):
    mint = smooth_bridge(t)
    inner_spin_velocity = No_slip_inner(mint)
    u_inner_spin.interpolate(inner_spin_velocity)
    bcu_inner = fem.dirichletbc(u_inner_spin, fem.locate_dofs_topological(V, fdim, facet_markers.find(in_marker)))
    bcu = [bcu_inner, bcu_outer] #boundary condtion for velocity
    return bcu 

#Weak form
from ufl import  dx, nabla_grad, div, inner, dot, lhs, rhs, curl, cross

def aa(u,v):
     return inner(nabla_grad(u),nabla_grad(v))
def a_sym(u,v):
    return inner(.5*(nabla_grad(u)+nabla_grad(u).T),nabla_grad(v))
def b(u,v,w):
    return .5*(inner(dot(u,nabla_grad(v)),w)-inner(dot(u,nabla_grad(w)),v))
def convect(u,v,w):
    return dot(dot(u, nabla_grad(v)), w)
def c(p,v):
    return inner(p,div(v))
def grad_sym(u):
    return 0.5 * (nabla_grad(u) + nabla_grad(u).T)
def frobenius_norm_squared(T):
    return ufl.inner(T, T)

dx = ufl.dx(domain=msh)
x = ufl.SpatialCoordinate(msh) 

from ufl import sqrt, algebra
#Wall normal distance
func_wall_n_dist = - algebra.Abs(sqrt(x[0]**2 + x[1]**2)-(R+r)/2)+(R-r)/2

############################################################
#Set up the Krylov Solver
from petsc4py import PETSc
MAX_ITERS = 10000
solver = PETSc.KSP().create(msh.comm)
solver.setType(PETSc.KSP.Type.GMRES)
solver.setTolerances(rtol=1e-6,max_it=MAX_ITERS)
# Set Hypre preconditioner
pc = solver.getPC()
pc.setType(PETSc.PC.Type.HYPRE)
pc.setHYPREType("boomeramg")
#set solve for k
solver_k = PETSc.KSP().create(msh.comm)
solver_k.setType(PETSc.KSP.Type.GMRES)
solver_k.setTolerances(rtol=1e-6,max_it=MAX_ITERS)
# Set Hypre preconditioner
pc = solver_k.getPC()
pc.setType(PETSc.PC.Type.HYPRE)
pc.setHYPREType("boomeramg")

def solve_u(a, L, t):
        bcu = update_bcu(t)
        A = fem.petsc.assemble_matrix(a, bcu)
        A.assemble()
        b= fem.petsc.create_vector(L)  
        with b.localForm() as loc_bu:
            loc_bu.set(0)  
        fem.petsc.apply_lifting(b, [a], bcs=[bcu])
        b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES,mode =PETSc.ScatterMode.REVERSE)
        fem.petsc.set_bc(b, bcu)
        solver.setOperators(A)
        solver.solve(b, w_.vector)
        w_.x.scatter_forward()
        unPlus1, pnPlus1 = w_.split()
        for i in range(len(unPlus1.x.array[:])):
            unPlus1.x.array[i] = unPlus1.x.array[i] -(1./3.)*(unPlus1.x.array[i] - 2*un.x.array[i] + unMinus1.x.array[i])
        unPlus1.x.scatter_forward()
        return unPlus1

#Let's calculate
for i in range(len(unPlus1.x.array[:])):
    unMinus1.x.array[i] = 0.0
    un.x.array[i]= 0.0
t = t_init
count = 0
for  jj in range (0, k_on): #Before cylinder is up to speed, only NSE  
    t = t + dt
    nu_t = 0 #initilize nu_t
    count += 1 
    a = form((1./dt)*inner(u,v)* dx + b(2*un-unMinus1,u,v)*dx + 2.*nu*a_sym(u,v)* dx -c(p,v)*dx + c(q,u)* dx)
    L = form((1./dt)*(inner(un,v))*dx)
    unPlus1 = solve_u(a, L, t)
    #update solutions
    for i in range(len(unPlus1.x.array[:])):
        unMinus1.x.array[i] = un.x.array[i]
        un.x.array[i] = unPlus1.x.array[i]
    unMinus1.x.scatter_forward()
    un.x.scatter_forward()

Output (Python)

No response

Version

main branch

DOLFINx git commit

No response

Installation

No response

Additional information

No response

[francesco-ballarin]

This question is better suited for the discourse forum at https://fenicsproject.discourse.group/ , as from the code it is not immediately clear that this an actual bug in the dolfinx codebase. We can still reopen in future if, from the discussion there, this will turn out to actually be a bug.

[jorgensd]

@ruf10 I would suggest you debug your BC by using something along the lines of:

bcu = [bcu_inner, bcu_outer]
wh = fem.Function(W)
[fem.petsc.set_bc(wh.vector, bc) for bc in bcu]
with dolfinx.io.VTXWriter(msh.comm, "u.bp", [wh.sub(0).collapse()]) as vtx:
    vtx.write(0.0)

and visualize u in Paraview, to see if it sets the bcs as you would expect