MPC - weird features on periodic boundary of anisotropic mesh

dolfinx
1

Dear Community,

I am inverting a blockmatrix (stemming from a transport problem) on a cuboidal mesh with periodic boundary conditions along the z-axis. The latter are enforced using the MPC package. Details:

  1. the domain is a box with the following measurements: (0.5,1)\times (0,0.5) \times (0,2)
  2. I use an anisotropic mesh, with a resolution of 150x150 elements in the x-y-plane, and 5 elements along the z-direction.

I use a direct solver, and a smooth RHS/initial value. The solution has strange features and wiggles, that are at the size of a single element, and which I would regard as nonphysical. Here is a screenshot:

image
image1389×772 202 KB

Here is a minimal example code, that reproduces the problem (I ran it on dolfinx 0.9.0):

import numpy as np
import ufl
from dolfinx.fem import (Function, dirichletbc, form, 
                         functionspace, locate_dofs_topological)
from dolfinx.mesh import *
from dolfinx.mesh import CellType
from ufl import dx, grad, inner
from mpi4py import MPI
from petsc4py import PETSc
from ufl import as_vector
from dolfinx_mpc import (MultiPointConstraint,assemble_vector_nest,
                         create_vector_nest,assemble_matrix_nest,
                         create_matrix_nest)


# mesh
Nx  = 150
Ny  = 150
Nz  = 5
Lx  = 0.5 
Lxx = 1.0
Ly  = 0.5
Lz  = 2
msh = create_box(MPI.COMM_WORLD, 
                    [[Lx, 0.0, 0.0], [Lxx, Ly, Lz]], 
                    [Nx,Ny,Nz],
                    cell_type=CellType.hexahedron)
x = ufl.SpatialCoordinate(msh)
Xp = functionspace(msh, ("Lagrange", 1))
Xu = functionspace(msh, ("Lagrange", 1))
offset = Xu.dofmap.index_map.size_local * Xu.dofmap.index_map_bs 
u = ufl.TrialFunction(Xu)
p = ufl.TrialFunction(Xp)
w = ufl.TestFunction(Xu)
q = ufl.TestFunction(Xp)
u_func = Function(Xu, name="u_func")
p_func = Function(Xp, name="p_func")

# vector field
def vec_field(x,y,z):
    eps = 0.4    
    vx = eps*ufl.sin(5*z)
    vy = eps*ufl.cos(5*z)
    vz = np.sqrt(1-eps**2)
    return (vx,vy,vz)
vfield = as_vector(np.array(vec_field(x[0],x[1],x[2])))

# bc
boundary_dim = (msh).topology.dim - 1
facets = locate_entities_boundary(msh, dim=boundary_dim, 
                                    marker=lambda x: np.logical_or.reduce((
                                        np.isclose(x[0], Lx),
                                        np.isclose(x[0], Lxx),
                                        np.isclose(x[1], 0.0),
                                        np.isclose(x[1], Ly))))
dofs = locate_dofs_topological(V=Xp, entity_dim=boundary_dim, entities=facets)
BCs = [dirichletbc(0.0, dofs=dofs, V=Xp)]
def periodic_boundary(x):
    return np.isclose(x[2], Lz)                
def periodic_relation(x):
    out_x = np.zeros_like(x)
    out_x[0] = x[0]
    out_x[1] = x[1]
    out_x[2] = 0
    return out_x
mpc_u = MultiPointConstraint(Xu)
mpc_u.create_periodic_constraint_geometrical(Xu, periodic_boundary, periodic_relation, BCs)
mpc_u.finalize()
mpc_p = MultiPointConstraint(Xp)
mpc_p.create_periodic_constraint_geometrical(Xp, periodic_boundary, periodic_relation, BCs)
mpc_p.finalize()

# ic
u_func.interpolate(lambda x: np.sin((0.5/Ly)*2*np.pi*x[1])\
                          * np.cos((1/Lz)*2*np.pi*x[2])\
                          + np.sin((2/Ly)*2*np.pi*x[1])\
                          * np.cos((1/Lz)*2*np.pi*x[2]))
p_func.interpolate(lambda x: 0*x[0])

# weak form
M = 0.0005*inner(u,w)*dx
E = -1 * inner(inner(vfield,grad(p)), w)*dx
N = -1 * inner(grad(p),grad(q))*dx
F = -1 * inner(inner(vfield,grad(q)), u)*dx
a = form([[M, E],[F, N]])
l = form([0.0005*inner(u_func,w)*dx \
          +1 * inner(inner(vfield,grad(p_func)), w)*dx,
          -1 * inner(grad(p_func), grad(q))*dx \
          +1 * inner(u_func,inner(vfield,grad(q)))*dx 
          ])

# assembly
A = create_matrix_nest(a, [mpc_u, mpc_p])
assemble_matrix_nest(A, a, [mpc_u, mpc_p], BCs)
A.assemble()
L = create_vector_nest(l,[mpc_u, mpc_p])
up_vec = L.copy()
assemble_vector_nest(L,l,[mpc_u, mpc_p])

# solve
solver = PETSc.KSP().create(msh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
solver.solve(L, up_vec)

# plot solution
from dolfinx import plot
import pyvista
u_func.x.array[:] = up_vec.array_r[:offset]
mpc_u.backsubstitution(u_func)
cells, types, x = plot.vtk_mesh(Xu)
grid = pyvista.UnstructuredGrid(cells, types, x)
grid.point_data["u"] = u_func.x.array.real
grid.set_active_scalars("u")
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.show_axes()
plotter.show()

I noticed three things:

  1. The wiggles appear only on (or come from) the periodic boundary
  2. These wiggles only occur when i choose an anisotropic mesh (does MPC somehow fail to match the DOFs from the too ends of the domain, if the mesh is two anisotropic?)
  3. If i translate the domain along the x-axis towards the origin, i.e. instead of (0.5,1)\times (0,0.5) \times (0,2), i use
    (0,0.5)\times (0,0.5) \times (0,2), the wiggles also do not appear anymore.

Any advice is highly appreciated