I am currently working on a FEM simulation in FEnicCSx09 to calculate the inductance of a PCB Coil. My current code and debug is attached below. The results for my coil shape are multiple orders of magnitude too high and I don`t know why. Can someone please help me find the error?
# solve_spiral_inductance_3d.py (verbose)
import argparse, json, time
import numpy as np
import ufl
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import fem, io
from dolfinx.io import gmshio
# ----------------------------
# Hilfsfunktionen
# ----------------------------
def build_rect_centerline(params):
"""Polyline der rechteckigen Spirale (2D, z wird separat gesetzt)."""
N = int(params["num_turns"])
din = float(params["inner_diameter"])
w = float(params["trace_width"])
s = float(params["spacing"])
directions = [(1,0), (0,1), (-1,0), (0,-1)]
seg_len = din # Start = inner_diameter
x = y = 0.0
pts = [(x,y)]
dir_idx = 0
for i in range(N * 4):
dx, dy = directions[dir_idx]
x += dx * seg_len
y += dy * seg_len
pts.append((x,y))
dir_idx = (dir_idx + 1) % 4
if i % 2 == 1:
seg_len += (w + s)
return pts
def cell_centroid(domain, cell):
"""Schwerpunkt einer 3D-Zelle (arithm. Mittel der Eckpunkte)."""
topo = domain.topology
tdim = topo.dim
if topo.connectivity(tdim, 0) is None:
topo.create_connectivity(tdim, 0)
c2v = topo.connectivity(tdim, 0)
vtx = c2v.links(cell)
X = domain.geometry.x[vtx]
return X.mean(axis=0) # (x,y,z)
def segments3d_from_path(centerline_xy, z_spiral_mid, p0, p1, z_top, via_h, z_bridge_mid):
"""
Erzeuge 3D-Segmente mit Typ-Labels:
- Spiral: rückwärts (p1->p0) in der Ebene z_spiral_mid
- Via0: von p0 z_top -> z_top+via_h
- Bridge: p0->p1 auf z_bridge_mid
- Via1: von p1 z_top+via_h -> z_top
Rückgabe: Liste [(A,B,t_hat,L,kind), ...] und Längen-Stats.
"""
segs = []
L_spiral = 0.0
# Spirale rückwärts
for i in range(len(centerline_xy)-1, 0, -1):
A = np.array([centerline_xy[i][0], centerline_xy[i][1], z_spiral_mid])
B = np.array([centerline_xy[i-1][0], centerline_xy[i-1][1], z_spiral_mid])
v = B - A; L = float(np.linalg.norm(v))
if L > 0:
segs.append((A, B, v/L, L, "spiral"))
L_spiral += L
# Via an p0: nach oben
A = np.array([p0[0], p0[1], z_top]); B = np.array([p0[0], p0[1], z_top + via_h])
v = B - A; L_via0 = float(np.linalg.norm(v))
segs.append((A, B, v/L_via0, L_via0, "via0"))
# Brücke: p0 -> p1, mittig in der Brücke
A = np.array([p0[0], p0[1], z_bridge_mid]); B = np.array([p1[0], p1[1], z_bridge_mid])
v = B - A; L_bridge = float(np.linalg.norm(v))
segs.append((A, B, v/L_bridge, L_bridge, "bridge"))
# Via an p1: nach unten
A = np.array([p1[0], p1[1], z_top + via_h]); B = np.array([p1[0], p1[1], z_top])
v = B - A; L_via1 = float(np.linalg.norm(v))
segs.append((A, B, v/L_via1, L_via1, "via1"))
stats = {
"L_spiral": L_spiral,
"L_via0": L_via0,
"L_bridge": L_bridge,
"L_via1": L_via1,
"L_total": L_spiral + L_via0 + L_bridge + L_via1
}
return segs, stats
def nearest_tangent_3d(segs3d, pxyz):
"""Tangente des nächstgelegenen Segmentes an Punkt pxyz."""
p = np.array(pxyz)
best_d2 = 1e300
best_t = np.array([1.0, 0.0, 0.0])
for A, B, t_hat, L, kind in segs3d:
v = B - A
t = np.dot(p - A, v) / (L*L)
t = min(1.0, max(0.0, t))
q = A + t*v
d2 = float(np.sum((p - q)**2))
if d2 < best_d2:
best_d2 = d2
best_t = t_hat
return best_t
# ----------------------------
# Main
# ----------------------------
def main():
# --- imports lokal halten, damit Drop-in einfach ist ---
import argparse, json, time
import numpy as np
import ufl
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import fem, io
from dolfinx.io import gmshio
parser = argparse.ArgumentParser()
parser.add_argument("--msh", default="rect_spiral_with_inner_d.msh",
help="Gmsh mesh (MSH4.1) (default: rect_spiral_with_inner_d.msh)")
parser.add_argument("--meta", default="spiral_meta.json",
help="JSON metadata file (default: spiral_meta.json)")
parser.add_argument("--I", type=float, default=1.0, help="Total current [A]")
parser.add_argument("--alpha", type=float, default=1e-4, help="Gauge parameter (div-div penalty)")
parser.add_argument("--mu_r_air", type=float, default=1.0)
parser.add_argument("--mu_r_cu", type=float, default=1.0)
parser.add_argument("--save-xdmf", type=int, default=1, help="Write XDMF with A and |B|")
parser.add_argument("--xdmf", default="solution_fields.xdmf")
parser.add_argument("--debug", type=int, default=1, help="Verbose prints (0/1)")
args = parser.parse_args()
comm = MPI.COMM_WORLD
rank = comm.rank
dbg = (args.debug == 1)
t0 = time.perf_counter()
if rank == 0 and dbg:
print(f"\n[dbg] start solve_spiral_inductance_3d: I={args.I}, alpha={args.alpha}\n")
# --- Mesh laden ---
t = time.perf_counter()
domain, cell_tags, facet_tags = gmshio.read_from_msh(args.msh, comm, rank)
if rank == 0 and dbg:
print(f"[dbg] read_from_msh: {time.perf_counter()-t:.3f} s")
tdim = domain.topology.dim
COIL = 1
AIR = 2
coil_cells = cell_tags.find(COIL)
air_cells = cell_tags.find(AIR)
if rank == 0:
total_cells = cell_tags.indices.size
print(f"[info] cells: total={total_cells}, coil={coil_cells.size}, air={air_cells.size}")
# --- Metadaten lesen ---
if rank == 0:
with open(args.meta, "r") as f:
meta = json.load(f)
else:
meta = None
meta = comm.bcast(meta, root=0)
params = meta.get("params", {})
num_layers = int(params.get("num_layers", 1))
layer_offset = float(params.get("layer_offset", 0.0))
thickness = float(params.get("thickness", 0.0))
w_trace = float(params.get("trace_width", 0.0))
N = int(params.get("num_turns", 0))
din = float(params.get("inner_diameter", 0.0))
s = float(params.get("spacing", 0.0))
if rank == 0:
print(f"[info] params: N={N}, din={din}, w={w_trace}, s={s}, layers={num_layers}, t={thickness}, offset={layer_offset}")
# Return-Loop Info (falls vorhanden, sonst Fallbacks)
rmeta = meta.get("return_loop", {})
via_h = float(rmeta.get("via_h", max(6.0*thickness, 0.2*w_trace)))
bridge_tz = float(rmeta.get("bridge_tz", max(2.0*thickness, 0.02*w_trace)))
z_top = float(rmeta.get("z_top", thickness))
z_bridge = float(rmeta.get("z_bridge", z_top + max(6.0*thickness, 0.2*w_trace)))
z_spiral_mid = thickness * 0.5
z_bridge_mid = z_bridge + 0.5 * bridge_tz
if rank == 0 and dbg:
print(f"[dbg] return_loop: via_h={via_h:g}, bridge_tz={bridge_tz:g}, z_top={z_top:g}, z_bridge={z_bridge:g}")
print(f"[dbg] z_spiral_mid={z_spiral_mid:g}, z_bridge_mid={z_bridge_mid:g}")
# --- Konstanten & Räume ---
mu0 = 4e-7 * np.pi
mu_air = mu0 * args.mu_r_air
mu_cu = mu0 * args.mu_r_cu
if rank == 0 and dbg:
print(f"[dbg] mu_air={mu_air:.6e}, mu_cu={mu_cu:.6e}")
V = fem.FunctionSpace(domain, ("N1curl", 1)) # H(curl)
Q = fem.FunctionSpace(domain, ("DG", 0)) # skalar DG0 (für invmu)
A = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
dx = ufl.Measure("dx", domain=domain, subdomain_data=cell_tags)
# --- 1/mu (DG0) ---
t = time.perf_counter()
invmu = fem.Function(Q)
invmu.x.array[:] = 0.0
dofs_air = fem.locate_dofs_topological(Q, tdim, air_cells)
dofs_coil = fem.locate_dofs_topological(Q, tdim, coil_cells)
invmu.x.array[dofs_air] = 1.0/mu_air
invmu.x.array[dofs_coil] = 1.0/mu_cu
invmu.x.scatter_forward()
if rank == 0 and dbg:
print(f"[dbg] invmu set (air/coils): {time.perf_counter()-t:.3f} s, "
f"n_dofs(Q)={Q.dofmap.index_map.size_local}")
# --- Stromdichte-Skalierung & Geometriepfad ---
I = float(args.I)
if w_trace <= 0 or thickness <= 0:
raise RuntimeError("trace_width or thickness missing/non-positive in meta.")
# 2D-Mittellinie & Ports (Fallbacks aus meta)
cl = build_rect_centerline({
"num_turns": N,
"inner_diameter": din,
"trace_width": w_trace,
"spacing": s
})
p0 = tuple(meta.get("port_plus_xy", [cl[0][0], cl[0][1]]))
p1 = tuple(meta.get("port_minus_xy", [cl[-1][0], cl[-1][1]]))
# 3D-Segmente & Stats
t = time.perf_counter()
segs3d, segstats = segments3d_from_path(cl, z_spiral_mid, p0, p1, z_top, via_h, z_bridge_mid)
if rank == 0 and dbg:
print(f"[dbg] built segments3d in {time.perf_counter()-t:.3f} s (n={len(segs3d)})")
print(f" L_spiral={segstats['L_spiral']:.6e}, "
f"L_via0={segstats['L_via0']:.6e}, "
f"L_bridge={segstats['L_bridge']:.6e}, "
f"L_via1={segstats['L_via1']:.6e}, "
f"L_total={segstats['L_total']:.6e}")
# Volumen der Spule (aus Mesh) und Pfadlänge (3D)
one = fem.Function(Q); one.x.array[:] = 1.0
V_coil_loc = fem.assemble_scalar(fem.form(one * dx(COIL)))
V_coil = domain.comm.allreduce(V_coil_loc, op=MPI.SUM)
L_path = segstats["L_total"]
# 2D-Pfadlänge + V_theo (nur Diagnose)
L_theo2d = 0.0
for i in range(1, len(cl)):
dx2 = cl[i][0] - cl[i-1][0]
dy2 = cl[i][1] - cl[i-1][1]
L_theo2d += float((dx2*dx2 + dy2*dy2)**0.5)
V_theo = L_theo2d * w_trace * thickness
if rank == 0:
print(f"[diag] V_coil(mesh) = {V_coil:.3e} m^3")
print(f"[diag] L_path(3D) = {L_path:.3e} m")
print(f"[diag] L_theo2d = {L_theo2d:.3e} m, V_theo ≈ {V_theo:.3e} m^3")
if V_coil <= 0 or L_path <= 0:
raise RuntimeError("V_coil oder L_path ist 0 – Mesh/Tags prüfen.")
# --- Outer Airbox Dirichlet-BC: A × n = 0 (magnetic insulation) ---
from dolfinx import mesh as dmesh
# Airbox aus Meta: [ax, ay, az, adx, ady, adz]
ax, ay, az, adx, ady, adz = [float(v) for v in meta.get("airbox", [0, 0, 0, 0, 0, 0])]
if (adx <= 0) or (ady <= 0) or (adz <= 0):
raise RuntimeError("Airbox-Daten fehlen oder sind ungültig in meta['airbox'].")
# Geometrische Toleranz für die Ebenen (1e-6 der Boxgröße ist robust)
tol = 1e-6 * max(adx, ady, adz)
# Marker: alle **äußeren** Facetten der Box (x==ax|ax+adx oder y==... oder z==...)
def on_outer_boundary(x):
near_xmin = np.isclose(x[0], ax, atol=tol)
near_xmax = np.isclose(x[0], ax + adx, atol=tol)
near_ymin = np.isclose(x[1], ay, atol=tol)
near_ymax = np.isclose(x[1], ay + ady, atol=tol)
near_zmin = np.isclose(x[2], az, atol=tol)
near_zmax = np.isclose(x[2], az + adz, atol=tol)
return near_xmin | near_xmax | near_ymin | near_ymax | near_zmin | near_zmax
facet_dim = tdim - 1
outer_facets = dmesh.locate_entities_boundary(domain, facet_dim, on_outer_boundary)
# DOFs des H(curl)-Raums auf diesen Facetten
outer_dofs = fem.locate_dofs_topological(V, facet_dim, outer_facets)
# Zero-Function als Wert (setzt die **tangentialen** Nédélec-DOFs auf 0)
A_zero = fem.Function(V)
A_zero.x.array[:] = 0.0
bc_outer = fem.dirichletbc(A_zero, outer_dofs)
if rank == 0 and dbg:
print(f"[dbg] outer BC: facets={outer_facets.size}, dofs={outer_dofs.size}")
# robuste J-Skalierung
Jmag = I * (L_path / V_coil)
J0_th = I / (w_trace * thickness)
if rank == 0:
print(f"[diag] Jmag = {Jmag:.3e} A/m^2 (skalierte J-Größe)")
print(f"[diag] J0_th (I/(w*t)) = {J0_th:.3e} A/m^2, ratio Jmag/J0_th = {Jmag/J0_th:.3f}")
# --- J als DG0-Vektorfeld (zellweise konstant entlang 3D-Tangente) ---
t = time.perf_counter()
Wv = fem.VectorFunctionSpace(domain, ("DG", 0))
J = fem.Function(Wv); J.x.array[:] = 0.0
jx_vals, jy_vals, jz_vals = [], [], []
for cell in coil_cells:
cx, cy, cz = cell_centroid(domain, int(cell))
t_hat = nearest_tangent_3d(segs3d, (cx, cy, cz))
jx_vals.append(Jmag * t_hat[0])
jy_vals.append(Jmag * t_hat[1])
jz_vals.append(Jmag * t_hat[2])
Jx = J.sub(0); Vx = Jx.function_space
Jy = J.sub(1); Vy = Jy.function_space
Jz = J.sub(2); Vz = Jz.function_space
dofs_x = fem.locate_dofs_topological(Vx, tdim, coil_cells)
dofs_y = fem.locate_dofs_topological(Vy, tdim, coil_cells)
dofs_z = fem.locate_dofs_topological(Vz, tdim, coil_cells)
Jx.x.array[dofs_x] = np.asarray(jx_vals, dtype=np.float64)
Jy.x.array[dofs_y] = np.asarray(jy_vals, dtype=np.float64)
Jz.x.array[dofs_z] = np.asarray(jz_vals, dtype=np.float64)
J.x.scatter_forward()
if rank == 0 and dbg:
nWv = Wv.dofmap.index_map.size_local * Wv.dofmap.index_map_bs
print(f"[dbg] filled J (DG0 vec) in {time.perf_counter()-t:.3f} s, n_dofs(Wv)={nWv}")
# zusätzliche J-Diagnosen
J_L2_loc = fem.assemble_scalar(fem.form(ufl.inner(J, J) * dx(COIL)))
J_L2 = domain.comm.allreduce(J_L2_loc, op=MPI.SUM)
if rank == 0 and dbg:
print(f"[dbg] ||J||_L2^2 over COIL = {J_L2:.6e}")
print(f"[dbg] Jx min/max on set DOFs = {np.min(Jx.x.array[dofs_x]):.3e} / {np.max(Jx.x.array[dofs_x]):.3e}")
print(f"[dbg] Jy min/max on set DOFs = {np.min(Jy.x.array[dofs_y]):.3e} / {np.max(Jy.x.array[dofs_y]):.3e}")
print(f"[dbg] Jz min/max on set DOFs = {np.min(Jz.x.array[dofs_z]):.3e} / {np.max(Jz.x.array[dofs_z]):.3e}")
# ===== Magnetisierungs-RHS: M(x) = 1/2 * (x - c_cell) x J_cell (DG0) =====
t = time.perf_counter()
Vsc = fem.FunctionSpace(domain, ("DG", 0))
Cx = fem.Function(Vsc); Cy = fem.Function(Vsc); Cz = fem.Function(Vsc)
Cx.x.array[:] = 0.0; Cy.x.array[:] = 0.0; Cz.x.array[:] = 0.0
dofs_sc = fem.locate_dofs_topological(Vsc, tdim, coil_cells) # DG0: 1 DOF pro Zelle
cx_list, cy_list, cz_list = [], [], []
for cell in coil_cells:
cx, cy, cz = cell_centroid(domain, int(cell))
cx_list.append(cx); cy_list.append(cy); cz_list.append(cz)
Cx.x.array[dofs_sc] = np.asarray(cx_list, dtype=np.float64)
Cy.x.array[dofs_sc] = np.asarray(cy_list, dtype=np.float64)
Cz.x.array[dofs_sc] = np.asarray(cz_list, dtype=np.float64)
Cx.x.scatter_forward(); Cy.x.scatter_forward(); Cz.x.scatter_forward()
x = ufl.SpatialCoordinate(domain)
C = ufl.as_vector((Cx, Cy, Cz))
r = ufl.as_vector((x[0] - C[0], x[1] - C[1], x[2] - C[2])) # x - c_cell
M = 0.5 * ufl.cross(r, J)
if rank == 0 and dbg:
print(f"[dbg] built magnetization field M in {time.perf_counter()-t:.3f} s")
# M-Diagnosen
M_L2_loc = fem.assemble_scalar(fem.form(ufl.inner(M, M) * dx(COIL)))
M_L2 = domain.comm.allreduce(M_L2_loc, op=MPI.SUM)
if rank == 0 and dbg:
print(f"[dbg] ||M||_L2^2 over COIL = {M_L2:.6e}")
# RHS-Tests: Norm der RHS-Vektoren für beide Formen
bM_vec = fem.petsc.assemble_vector(fem.form(ufl.inner(M, ufl.curl(v)) * dx(COIL)))
bM_norm2 = bM_vec.norm(PETSc.NormType.NORM_2)
bJ_vec = fem.petsc.assemble_vector(fem.form(ufl.inner(J, v) * dx(COIL)))
bJ_norm2 = bJ_vec.norm(PETSc.NormType.NORM_2)
if rank == 0:
print(f"[diag] ||b_M||_2 = {bM_norm2:.6e} (RHS for M·curl(v))")
print(f"[diag] ||b_J||_2 = {bJ_norm2:.6e} (RHS for J·v)")
use_Mcurlv = True
if not np.isfinite(bM_norm2) or bM_norm2 < 1e-14:
if rank == 0:
print("[warn] b_M is ~0 → switching RHS to ∫ J·v dx(COIL) for this run.")
use_Mcurlv = False
# ----------------------------
# Variationsform & Lösung
# ----------------------------
curlA = ufl.curl(A); curlv = ufl.curl(v)
divA = ufl.div(A); divv = ufl.div(v)
a_form = ufl.inner(invmu * curlA, curlv) * dx + args.alpha * ufl.inner(divA, divv) * dx
L_form = ufl.inner(M, curlv) * dx(COIL) if use_Mcurlv else ufl.inner(J, v) * dx(COIL)
petsc_opts = {
"ksp_type": "cg",
"pc_type": "hypre",
"ksp_rtol": 1e-8,
"ksp_max_it": 50000
}
if dbg:
petsc_opts.update({"ksp_monitor_true_residual": None})
# (optional) Matrix- und RHS-Diagnose
if dbg:
A_mat = fem.petsc.assemble_matrix(fem.form(a_form)); A_mat.assemble()
b_vec = fem.petsc.assemble_vector(fem.form(L_form))
if rank == 0:
m, n = A_mat.getSize()
print(f"[dbg] assembled A: {m}x{n}, ||b||_2 (used) = {b_vec.norm(PETSc.NormType.NORM_2):.6e}")
t = time.perf_counter()
problem = fem.petsc.LinearProblem(a_form, L_form, bcs=[bc_outer], petsc_options=petsc_opts)
A_sol = problem.solve()
if rank == 0 and dbg:
print(f"[dbg] linear solve done in {time.perf_counter()-t:.3f} s")
# ----------------------------
# Energie & Induktivität
# ----------------------------
t = time.perf_counter()
energy_density = 0.5 * ufl.inner(invmu * ufl.curl(A_sol), ufl.curl(A_sol))
Wm_local = fem.assemble_scalar(fem.form(energy_density * dx))
Wm = domain.comm.allreduce(Wm_local, op=MPI.SUM)
L_val = 2.0 * Wm / (I * I)
if rank == 0:
print("---------- Ergebnisse ----------")
print(f"I = {I:.6g} A")
print(f"W_m = {Wm:.6e} J")
print(f"L (Num) = {L_val:.6e} H ({L_val*1e6:.3f} µH)")
print(f"alpha (gauge) = {args.alpha:g}")
print("Hinweis: alpha zu klein -> langsamer; zu groß -> leichte Verfälschung (1e-5..3e-4 üblich).")
if rank == 0 and dbg:
print(f"[dbg] energy & L computed in {time.perf_counter()-t:.3f} s")
print(f"[dbg] total wall time: {time.perf_counter()-t0:.3f} s\n")
# ----------------------------
# Feldausgabe (optional)
# ----------------------------
if args.save_xdmf:
t = time.perf_counter()
with io.XDMFFile(comm, args.xdmf, "w") as xdmf:
xdmf.write_mesh(domain)
A_sol.name = "A"
xdmf.write_function(A_sol)
# |B| ≡ |curl A| als DG0-Projektion
V0 = fem.FunctionSpace(domain, ("DG", 0))
bmag_trial = ufl.TrialFunction(V0)
bmag_test = ufl.TestFunction(V0)
b_abs = ufl.sqrt(ufl.inner(ufl.curl(A_sol), ufl.curl(A_sol)))
Mmat = fem.petsc.assemble_matrix(fem.form(ufl.inner(bmag_trial, bmag_test) * dx))
Mmat.assemble()
rhs = fem.petsc.assemble_vector(fem.form(ufl.inner(b_abs, bmag_test) * dx))
kspM = PETSc.KSP().create(comm)
kspM.setOperators(Mmat)
kspM.setType("preonly")
kspM.getPC().setType("lu")
Bmag = fem.Function(V0)
kspM.solve(rhs, Bmag.vector)
kspM.destroy(); Mmat.destroy()
Bmag.name = "B_mag"
xdmf.write_function(Bmag)
if rank == 0 and dbg:
print(f"[dbg] wrote XDMF '{args.xdmf}' in {time.perf_counter()-t:.3f} s")
if __name__ == "__main__":
main()
DEBUG:
Info : Reading 'rect_spiral_with_inner_d.msh'...
Info : 354 entities
Info : 87399 nodes
Info : 514648 elements
Info : Done reading 'rect_spiral_with_inner_d.msh'
[dbg] read_from_msh: 2.935 s
[info] cells: total=514172, coil=64634, air=449538
[info] params: N=5, din=0.02, w=0.001, s=0.000305, layers=1, t=3.5e-05, offset=0.1
[dbg] return_loop: via_h=0.00021, bridge_tz=7e-05, z_top=3.5e-05, z_bridge=0.000245
[dbg] z_spiral_mid=1.75e-05, z_bridge_mid=0.00028
[dbg] mu_air=1.256637e-06, mu_cu=1.256637e-06
[dbg] invmu set (air/coils): 0.011 s, n_dofs(Q)=514172
[dbg] built segments3d in 0.000 s (n=23)
L_spiral=5.174500e-01, L_via0=2.100000e-04, L_bridge=9.227743e-03, L_via1=2.100000e-04, L_total=5.270977e-01
[diag] V_coil(mesh) = 1.840e-08 m^3
[diag] L_path(3D) = 5.271e-01 m
[diag] L_theo2d = 5.175e-01 m, V_theo ≈ 1.811e-08 m^3
[dbg] outer BC: facets=15438, dofs=23157
[diag] Jmag = 2.865e+07 A/m^2 (skalierte J-Größe)
[diag] J0_th (I/(w*t)) = 2.857e+07 A/m^2, ratio Jmag/J0_th = 1.003
[dbg] filled J (DG0 vec) in 6.484 s, n_dofs(Wv)=1542516
[dbg] ||J||_L2^2 over COIL = 1.509898e+07
[dbg] Jx min/max on set DOFs = -2.865e+07 / 2.865e+07
[dbg] Jy min/max on set DOFs = -2.865e+07 / 2.865e+07
[dbg] Jz min/max on set DOFs = -2.865e+07 / 2.865e+07
[dbg] built magnetization field M in 0.478 s
[dbg] ||M||_L2^2 over COIL = 6.966018e-03
[diag] ||b_M||_2 = 0.000000e+00 (RHS for M·curl(v))
[diag] ||b_J||_2 = 1.168467e+01 (RHS for J·v)
[warn] b_M is ~0 → switching RHS to ∫ J·v dx(COIL) for this run.
[dbg] assembled A: 609289x609289, ||b||_2 (used) = 1.168467e+01
Residual norms for dolfinx_solve_140392748349328 solve.
0 KSP preconditioned resid norm 2.391738792685e-09 true resid norm 1.168467430546e+01 ||r(i)||/||b|| 1.000000000000e+00
1 KSP preconditioned resid norm 2.115977861520e-09 true resid norm 3.025691613538e+01 ||r(i)||/||b|| 2.589453102791e+00
2 KSP preconditioned resid norm 2.889406258169e-09 true resid norm 2.725998075772e+01 ||r(i)||/||b|| 2.332968814115e+00
3 KSP preconditioned resid norm 3.558588537718e-09 true resid norm 2.740711052062e+01 ||r(i)||/||b|| 2.345560501230e+00
.
.
.
93 KSP preconditioned resid norm 2.078192668759e-05 true resid norm 8.471668680144e+04 ||r(i)||/||b|| 7.250239466397e+03
94 KSP preconditioned resid norm 2.288764633062e-05 true resid norm 9.188390810838e+04 ||r(i)||/||b|| 7.863625951940e+03
95 KSP preconditioned resid norm 2.554314848967e-05 true resid norm 1.022047265633e+05 ||r(i)||/||b|| 8.746904183328e+03
[dbg] linear solve done in 17.914 s
---------- Ergebnisse ----------
I = 1 A
W_m = 6.438404e+00 J
L (Num) = 1.287681e+01 H (12876807.870 µH)
alpha (gauge) = 0.0001
Hinweis: alpha zu klein -> langsamer; zu groß -> leichte Verfälschung (1e-5..3e-4 üblich).
[dbg] energy & L computed in 0.085 s
[dbg] total wall time: 29.302 s
[dbg] wrote XDMF 'solution_fields.xdmf' in 0.249 s
Process finished with exit code 0