I am trying to run a magnetostatic simulation on a simple cylinder domain where both ends are kept at a constant charge. I can successfully solve the PDE to calculate the potential field and from that get the electric vector field but I’m having trouble with adding the magnetic field.
I calculate the current density directly from the electric field resulting in a vector field, the J field.
J = \sigma E
where \sigma is the conductivity. Then my intention is to use the Poisson PDE to solve for the magnetic vector potential A.
\nabla^2A = -\mu J
where \mu is the permeability. My intention was to use the following to solve for A but this fails with a UFL error: Can only integrate scalar expressions. The integrand is a tensor expression with value shape (3,) and free indices with labels ().
# Construct the PDE and solve for the magnetic field
# Set up the functions
V = FunctionSpace(domain, ("CG", 1))
u = TrialFunction(V)
v = TestFunction(V)
a = dot(nabla_grad(u), nabla_grad(v)) * dx
L = mu * J_field * v * dx
# Solve the linear problem
ah = Function(V)
problem = LinearProblem(a, L, u=ah, bcs=bcs)
print("Now solving the magnetic field PDE...")
problem.solve()
The full code example is given below.
from mpi4py import MPI
from dolfinx.io import gmshio
from dolfinx.io.utils import XDMFFile
import numpy as np
import gmsh
from ufl import TestFunction, TrialFunction, dot, dx, nabla_grad, curl
from dolfinx.fem import (Constant, Function, FunctionSpace, Expression, VectorFunctionSpace,
dirichletbc, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from petsc4py.PETSc import ScalarType
import matplotlib.pyplot as plt
import pyvista
from dolfinx.plot import create_vtk_mesh
# specified in mm
Tube_hight = 250
Tube_radius = 50
tube_perm_abs = 8.85e-12 # 4 * np.pi*1e-7 # Vacuum
tube_perm = 1 # tube_perm_abs*70 # 4 * np.pi*1e-7 # Vacuum
tube_cond = 1 # Conductivity of medium inside the tube (milk)
tube_permeability = 1 # Permiability of medium inside the tube (milk)
obj1_perm = 100 # Object with large permittivity # tube_perm_abs*140
obj1_cond = 1e-4 # small conductivity
obj1_permeability = 0.01 # small permiability
obj2_perm = 0.01 # Object with small permittivity # tube_perm_abs*10
obj2_cond = 1e4 # large conductivity
obj2_permeability = 100 # large permiability
gmsh.initialize()
max_mesh = Tube_radius * 0.1
min_mesh = Tube_radius * 0.05
min_dist = Tube_radius * 0.1
max_dist = Tube_radius * 1.5
# Create model
main_model = gmsh.model()
main_model.add('main_mesh')
main_model.setCurrent('main_mesh')
# Tags:
tube_interior_tag = 10 # 10 = tube interior
tude_sides_tag = 101 # 101 = tube sides
top_disk_tag = 102 # 102 = top disk
bottom_disk_tag = 100 # 100 = bottom disk
obj1_tag = 103 # 103 = obj1 (large perm)
obj2_tag = 104 # 104 = obj2 (small perm)
cylinder = main_model.occ.addCylinder(0, 0, 0, 0, 0, Tube_hight, Tube_radius)
# Add an object with large permitivity
obj1 = main_model.occ.addSphere(0, 0, Tube_hight // 2 - 50, 5)
# Add an object with large permitivity
obj2 = main_model.occ.addSphere(0, 0, Tube_hight // 2 + 50, 5)
# Use fragment to make sure that all interfaces are conformal.
cell_fragments = main_model.occ.fragment([(3, cylinder)], [(3, obj1), (3, obj2)])
# Tag the physical groups
main_model.occ.synchronize()
important_areas = []
important_surfaces = []
# Find the volumes based on COM
for vol in main_model.getEntities(dim=3):
com = main_model.occ.getCenterOfMass(vol[0], vol[1])
if np.allclose(com, [0, 0, Tube_hight // 2 - 50]):
main_model.addPhysicalGroup(3, [vol[1]], tag=obj1_tag, name="obj1")
important_areas.append(vol)
elif np.allclose(com, [0, 0, Tube_hight // 2 + 50]):
main_model.addPhysicalGroup(3, [vol[1]], tag=obj2_tag, name="obj2")
important_areas.append(vol)
else:
main_model.addPhysicalGroup(3, [vol[1]], tag=tube_interior_tag, name="volume")
# Find the facets based on COM
for vol in main_model.getEntities(dim=2):
com = main_model.occ.getCenterOfMass(vol[0], vol[1])
if np.allclose(com, [0, 0, 0]):
main_model.addPhysicalGroup(2, [vol[1]], tag=bottom_disk_tag, name="bottom_disk")
important_surfaces.append(vol)
elif np.allclose(com, [0, 0, Tube_hight]):
main_model.addPhysicalGroup(2, [vol[1]], tag=top_disk_tag, name="top_disk")
important_surfaces.append(vol)
elif np.allclose(com, [0, 0, Tube_hight / 2]):
main_model.addPhysicalGroup(2, [vol[1]], tag=tude_sides_tag, name="side_surface")
# Set adaptive mesh resolutions
surfaces = main_model.getBoundary(important_areas, oriented=False)
curves = main_model.getBoundary(important_surfaces, oriented=False)
main_model.mesh.field.add("Distance", 1)
main_model.mesh.field.setNumber(1, "Sampling", 200)
# main_model.mesh.field.setNumbers(1, "CurvesList", [e[1] for e in curves])
main_model.mesh.field.setNumbers(1, "SurfacesList", [e[1] for e in surfaces])
main_model.mesh.field.add("Threshold", 2)
main_model.mesh.field.setNumber(2, "InField", 1)
main_model.mesh.field.setNumber(2, "LcMin", min_mesh)
main_model.mesh.field.setNumber(2, "LcMax", max_mesh)
main_model.mesh.field.setNumber(2, "DistMin", min_dist)
main_model.mesh.field.setNumber(2, "DistMax", max_dist)
# The background mesh size is the minimum of all mesh fields
main_model.mesh.field.add("Min", 5)
main_model.mesh.field.setNumbers(5, "FieldsList", [2])
main_model.mesh.field.setAsBackgroundMesh(5)
main_model.occ.synchronize()
main_model.mesh.generate(3)
main_model.mesh.optimize("Netgen")
gmsh.write("mesh3D.msh")
gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
domain, ct, ft = gmshio.model_to_mesh(main_model, mesh_comm, gmsh_model_rank, gdim=3)
gmsh.finalize()
# Set up the functions
V = FunctionSpace(domain, ("CG", 1))
u = TrialFunction(V)
v = TestFunction(V)
# Define the material parameters
tdim = domain.topology.dim
fdim = tdim - 1
Q_eps = FunctionSpace(domain, ("DG", 0))
Q_mu = FunctionSpace(domain, ("DG", 0))
Q_cond = FunctionSpace(domain, ("DG", 0))
material_tags = np.unique(ct.values)
eps = Function(Q_eps)
mu = Function(Q_mu)
cond = Function(Q_cond)
# As we only set some values in eps, initialize all as vacuum
eps.x.array[:] = tube_perm
cond.x.array[:] = tube_cond
mu.x.array[:] = tube_permeability
for tag in material_tags:
cells = ct.find(tag)
if tag == obj1_tag:
eps_ = obj1_perm
cond_ = obj1_cond
mu_ = obj1_permeability
elif tag == obj2_tag:
eps_ = obj2_perm
cond_ = obj2_cond
mu_ = obj2_permeability
else:
eps_ = tube_perm
cond_ = tube_cond
mu_ = tube_permeability
eps.x.array[cells] = np.full_like(cells, eps_, dtype=ScalarType)
cond.x.array[cells] = np.full_like(cells, cond_, dtype=ScalarType)
mu.x.array[cells] = np.full_like(cells, mu_, dtype=ScalarType)
domain.topology.create_connectivity(tdim, fdim)
# Set the boundary conditions
# The top boundary is set to 10V
top_disk_facets = ft.indices[ft.values == top_disk_tag]
top_disk_dofs = locate_dofs_topological(V, fdim, top_disk_facets)
u_T = Constant(domain, ScalarType(10))
top_bc = dirichletbc(u_T, top_disk_dofs, V)
# The bottom boundary is set to -10V
bottom_disk_facets = ft.indices[ft.values == bottom_disk_tag]
bottom_disk_dofs = locate_dofs_topological(V, fdim, bottom_disk_facets)
u_B = Constant(domain, ScalarType(-10))
bottom_bc = dirichletbc(u_B, bottom_disk_dofs, V)
bcs = [top_bc, bottom_bc] #, sides_bc]
# Construct the PDE and solve for the electric potential
a = dot(nabla_grad(u), nabla_grad(v)) * eps * dx
L = Constant(domain, ScalarType(1e-16)) * v * dx
# Solve the linear problem
uh = Function(V)
problem = LinearProblem(a, L, u=uh, bcs=bcs)
print("Now solving the electric potential PDE...")
problem.solve()
# calculate the electric field E and the current density J
E_field = -nabla_grad(uh)
J_field = cond * E_field
# project to mesh
V_E_field = VectorFunctionSpace(domain, ("DG", 0))
V_J_field = VectorFunctionSpace(domain, ("DG", 0))
E_field_expr = Expression(E_field, V_E_field.element.interpolation_points())
J_field_expr = Expression(J_field, V_J_field.element.interpolation_points())
E_field_projected = Function(V_E_field)
J_field_projected = Function(V_J_field)
E_field_projected.interpolate(E_field_expr)
J_field_projected.interpolate(J_field_expr)
# Construct the PDE and solve for the magnetic field
# Set up the functions
V = FunctionSpace(domain, ("CG", 1))
u = TrialFunction(V)
v = TestFunction(V)
a = dot(nabla_grad(u), nabla_grad(v)) * dx
L = mu * J_field * v * dx
# Solve the linear problem
ah = Function(V)
problem = LinearProblem(a, L, u=ah, bcs=bcs)
print("Now solving the magnetic field PDE...")
problem.solve()
b_field = curl(ah)
V_B_field = VectorFunctionSpace(domain, ("DG", 0))
B_field_expr = Expression(b_field, V_B_field.element.interpolation_points())
B_field_projected = Function(V_B_field)
B_field_projected.interpolate(B_field_expr)
What is the correct way to deal with vector values in the integrand? Should I try to solve each component separately?