I was recently working to get a solenoid model working in dolfinx using the magnetostatic equations from the magnetostatics example of the FEniCs tutorial. I was able to get the example working in two dimensions, however adapting it to 3D has led to some issues. Although I am able to reproduce theoretical results for the B field in some examples (the Helmholtz Coil) I have noticed discontinuities in my solution for the vector potential (\vec{A}) that do not appear when I run similar code in dolfin.
To test this, I made a simple solenoid mesh using pygmsh and meshio and used the same mesh to solve for A_x in both dolfin and dolfinx. Comparing the results in paraview, the peak values are different (1.5e-7 vs 2.2e-7) and the dolfinx solution shows noticeable visual artifacts that are also visible clearly when I do a Plot Over Line in paraview. I have pasted two images below comparing the two solutions. I used the stable:current docker image for dolfin and the latest dolfinx docker image.
Here is the dolfin solution image:
And here is the dolfinx solution image:
Here is the MWE code for mesh generation, dolfin, and dolfinx solutions:
Mesh Generation:
from pygmsh.opencascade import Geometry
from pygmsh import generate_mesh
import meshio
import numpy as np
domain_radius = 20.0
length = 15.0
rect1_r = 2.0
rect2_r = .5
geom = Geometry()
rect1 = geom.add_box([-rect1_r/2,-rect1_r/2,-length/2], [rect1_r,rect1_r,length], char_length=0.2)
rect2 = geom.add_box([-rect2_r/2,-rect2_r/2,-length/2], [rect2_r,rect2_r,length], char_length=0.2)
solenoid = geom.boolean_difference([rect1], [rect2])
sphere = geom.add_ball([0,0,0], domain_radius, char_length = 4.0)
sphere_domain = geom.boolean_difference([sphere], [solenoid], delete_other=False)
geom.boolean_fragments([solenoid], [sphere_domain], delete_first = False, delete_other=False)
geom.add_raw_code("Physical Volume(1) = {{{}}};".format(solenoid.id)) # Torus 1
geom.add_raw_code("Physical Volume(2) = {{{}}};".format(sphere_domain.id))
geom.add_raw_code("Physical Surface(3) = {{{}}};".format("17")) # Sphere border
# Build the .geo file.
mesh = generate_mesh(geom, geo_filename="mesh.geo", )
# Mesh Conversion
tetra_cells = None
tetra_data = None
triangle_cells = None
triangle_data = None
print(mesh.cell_data_dict["gmsh:physical"].keys())
for key in mesh.cell_data_dict["gmsh:physical"].keys():
if key == "triangle":
triangle_data = mesh.cell_data_dict["gmsh:physical"][key]
elif key == "tetra":
tetra_data = mesh.cell_data_dict["gmsh:physical"][key]
for cell in mesh.cells:
if cell.type == "tetra":
if tetra_cells is None:
tetra_cells = cell.data
else:
tetra_cells = np.vstack([tetra_cells, cell.data])
elif cell.type == "triangle":
if triangle_cells is None:
triangle_cells = cell.data
else:
triangle_cells = np.vstack([triangle_cells, cell.data])
tetra_mesh = meshio.Mesh(points=mesh.points,
cells={"tetra": tetra_cells},
cell_data={"name_to_read":[tetra_data]})
triangle_mesh =meshio.Mesh(points=mesh.points,
cells=[("triangle", triangle_cells)],
cell_data={"name_to_read":[triangle_data]})
meshio.xdmf.write("tetra_mesh.xdmf", tetra_mesh)
meshio.xdmf.write("triangle_mesh.xdmf", triangle_mesh)
Dolfin solution:
from fenics import *
from mshr import *
from mpi4py import MPI
import dolfin
import numpy as np
r = 1 #Approx radius
J = 1 #Current density in A/m^2
mesh = Mesh()
with XDMFFile("tetra_mesh.xdmf") as infile:
infile.read(mesh)
mvc_tetra = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("tetra_mesh.xdmf") as infile:
infile.read(mvc_tetra, "name_to_read")
markers = dolfin.cpp.mesh.MeshFunctionSizet(mesh, mvc_tetra)
mvc_triangle = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("triangle_mesh.xdmf") as infile:
infile.read(mvc_triangle, "name_to_read")
mf_triangle = dolfin.cpp.mesh.MeshFunctionSizet(mesh, mvc_triangle)
# Define function space and integration measure
V = FunctionSpace(mesh, "P", 1)
dx = Measure('dx', domain=mesh, subdomain_data=markers)
# Define Dirichlet boundary condition
bc = DirichletBC(V, Constant(0), mf_triangle, 3)
# define current densities as function of subdomains
class CurrentDensity_x(UserExpression):
def __init__(self, markers, **kwargs):
super().__init__(**kwargs)
self.markers = markers
def eval_cell(self, values, x, cell):
if (self.markers[cell.index] == 1):
values[0] = -x[1]/r*J
else:
values[0] = Constant(0.0)
mu = Constant(4*np.pi*10**-7)
# initialise current densities
Jx = CurrentDensity_x(markers, degree = 1)
# Define Variational Problem
Ax = TrialFunction(V)
vx = TestFunction(V)
# left hand sides of variational problems
ax = (1/mu)*dot(grad(Ax), grad(vx))*dx
# right hand sides of variational problems
Lx = Jx*vx*dx
# solve variational problems for magnetic Ax (in Vs/mm)
Ax = Function(V)
solve(ax == Lx, Ax, bc)
with XDMFFile(MPI.COMM_WORLD,"A_x_dolfin.xdmf") as outfile:
outfile.write(mesh)
outfile.write_checkpoint(Ax, 'Ax')
Dolfinx Solution:
import dolfinx
import numpy as np
from dolfinx import Mesh, geometry, Function, DirichletBC, Constant, solve, FunctionSpace, VectorFunctionSpace
from dolfinx.io import XDMFFile
import ufl
from ufl import TrialFunction, TestFunction, dot, grad, dx, inner, Measure
from mpi4py import MPI
r = 1 #Approx radius, coil
J = 1 #Current density in A/m^2
with XDMFFile(MPI.COMM_WORLD,"tetra_mesh.xdmf",'r') as infile:
mesh = infile.read_mesh(dolfinx.cpp.mesh.GhostMode.none, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
with XDMFFile(MPI.COMM_WORLD, "tetra_mesh.xdmf", "r") as xdmf_infile:
mt = xdmf_infile.read_meshtags(mesh, name="Grid")
with XDMFFile(MPI.COMM_WORLD, "triangle_mesh.xdmf", "r") as xdmf_infile:
mt_triangle = xdmf_infile.read_meshtags(mesh, name="Grid")
V = FunctionSpace(mesh, ("CG", 1))
# Define Dirichlet boundary condition
outer_facets = mt_triangle.indices[mt_triangle.values == 3]
dofs_outer = dolfinx.fem.locate_dofs_topological(V, 1, outer_facets)
ubcOuter = Function(V)
ubcOuter.vector.set(0)
bc = DirichletBC(ubcOuter, dofs_outer)
# Define subdomain markers and integration measure
dx = Measure('dx', domain=mesh, subdomain_data=mt)
#Define Current Density and Mu0
x = ufl.SpatialCoordinate(mesh)
J_x = -x[1]/r*J
mu = Constant(mesh, 4*np.pi*10**-7)
#Define Variational Problem
A_x = TrialFunction(V)
vx = TestFunction(V)
ax = (1 / mu)*inner(grad(A_x), grad(vx))*dx
Lx = inner(J_x ,vx)*dx(1)
#Solve Variational Problem
A_x = Function(V)
solve(ax == Lx, A_x, bc)
#Store solution
with XDMFFile(MPI.COMM_WORLD, "A_x_dolfinx.xdmf", "w") as file:
file.write_mesh(mesh)
file.write_function(A_x)
I have tried experimenting with different preconditioners and ksp types but have not seen any differences. I have also tried many different methods to assign J in dolfinx thinking that it might be an issue with the discontinuous nature of J and found the way shown above to be the shortest code-wise but other methods (vector.set_local(i,___), etc) all yielded the same result. These issues with A_x yield further issues when taking the curl later to calculate the magnetic field, yielding even more numerical discontinuities. Are numerical discontinuities a known issue in Dolfinx, and what steps can I take to mitigate this?