3D magnetostatics - Curl-conform field

General
1

Dear all,

I have finally finished to build the 3D magnetostatics model; a coil with impressed current density surrounded by “air”.

I am using the pyvistaqt library to plot the 3D field.

Here is the code for those who want to play with and likely to improve it. Please share if you improve. In particular,I found rather difficult to converge (see the options for PETSC). I am using GMRES with a composite preconditioner that gave the “best” convergence.

# Install gmsh via "pip install gmsh"

import time
import tracemalloc

tracemalloc.start()
st = time.time()

import os, psutil

from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.io import gmshio 
from dolfinx.io import XDMFFile
from dolfinx.fem import (Expression, Function, FunctionSpace, VectorFunctionSpace)
from dolfinx.fem import (dirichletbc, locate_dofs_topological)
from dolfinx import *
from dolfinx import io

import ufl
from ufl import as_vector, dot, dx, grad, curl, inner, sqrt

from petsc4py.PETSc import ScalarType
import numpy as np

rank = MPI.COMM_WORLD.rank

import gmsh
import sys
import math

gmsh.initialize()

gmsh.model.add("magnetostatics3D")

gmsh.option.setColor("Mesh.Color.Lines", 0, 0, 0)

gmsh.option.setNumber('Mesh.Optimize', 1)
gmsh.option.setNumber('Mesh.OptimizeNetgen', 1)
gmsh.option.setNumber('Mesh.Algorithm3D', 10)

model_rank = 0
mesh_comm = MPI.COMM_WORLD
gdim = 3

flag_visu = 0
flag_tet = 0

Ri = 0.1
th = 0.05
Ra = 0.3
la = 0.5

nl = 30
lc = th / 4

if mesh_comm.rank == model_rank:
	cascade = gmsh.model.occ
	meshing = gmsh.model.mesh

	cascade.addPoint(Ri, 0, -th/2, lc, 1)
	cascade.addPoint(Ri+th, 0, -th/2, lc, 2)
	cascade.addPoint(Ri+th, 0, th/2, lc, 3)
	cascade.addPoint(Ri, 0, th/2, lc, 4)
	cascade.addLine(1, 2, 1)
	cascade.addLine(2, 3, 2)
	cascade.addLine(3, 4, 3)
	cascade.addLine(4, 1, 4)
	cascade.addCurveLoop([1, 2, 3, 4], 1)
	cascade.addPlaneSurface([1], 1)
	
	cascade.addPoint(0., 0., -la/2, lc/2, 5)
	cascade.addPoint(0., 0., la/2, lc/2, 6)
	cascade.addLine(5, 6, 5)
	
	cascade.synchronize()

	ls_1 = cascade.revolve([(gdim-1, 1)], 0, 0, 0, 0, 0, 1, -math.pi, [nl], recombine = False)
	ls_2 = cascade.revolve([(gdim-1, 1)], 0, 0, 0, 0, 0, 1, math.pi, [nl], recombine = False)
	
	if (flag_tet == 1):
		ls_3 = cascade.fuse([(gdim, 1)], [(gdim, 2)], 3)
		coil_tag = [3]
		coil_dimtag = [(gdim, 3)]
		air_tag = [4]
	else:
		coil_tag = [1, 2]
		coil_dimtag = [(gdim, 1), (gdim, 2)]
		air_tag = [3]
		
	cascade.fragment(gmsh.model.occ.getEntities(gdim), [])
	
	cascade.addSphere(0, 0, 0, Ra, 4)
	# Cut out the coil from the sphere keeping a unique interface for both
	coil_dimtag += [(gdim-2, 5)]
	
	cascade.fragment([(gdim, 4)], coil_dimtag)
	
	cascade.removeAllDuplicates
	
	cascade.synchronize()
	
	up, down = gmsh.model.getAdjacencies(gdim, air_tag[0])
	coilBoundary_dimtag = gmsh.model.getBoundary([(gdim, 3)], oriented=False)
	coilBoundary_tags = [tup[1] for tup in coilBoundary_dimtag]
	if (flag_tet == 0):
		del coilBoundary_tags[-1]

	gmsh.model.addPhysicalGroup(gdim, coil_tag, 1, name="Coil")
	gmsh.model.setColor(ls_1, 255, 0, 0)
	gmsh.model.setColor(ls_2, 255, 0, 0)
	gmsh.model.addPhysicalGroup(gdim, air_tag, 2, name="Air")
	gmsh.model.setColor([(gdim, 3)], 0, 0, 255)
	gmsh.model.addPhysicalGroup(gdim-1, [down[0]], 3, name="AirBoundary")
	gmsh.model.setColor([(gdim-1, down[0])], 0, 0, 255)
	gmsh.model.addPhysicalGroup(gdim-1, coilBoundary_tags, 4, name="CoilBoundary")
	gmsh.model.setColor(coilBoundary_dimtag, 0, 0, 255)
	gmsh.model.addPhysicalGroup(1, [5], 5, name="Axis")
	
	meshing.removeDuplicateNodes
	meshing.removeDuplicateElements

	meshing.generate(gdim)

	gmsh.write("magnetostatics3D.msh")

# Read in by dolfinx
mesh, cell_tags, facet_tags = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
gmsh.finalize()

print("### MESH: DONE ###")

# Output DOLFINx meshes to file
#with XDMFFile(MPI.COMM_WORLD, "meshCheck.xdmf", "w") as xdmf:
#    xdmf.write_mesh(mesh)
#    xdmf.write_meshtags(facet_tags)
#    xdmf.write_meshtags(cell_tags)


# Function space: edge elements
V0 = FunctionSpace(mesh, ("N1curl", 1))

# Create facet to cell connectivity required to determine boundary facets
dofs_facets = locate_dofs_topological(V0, mesh.topology.dim-1, facet_tags.find(3))

# Working boundary condition
u0 = Function(V0)
u0.x.set(0.)
bc = dirichletbc(u0, dofs_facets)

# Trial and test functions
u = ufl.TrialFunction(V0)
v = ufl.TestFunction(V0)

##Define Current Density (constant discontinuous)
V1 = VectorFunctionSpace(mesh, ("DG", 0))
J = Function(V1, dtype=ScalarType, name='Current_density')
cells1 = cell_tags.find(1)

J0 = 1e8
def Je(x, J0):
    return np.vstack([-J0*np.sin(np.arctan2(x[1],x[0]))*np.ones(x.shape[1]), J0*np.cos(np.arctan2(x[1],x[0]))*np.ones(x.shape[1]), np.zeros(x.shape[1])])

J.interpolate(lambda x: Je(x, J0), cells1)
J.x.scatter_forward()

nu = 1e7 / (4*np.pi)

a = nu * inner(curl(u), curl(v)) * dx
l = inner(J, v) * dx

A = Function(V0)
options = [ \
	{"ksp_type": "preonly", "pc_type": "lu", "ksp_monitor": None, "ksp_view": None}, \
	
	{"ksp_type": "gmres", "ksp_gmres_restart":30, "ksp_rtol":1e-6, "ksp_atol":1e-12, "ksp_max_it": 1000, "pc_type": "composite", "pc_composite_pcs": "sor,icc", "ksp_monitor": None, "ksp_monitor_true_residual": None, "ksp_view": None, "monitor_convergence": True}, \
	
	{"ksp_type": "gmres", "ksp_rtol":1e-6, "ksp_atol":1e-10, "ksp_max_it": 1500, "pc_type": "sor", "ksp_monitor": None, "ksp_view": None, "ksp_monitor_true_residual": None, "monitor_convergence": True}, \
	{"ksp_type": "lgmres", "ksp_rtol":1e-6, "ksp_atol":1e-10, "ksp_max_it": 1500, "pc_type": "ilu", "ksp_monitor": None, "ksp_view": None, "monitor_convergence": True}, \
	{"ksp_type": "pipefgmres", "ksp_rtol":1e-6, "ksp_atol":1e-10, "ksp_max_it": 1500, "pc_type": "ilu", "ksp_monitor": None, "ksp_view": None, "monitor_convergence": True}, \
	{"ksp_type": "bcgsl", "ksp_rtol":1e-6, "ksp_atol":1e-10, "ksp_max_it": 1000, "ksp_monitor": None, "ksp_view": None, "monitor_convergence": True}, \
	{"ksp_type": "ibcgs", "ksp_rtol":1e-6, "ksp_atol":1e-10, "ksp_max_it": 1000, "ksp_monitor": None, "ksp_view": None, "monitor_convergence": True}, \
	{"ksp_type": "cg", "ksp_rtol":1e-6, "ksp_atol":1e-10, "ksp_max_it": 1000, "ksp_monitor": None, "ksp_view": None, "monitor_convergence": True}, \
	{"ksp_type": "bicg", "ksp_rtol":1e-6, "ksp_atol":1e-10, "ksp_max_it": 1000, "ksp_monitor": None, "ksp_view": None, "monitor_convergence": True} \
	]
#"pc_factor_mat_solver_type": "umfpack", or mumps

problem = fem.petsc.LinearProblem(a, l, bcs = [bc], petsc_options = options[1])
A = problem.solve()

solver = problem.solver
solver.setFromOptions()
print("### SOLVE: DONE ###")
num_dofs_global = V0.dofmap.index_map.size_global * V0.dofmap.index_map_bs
print(f"Number of dofs: {num_dofs_global}")
print('Iteration number:', solver.getIterationNumber())
print('Iteration residual:', solver.getResidualNorm())

V2 = VectorFunctionSpace(mesh, ("DG", 1))
B = Function(V2, dtype=ScalarType, name='B_field')
B_expr = Expression(ufl.curl(A), V2.element.interpolation_points())
B.interpolate(B_expr)

V3 = FunctionSpace(mesh, ("DG", 1))
Bn = Function(V3, dtype=ScalarType, name='B_norm')
Bn_expr = Expression(ufl.sqrt(ufl.dot(B, B)), V3.element.interpolation_points())
Bn.interpolate(Bn_expr)

print("Maximum B field (T)", Bn.x.array.max())


print("### POSTPROCESS: DONE ###")

#with io.XDMFFile(mesh.comm, "J.xdmf", "w") as xdmf:
#    xdmf.write_mesh(mesh)
#    xdmf.write_meshtags(cell_tags) 
#    xdmf.write_function(J)

#with io.XDMFFile(mesh.comm, "A.xdmf", "w") as xdmf:
#    xdmf.write_mesh(mesh)
#    xdmf.write_meshtags(cell_tags) 
#    xdmf.write_function(A)
#    
#with io.XDMFFile(mesh.comm, "B.xdmf", "w") as xdmf:
#    xdmf.write_mesh(mesh)
#    xdmf.write_meshtags(cell_tags) 
#    xdmf.write_function(B)


nbp = 101
z = np.linspace(-la/2, la/2, nbp)
points = np.zeros((gdim, nbp))
points[2] = z
Bn_values = []

from dolfinx import geometry
bb_tree = geometry.BoundingBoxTree(mesh, mesh.topology.dim)

cells = []
points_on_proc = []
# Find cells whose bounding-box collide with the the points
cell_candidates = geometry.compute_collisions(bb_tree, points.T)
# Choose one of the cells that contains the point
colliding_cells = geometry.compute_colliding_cells(mesh, cell_candidates, points.T)
for i, point in enumerate(points.T):
	if len(colliding_cells.links(i)) > 0:
		points_on_proc.append(point)
		cells.append(colliding_cells.links(i)[0])

points_on_proc = np.array(points_on_proc, dtype=np.float64)
Bn_values = Bn.eval(points_on_proc, cells)

import matplotlib.pyplot as plt
fig = plt.figure()
plt.plot(points_on_proc[:,2], Bn_values, "r", linewidth=2, label="|$B$| on axis")
plt.grid(True)
plt.ylabel("|$B$| (T)")
plt.xlabel("$z$ (m)")
plt.legend()
plt.show(block=False)
# If run in parallel as a python file, we save a plot per processor
plt.savefig(f"BFieldOnAxis.png")
if (flag_tet == 1):
	np.savetxt('fenicsxTetra.dat', np.c_[points_on_proc[:,2], Bn_values])
else:
	np.savetxt('fenicsxHexaDiag.dat', np.c_[points_on_proc[:,2], Bn_values])

if (flag_visu == 1):
	import pyvista
	from pyvistaqt import BackgroundPlotter
	import os
	os.environ['PYVISTA_OFF_SCREEN'] = 'false'

	# If environment variable PYVISTA_OFF_SCREEN is set to true save a png
	# otherwise create interactive plot
	if pyvista.OFF_SCREEN:
		pyvista.start_xvfb(wait=0.1)

	# Set some global options for all plots
	transparent = True
	figsize_x = 1200
	figsize_y = 800
	pyvista.rcParams["background"] = [1., 1., 1.]

	# Create a plotter with two subplots, and add mesh tag plot to the
	# first sub-window
	subplotter = BackgroundPlotter(shape=(1, 1))

	subplotter.subplot(0, 0)
	subplotter.add_text("B field", font_size=14, color="black", position="upper_edge")
	# Next, we create a pyvista.UnstructuredGrid based on the mesh
	mesh_cells, cell_types, x = plot.create_vtk_mesh(mesh)
	gridmesh = pyvista.UnstructuredGrid(mesh_cells, cell_types, x)
	# Attach the cells tag data to the pyvita grid
	gridmesh.cell_data["Marker"] = cell_tags.values
	gridmesh.set_active_scalars("Marker")
	# Next, we create a pyvista.UnstructuredGrid based on the mesh
	cells, cell_types, x = plot.create_vtk_mesh(V2)
	grid = pyvista.UnstructuredGrid(cells, cell_types, x)
	# Add this grid (as a wireframe) to the plotter
	subplotter.add_mesh(gridmesh, style="wireframe", line_width=2, color="black")
	# Create point cloud of vertices, and add the vertex values to the cloud
	grid.point_data["B"] = B.x.array.reshape(x.shape[0], V2.dofmap.index_map_bs)
	glyphs1 = grid.glyph(orient = "B", factor = 0.04)
	sargs = dict(
    title_font_size=26,
    label_font_size=20,
    shadow=False,
    n_labels=6,
    italic=True,
    fmt="%.2f",
    font_family="arial",
    title='B (T)',
    interactive=True,
    color="black")
	subplotter.add_mesh(glyphs1, scalar_bar_args=sargs)
	
#	subplotter.subplot(0, 1)
#	subplotter.add_text("B norm on coil", font_size=14, color="black", position="upper_edge")
#	cells, types, x = plot.create_vtk_mesh(mesh, entities=cells1)
#	sub_mesh = pyvista.UnstructuredGrid(cells, types, x)
#	subplotter.add_mesh(sub_mesh, style="wireframe", line_width=1, color="gray")

#	facets1 = facet_tags.find(4)
#	cells, cell_types, x = plot.create_vtk_mesh(V3, entities=facets1)
#	grid2 = pyvista.UnstructuredGrid(cells, cell_types, x)
#	grid2.point_data["Bn"] = Bn.x.array
#	grid2.set_active_scalars("Bn")
#	subplotter.add_mesh(grid2, scalar_bar_args={'title': '|B| (T)'}, scalars="Bn")
#	subplotter.view_xy()

#Does not plot (DG, 0)
#	subplotter.subplot(0, 2)
#	subplotter.add_text("Current density in coil", font_size=14, color="black", position="upper_edge")
#	cells, types, x = plot.create_vtk_mesh(mesh, entities=cell_tags.find(1))
#	sub_mesh = pyvista.UnstructuredGrid(cells, types, x)
#	subplotter.add_mesh(sub_mesh, style="wireframe", line_width=2, color="black")

#	cells, cell_types, x = plot.create_vtk_mesh(V1)
#	grid3 = pyvista.UnstructuredGrid(cells, cell_types, x)
#	grid3.point_data["Je"] = J.x.array.reshape(x.shape[0], V1.dofmap.index_map_bs)
#	glyphs3 = grid3.glyph(orient = "Je", factor = 0.5e-9)
#	subplotter.add_mesh(glyphs3, scalar_bar_args={'title': 'Je (A/m²)'})

	# displaying the memory
	print("Current and peak memory (Mb)", 1e-6*tracemalloc.get_traced_memory()[0], 1e-6*tracemalloc.get_traced_memory()[1])
 
# stopping the library
tracemalloc.stop()

process = psutil.Process(os.getpid())
print("Total memory (Gb)", 1e-9*process.memory_info().rss)  # in bytes 

et = time.time()
print("CPU time (s) =", et-st)

if (flag_visu == 1):
	# Save as png if we are using a container with no rendering
	if pyvista.OFF_SCREEN:
		subplotter.screenshot("fields.png", transparent_background=transparent, window_size=[figsize_x, figsize_y])
	else:
		subplotter.show()

	subplotter.screenshot("fields.png", transparent_background=transparent, window_size=[figsize_x, figsize_y])

	subplotter.app.exec_()

Here is a comparison of the field on axis with GetDP and elmerfem. GetDP and elmerfem use the exact same mesh.

resultBOnAxis
resultBOnAxis825×600 38.9 KB

Best,

Frederic

3 Likes
2

I just want to say a big thank you, Frédéric. I had tried to compute the magnetic field in 3D of a permanent magnet, without getting a sensical answer (here was my attempt: Magnetic field in 3D of a "permanent magnet"). When I get the time, I’ll try to use your code as a guide, to fix my code.
Have a good day.

3

Hi raboynics,

You are welcome.

I am no expert in FEM, I would say I am advanced end-user, but I have some interest in using opensource software to deal with multphysics modelling. There is a clear need for non-commercial options.

Some colleagues of mine and I are currently exploring different opensource software to model coils. The idea is to develop some multiphysics FE tools providing alternative options end-users, in our specific case, for the designers of superconducting magnets.

We are currently benchmarking against a couple of commercial software with a colleague, expert in the modelling of superconducting magnet. The code can be used with some level of confidence. Nonetheless, there are still some details to fix. In particular, we found out that the “air” radius should be increased to get a better results on the field.

I am planning to do a nonlinear case in the near future. I will post an updated version when ready.

Best,

Frederic

1 Like
4

Hello, friend!
I tried to run your code, but there was a import error regarding the VectorFunctionSpace.
Do you what can be used to replace this?
Thank you!

5

Please search the forum, using VectorFunctionSpace as keyword and sorting the results by date. There was a backward incompatible change in dolfinx 0.8.0 that removed that class, but the question has been asked several times in the recent past, and we won’t be answering it again.

1 Like