Dear @dokken and all,
I’m working with DOLFINx 0.9 to simulate the time evolution of eddy currents in a 2D axisymmetric electromagnetic model. The model is driven by a decaying toroidal current, and the goal is to evaluate how currents are induced in passive conductors (e.g., copper coils).
I have successfully implemented the transient formulation using:
CG1for the poloidal magnetic flux \Psi(r, z)DG0for material properties \mu, \sigma
The governing equation is:
\nabla \cdot \left( \frac{1}{2\pi r} \mu^{-1} \nabla \Psi \right) + \frac{\sigma}{2\pi r} \frac{\partial \Psi}{\partial t} = J_\phi
In the attached example:
- The source current J_\phi is applied over a small circular region (
wire) - Two passive rectangular copper coils (
passive1,passive2) experience induced currents via the \sigma \frac{\partial \Psi}{\partial t} term - I use
Expressionandassemble_scalarto compute the total induced current in each region at every time step
Goal
I would now like to upgrade the model by explicitly introducing electrical constraints among passive components, modeled as 0D circuits. Specifically, I want to couple each group of passive conductors (defined by physical groups) to a 0D Faraday’s law:
L \frac{di}{dt} + R i = V_0 - \frac{d \bar{\Psi}}{dt}
Where:
- i(t): total current in the passive circuit
- \bar{\Psi}: average poloidal flux over the materials of the circuit
- L, R: known inductance and resistance values
- V_0: external driving voltage (assumed 0 for now)
In the current example, I would like to connect the two passive coils in series so as to have the same Amps of eddy currents on both conductors.
Objectives
- Define one or more circuits, each composed of a subset of physical groups
- Couple the 2D FEM problem with ODEs for these circuits
- Ensure the flux consistency between the average magnetic flux and the circuit equation
- Obtain a numerically stable implementation, suitable for use inside nonlinear iterations (e.g., Picard iterations)
Questions
- How would you suggest implementing this FEM–circuit coupling in DOLFINx?
- Is it feasible to represent the 0D current(s) as additional unknown(s) in the same linear system? (e.g. using a mixed formulation or block system)
- Are there any recommended strategies to ensure numerical stability and scalability of this approach?
Any insights, references, or example structures would be extremely appreciated.
Thanks in advance for your help!
import numpy as np
import gmsh
from mpi4py import MPI
from dolfinx import default_scalar_type
from dolfinx.io.gmshio import model_to_mesh
from dolfinx.fem import (
dirichletbc,
Expression,
Function,
functionspace,
locate_dofs_topological,
form,
assemble_scalar,
)
from dolfinx.fem.petsc import (
apply_lifting,
assemble_matrix,
assemble_vector,
create_vector,
set_bc,
)
from dolfinx.mesh import locate_entities_boundary
from ufl import (
TestFunction,
TrialFunction,
dot,
dx,
grad,
SpatialCoordinate,
Measure,
)
from petsc4py import PETSc
import matplotlib.pyplot as plt
model_rank = 0
mesh_comm = MPI.COMM_WORLD
# Geometry parameters
R = 10
r_wire = 0.005
rw, zw = 0.5, 0.0
r_cu1_offset = -0.1
z_cu1_offset = 0.5
r_cu2_offset = -0.05
z_cu2_offset = 0.25
h_cu = 0.2
w_cu = 0.2
lc_air = 1
lc_wire = 0.001
lc_cu = 0.25
# Tags
tag_air = 1
tag_wire = 2
tag_passive1 = 3
tag_passive2 = 4
rho_cu = 1.68e-8
sigma_cu = 1 / rho_cu
Iw0 = 1e6
deltat = 0.05
dt = 1e-3
def Iw(t):
return Iw0 * (1 - t / deltat)
def integrate_f(f, mesh, ct=None, tag=None) -> float:
dx_obj = dx if ct is None else Measure("dx", subdomain_data=ct, domain=mesh)
return (
assemble_scalar(form(f * dx_obj))
if tag is None
else assemble_scalar(form(f * dx_obj(tag)))
)
Aw = r_wire * np.pi**2
gmsh.initialize()
gmsh.model.add("Transient with circuits")
# --- Air region (half disk + rectangle to close D shape) ---
p_center = gmsh.model.occ.addPoint(0.0, 0.0, 0, lc_air)
p_top = gmsh.model.occ.addPoint(0.0, R, 0, lc_air)
p_right = gmsh.model.occ.addPoint(R, 0.0, 0, lc_air)
p_bottom = gmsh.model.occ.addPoint(0.0, -R, 0, lc_air)
arc1 = gmsh.model.occ.addCircleArc(p_top, p_center, p_right)
arc2 = gmsh.model.occ.addCircleArc(p_right, p_center, p_bottom)
line1 = gmsh.model.occ.addLine(p_bottom, p_top)
loop_air = gmsh.model.occ.addCurveLoop([arc1, arc2, line1])
gmsh.model.occ.synchronize()
# --- Circular wire (full disk) ---
wire = gmsh.model.occ.addDisk(rw, zw, 0.0, r_wire, r_wire, tag=tag_wire)
gmsh.model.occ.synchronize()
# --- Passive coil 1 (rectangle) ---
z1_bot = zw + z_cu1_offset - h_cu / 2
r1_left = rw + r_cu1_offset - w_cu / 2
passive1 = gmsh.model.occ.addRectangle(
r1_left, z1_bot, 0.0, w_cu, h_cu, tag=tag_passive1
)
gmsh.model.occ.synchronize()
# --- Passive coil 2 (rectangle) ---
z2_bot = zw - z_cu2_offset - h_cu / 2
r2_left = rw + r_cu2_offset - w_cu / 2
passive2 = gmsh.model.occ.addRectangle(
r2_left, z2_bot, 0.0, w_cu, h_cu, tag=tag_passive2
)
air = gmsh.model.occ.addPlaneSurface([loop_air, wire, passive1, passive2], tag=tag_air)
gmsh.model.occ.synchronize()
# --- Physical surfaces ---
gmsh.model.addPhysicalGroup(2, [air], tag_air)
gmsh.model.setPhysicalName(2, tag_air, "Air")
gmsh.model.addPhysicalGroup(2, [wire], tag_wire)
gmsh.model.setPhysicalName(2, tag_wire, "Wire")
gmsh.model.addPhysicalGroup(2, [passive1], tag_passive1)
gmsh.model.setPhysicalName(2, tag_passive1, "Passive1")
gmsh.model.addPhysicalGroup(2, [passive2], tag_passive2)
gmsh.model.setPhysicalName(2, tag_passive2, "Passive2")
# --- Mesh options ---
# Mesh size control options
# Gather all refined surfaces
refined_surfaces = [(2, wire), (2, passive1), (2, passive2)]
# Get boundary edges of the refined surfaces
gmsh.model.occ.synchronize()
edges = gmsh.model.getBoundary(refined_surfaces, oriented=False)
# Field 1: distance from these edges
gmsh.model.mesh.field.add("Distance", 1)
gmsh.model.mesh.field.setNumbers(1, "EdgesList", [e[1] for e in edges])
# Field 2: threshold-based size depending on distance
gmsh.model.mesh.field.add("Threshold", 2)
gmsh.model.mesh.field.setNumber(2, "IField", 1)
gmsh.model.mesh.field.setNumber(2, "LcMin", 0.1) # finest near coils
gmsh.model.mesh.field.setNumber(2, "LcMax", 1.0) # coarser away
gmsh.model.mesh.field.setNumber(2, "DistMin", 0.01) # radius within which it's LcMin
gmsh.model.mesh.field.setNumber(2, "DistMax", 0.2) # after which it's LcMax
# Combine (you can add more fields if needed)
gmsh.model.mesh.field.add("Min", 5)
gmsh.model.mesh.field.setNumbers(5, "FieldsList", [2])
gmsh.model.mesh.field.setAsBackgroundMesh(5)
gmsh.option.setNumber(
"Mesh.CharacteristicLengthFactor", 0.25
) # Global scaling factor for all mesh sizes (refines everything)
gmsh.option.setNumber(
"Mesh.MeshSizeFromPoints", 1
) # Use lc values provided at each Point definition
gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 2) # Smoo
gmsh.model.mesh.generate(2)
gmsh.fltk.run() # uncomment to preview mesh
# --- Convert to DolfinX
mesh, ct, _ = model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=2)
gmsh.finalize()
# %%
Q = functionspace(mesh, ("DG", 0))
mu = Function(Q)
J = Function(Q)
sigma = Function(Q)
J.x.array[:] = 0
mu.x.array[:] = 4 * np.pi * 1e-7
passive1_cells = ct.find(tag_passive1)
passive2_cells = ct.find(tag_passive2)
sigma.x.array[passive1_cells] = sigma_cu
sigma.x.array[passive2_cells] = sigma_cu
wire_cells = ct.find(tag_wire)
V = functionspace(mesh, ("Lagrange", 1))
tdim = mesh.topology.dim
facets = locate_entities_boundary(mesh, tdim - 1, lambda x: np.full(x.shape[1], True))
dofs = locate_dofs_topological(V, tdim - 1, facets)
bc = dirichletbc(default_scalar_type(0), dofs, V)
Psi = Function(V)
Psi_old = Function(V)
u = TrialFunction(V)
v = TestFunction(V)
x = SpatialCoordinate(mesh)
a = (1 / mu) * (1 / (2 * np.pi * x[0])) * dot(grad(u), grad(v)) * dx + (
1 / (2 * np.pi * x[0])
) * (sigma / dt) * u * v * dx
L = J * v * dx + (1 / (2 * np.pi * x[0])) * (sigma / dt) * Psi_old * v * dx
print("Assembling A matrix....")
a_form = form(a)
A = assemble_matrix(a_form, bcs=[bc])
A.assemble()
print("A matrix assembled")
print("Creating solver....")
solver = PETSc.KSP().create(mesh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.CG)
solver.getPC().setType(PETSc.PC.Type.SOR)
print("Solver created")
L_form = form(L)
Je_expr = Expression(
(-sigma / (2 * np.pi * x[0]) * (Psi - Psi_old) / dt),
Q.element.interpolation_points(),
)
Je = Function(Q)
t = np.arange(0.0, deltat, dt)
Ie = np.zeros(shape=(len(t), 2), dtype=default_scalar_type)
for i, ti in enumerate(t):
print(f"****** t = {ti:.4f} s ******")
J.x.array[wire_cells] = np.full_like(
wire_cells, Iw(ti) / Aw, dtype=default_scalar_type
)
print("Assembling b vector....")
b = create_vector(L_form)
assemble_vector(b, L_form)
apply_lifting(b, [a_form], [[bc]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, [bc])
print("b vector assembled")
print("Solving...")
solver.solve(b, Psi.x.petsc_vec)
Psi.x.petsc_vec.ghostUpdate(
addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD
)
print("Done")
Je.interpolate(Je_expr)
for j, tag in enumerate([passive1, passive2]):
Ie[i, j] = integrate_f(Je, mesh, ct, tag)
Psi_old.x.petsc_vec.copy(Psi.x.petsc_vec)
# %%
fig, ax1 = plt.subplots()
color1 = "tab:blue"
ax1.set_xlabel("Time [s]")
ax1.set_ylabel("Driving current $I_w(t)$ [A]", color=color1)
line1, = ax1.plot(t, Iw(t), color=color1, label="$I_w(t)$")
ax1.tick_params(axis="y", labelcolor=color1)
ax2 = ax1.twinx()
color2 = "tab:red"
line2, = ax2.plot(t, Ie[:, 0], color="tab:red", linestyle="--", label="$I_{e1}(t)$")
line3, = ax2.plot(t, Ie[:, 1], color="tab:green", linestyle=":", label="$I_{e2}(t)$")
ax2.set_ylabel("Induced currents $I_e(t)$ [A]", color=color2)
ax2.tick_params(axis="y", labelcolor=color2)
lines = [line1, line2, line3]
labels = [line.get_label() for line in lines]
ax1.legend(lines, labels, loc="upper right")
plt.title("Driving and Induced Currents Over Time")
plt.grid(True)
plt.tight_layout()
plt.show()
,
Questions
Example