I am solving Laplace’s equation \Delta \phi = 0 on a box. The boundaries have values of 0 for the front y-face, back y-face, and bottom. \phi = 1 on the top. The x faces are periodic. The bottom face of the box is located at z = 0.
When the box has a height of 2, I see an edge artifact (negative value ~8.3e-5) from a single point on the edge parallel with the z axis, at ~z=1.5. It is on the upper, Dirichlet y-face, next to the lower, periodic x face.
However, when I run with a height zlen = 1, the artifact disappears and the lowest global value of potential is 0. I am using dolfin-x version 0.9.0 with conda. The program is ran with 1 processor.
# This demo solves Laplace's Eqn on a simple box
# - top potential = 1
# - bottom potential = 0
# - front and back y faces = 0
# - x faces are periodic
from mpi4py import MPI
from petsc4py import PETSc
import dolfinx.fem as fem
import numpy as np
from dolfinx.fem import Function
from dolfinx.io import XDMFFile, gmshio
from petsc4py.PETSc import ScalarType # type: ignore
from ufl import (
TestFunction,
TrialFunction,
dx,
grad,
inner,
Measure,
)
import dolfinx_mpc
import gmsh
def create_box(xlen, ylen, zlen):
# generate mesh
gmsh.initialize()
# name the model
gmsh.model.add("square")
# Create box using occ
# bottom centered on (0, 0, 0)
# top centered on (0, 0, zlen)
gmsh.model.occ.addBox(x=-0.5 * xlen, y=-0.5*ylen, z=0, dx=xlen, dy=ylen, dz=zlen)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.2)
gmsh.model.occ.synchronize()
return
def add_groups(xlen, ylen, zlen):
#add physical groups
top = []
bottom = []
sides = []
# add boundaries as physical groups
boundaries = gmsh.model.occ.getEntities(dim=2)
# loop over boundaries
for boundary in boundaries:
com = gmsh.model.occ.getCenterOfMass(boundary[0], boundary[1])
if np.allclose(com[2], 0):
bottom.append(boundary[1])
elif np.allclose(com[2], zlen):
top.append(boundary[1])
else:
sides.append(boundary[1])
for qty in [top, bottom, sides]:
#top added first, bottom added second
gmsh.model.addPhysicalGroup(2, qty)
# add entire domain to groups
layers = gmsh.model.occ.getEntities(dim=3)
gmsh.model.addPhysicalGroup(layers[0][0], [layers[0][1]])
gmsh.model.mesh.generate(3)
# save mesh to file
gmsh.write("square.msh")
gmsh.finalize()
return
def assign_periodic_mesh(xlen, ylen, zlen):
translation_x = [1, 0, 0, xlen,
0, 1, 0, 0,
0, 0, 1, 0,
0, 0, 0, 1]
eps = 1e-5
# Get left surface:
sxmin = gmsh.model.getEntitiesInBoundingBox(
xmin=-0.5 * xlen - eps,
ymin=-0.5 * ylen - eps,
zmin=-eps,
xmax=-0.5 * xlen + eps,
ymax=0.5 * ylen + eps,
zmax=zlen + eps,
dim=2,
)
for i_x in sxmin:
# Get the bounding box of the left surface
xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.getBoundingBox(
i_x[0], i_x[1]
)
# Translate the bounding box to the right and look for surfaces inside it:
sxmax = gmsh.model.getEntitiesInBoundingBox(
xmin - eps + xlen,
ymin - eps,
zmin - eps,
xmax + eps + xlen,
ymax + eps,
zmax + eps,
2,
)
# For all the matches, we compare the corresponding bounding boxes...
for j_x in sxmax:
xmin2, ymin2, zmin2, xmax2, ymax2, zmax2 = gmsh.model.getBoundingBox(
j_x[0], j_x[1]
)
xmin2 -= xlen
xmax2 -= xlen
# ...and if they match, we apply the periodicity constraint
if (
abs(xmin2 - xmin) < eps
and abs(xmax2 - xmax) < eps
and abs(ymin2 - ymin) < eps
and abs(ymax2 - ymax) < eps
and abs(zmin2 - zmin) < eps
and abs(zmax2 - zmax) < eps
):
print('x periodic success!')
gmsh.model.mesh.setPeriodic(1, [j_x[1]], [i_x[1]], translation_x)
return
def add_dirichlet_BCs(V, ylen, zlen):
# front y face
geo_dofs = fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[1], 0.5 * ylen))
bc_front_y = fem.dirichletbc(ScalarType(0.0), geo_dofs, V)
# back y face
geo_dofs = fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[1], -0.5 * ylen))
bc_back_y = fem.dirichletbc(ScalarType(0.0), geo_dofs, V)
# bottom
geo_dofs = fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[2], 0.0))
bc_low = fem.dirichletbc(ScalarType(0.0), geo_dofs, V)
# top
geo_dofs = fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[2], zlen))
bc_hi = fem.dirichletbc(ScalarType(1.0), geo_dofs, V)
bcs = [bc_low, bc_hi, bc_front_y, bc_back_y]
return bcs
def periodic_boundaries(bcs, xlen):
def periodic_boundary(x):
return np.isclose(x[0], 0.5 * xlen)
def periodic_relation(x):
out_x = np.zeros_like(x)
out_x[0] = x[0] - xlen
out_x[1] = x[1]
out_x[2] = x[2]
return out_x
mpc = dolfinx_mpc.MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(V, periodic_boundary, periodic_relation, bcs)
mpc.finalize()
return mpc
### Main ###
# create mesh
# dimensions of box
xlen = 1
ylen = 1
zlen = 2
rank = MPI.COMM_WORLD.rank
if rank == 0:
create_box(xlen, ylen, zlen)
assign_periodic_mesh(xlen, ylen, zlen)
add_groups(xlen, ylen, zlen)
# load in mesh
mesh, cell_markers, facet_markers = gmshio.read_from_msh(
"square.msh", MPI.COMM_WORLD, gdim=3
)
# begin BCs
tdim = mesh.geometry.dim
V = fem.functionspace(mesh, ("Lagrange", 1))
bcs = add_dirichlet_BCs(V, ylen, zlen)
mpc = periodic_boundaries(bcs, xlen)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
dx = Measure("dx", domain=mesh)
a = inner(grad(u), grad(v)) * dx
bilinear_form = fem.form(a)
# set homogeneous source term for Laplace Eqn.
f = fem.Constant(mesh, ScalarType(0.0))
L = inner(f, v) * dx
linear_form = fem.form(L)
A_org = dolfinx_mpc.assemble_matrix(bilinear_form, mpc, bcs=bcs)
A_org.assemble()
L_org = dolfinx_mpc.assemble_vector(linear_form, mpc)
dolfinx_mpc.apply_lifting(L_org, [bilinear_form], [bcs], mpc)
fem.petsc.set_bc(L_org, bcs)
# solver settings
ksp = PETSc.KSP().create(mesh.comm)
ksp.setOperators(A_org)
ksp.setType(PETSc.KSP.Type.PREONLY)
pc = ksp.getPC()
pc.setType("lu")
pc.setHYPREType("mumps")
pc.setUp()
solution = Function(mpc.function_space)
ksp.solve(L_org, solution.x.petsc_vec)
processors = MPI.COMM_WORLD.size
name = 'solution_%i_procs'%(processors)
print("Writing solution to", name)
with XDMFFile(MPI.COMM_WORLD, name + '.xdmf', "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(solution)
