From 5694a4e59900002fdc6aabca33a785361a5c4648 Mon Sep 17 00:00:00 2001 From: Max Date: Tue, 3 Nov 2026 21:02:57 +0100 Subject: [PATCH 1/3] Add coils of real wire (Feed3D), the slotted-cylinder tests and R_net of the BUG driver. --- .gitignore | 1 + CHANGELOG.md | 16 ++ benchmarks/rnet_driver.json | 75 +++++ benchmarks/rnet_driver.py | 259 ++++++++++++++++++ docs/src/content/docs/development/field-3d.md | 64 ++++- src/wellbalanced/spatial/__init__.py | 2 + src/wellbalanced/spatial/inductive3d.py | 113 +++++++- src/wellbalanced/spatial/iterative.py | 8 +- src/wellbalanced/tests/unit/test_cutcell.py | 69 +++++ src/wellbalanced/tests/unit/test_feed3d.py | 145 ++++++++++ tutorials/field_3d.ipynb | 109 ++++++++ 11 files changed, 856 insertions(+), 5 deletions(-) create mode 100644 benchmarks/rnet_driver.json create mode 100644 benchmarks/rnet_driver.py create mode 100644 src/wellbalanced/tests/unit/test_feed3d.py diff --git a/.gitignore b/.gitignore index bee74c7..2045252 100644 --- a/.gitignore +++ b/.gitignore @@ -174,3 +174,4 @@ junit/ # Quarto .quarto/ README_files/ +.rnet_cache/ diff --git a/CHANGELOG.md b/CHANGELOG.md index fa54de1..facd293 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,19 @@ format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). ### Added +- `Feed3D`: coils of real wire in the 3D field. The coil current flows along grid edges + across a gap between the terminals of an `ImpedanceConductor3D` (a delta-gap source) and + returns on its surface, so the impedance includes the coil's own loss and proximity + effect; terminals not joined through one piece of metal are refused. Tested on a coaxial + line against its closed form (reactance to 0.13 %, resistance to 4 % at a / h = 4) and on + a ring of wire. +- The slotted-cylinder test of step R4 as unit tests (`test_cutcell.py`): the z-invariant + slotted tube against NGSolve to 0.1 % with cut cells (staircase 20 % off), the finite tube + with closed and through slits to 6 %. +- `benchmarks/rnet_driver.py`: R_net of the BUG driver (helical coil, Faraday shield with + 24 slits, optional EMS with supports) on reduced grids, with the cut cells cached. +- `InductiveField3D(cut_cell_samples=...)`: the sampling of the cut cells' crossings. + - Method cards (methodcard) on the equations of the models: the reactions, the global model (cusps, wall areas, gas temperature, the steady state), the OpenCR rates, the neutralizer (beam, Knudsen profile, pumping, fractions, sources, fixed point), the 0D RF driver @@ -264,6 +277,9 @@ format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). ### Fixed +- The auxiliary-space preconditioner's coarsest multigrid level could be exactly singular + (two coarse nodes reaching one fine node through a gap in a conductor); a relative + diagonal shift of 1e-10 keeps it regular. - The docstrings' `:func:`, `:class:` and `:meth:` roles are code spans, which the API pages render. diff --git a/benchmarks/rnet_driver.json b/benchmarks/rnet_driver.json new file mode 100644 index 0000000..b7f4a81 --- /dev/null +++ b/benchmarks/rnet_driver.json @@ -0,0 +1,75 @@ +{ + "results": [ + { + "variant": "bug", + "h": 0.004, + "r_net": 0.2664038792549869, + "resistances": { + "coil": 0.07089950205981722, + "shield": 0.1955043360682568 + }, + "inductance": 4.11169189604374e-06, + "power_balance": -1.5437805556128836e-07, + "unknowns": 1116650, + "iterations": 758, + "setup_seconds": 443.21280483300507, + "seconds": 382.2669253750064, + "tolerance": 0.0001, + "driver": { + "plasma_radius": 0.1185, + "length": 0.17, + "shield_wall": 0.008, + "slits": 24, + "slit_width": 0.002, + "slit_margin": 0.01, + "coil_radius": 0.1385, + "turns": 6, + "pitch": 0.02, + "wire_radius": 0.004, + "lead_radius": 0.155, + "lead_side": 0.01, + "feed_gap": 0.01, + "ems_distance": null, + "ems_wall": 0.008, + "box_half": 0.17, + "box_end": 0.22, + "frequency": 1000000.0 + } + }, + { + "variant": "coil", + "h": 0.004, + "r_net": 0.07035406989410348, + "resistances": { + "coil": 0.07035458043107835 + }, + "inductance": 6.841778988157454e-06, + "power_balance": 7.256680041978214e-06, + "unknowns": 1152314, + "iterations": 638, + "setup_seconds": 0.4532764999894425, + "seconds": 687.7985254580126, + "tolerance": 1e-05, + "driver": { + "plasma_radius": 0.1185, + "length": 0.17, + "shield_wall": 0.008, + "slits": 24, + "slit_width": 0.002, + "slit_margin": 0.01, + "coil_radius": 0.1385, + "turns": 6, + "pitch": 0.02, + "wire_radius": 0.004, + "lead_radius": 0.155, + "lead_side": 0.01, + "feed_gap": 0.01, + "ems_distance": null, + "ems_wall": 0.008, + "box_half": 0.17, + "box_end": 0.22, + "frequency": 1000000.0 + } + } + ] +} diff --git a/benchmarks/rnet_driver.py b/benchmarks/rnet_driver.py new file mode 100644 index 0000000..b3d334a --- /dev/null +++ b/benchmarks/rnet_driver.py @@ -0,0 +1,259 @@ +"""R_net of an RF driver: the losses of its coil and Faraday shield, computed in 3D. + +Step R4's output. The network resistance ``R_net`` in ``eta = P_plasma / (P_plasma + +1/2 R_net I^2)`` is measured without plasma (0.6 Ohm at BUG; 1.28 Ohm for ELISE's two +upper drivers in series, Zielke et al., NIBS'22, arXiv:2211.10143); at BUG 25 % of it +is the coil's Joule heating and 75 % the Faraday shield's eddy currents (Zielke et al., +J. Phys. D 54, 155202, 2021). Here it is the resistance the coil's terminals see without +plasma, ``Re Z = 2 P / |I|^2`` from `InductiveField3D` with: + +- the coil: a copper `simgeometry.Helix` with two radial leads and an axial return bar, + all `ImpedanceConductor3D`, driven by a `Feed3D` across a gap in the bar (its own + Joule loss and proximity effect included); +- the Faraday shield: a copper tube with axial slits closed at both ends, on the back + plate; +- the housing: a perfectly conducting box (the back plate at z = 0, side walls, a wall + beyond the driver's open end standing in for the expansion chamber), lossless; +- optionally an electromagnetic shield (EMS, ELISE): a copper tube around the coil + with four axial support bars to the back plate. + +Geometry (BUG; ``*`` marks an assumption, to be replaced by IPP's drawings): plasma +radius 0.1185 m, driver length 0.17 m (Zielke et al.); the shield's inner radius at the +plasma, its wall 8 mm* thick (real shields are thinner: a wall below two cells leaks +the field, see below), 24* slits 2 mm* wide ending 10 mm* before its ends; six turns of +radius 0.1385 m (2 cm outside the plasma*, as ``illustrative_hydrogen_driver``), pitch +2 cm*, an 8 mm* copper tube; copper at 5.8e7 S/m, 1 MHz. + +Resolution: the solver needs about two cells across every metal part (an edge longer +than the wall reaches through it, and a feed needs its terminals joined through the +metal), and the losses converge at first order in h at the 3D corners (the slits' ends; +``docs/.../field-3d.md``). The grids here (5, 4 and 3.33 mm) are the reduced ones; the +converged R_net, with the real shield's thickness (about 2 mm) and h about 1 mm (some +10^8 edges), needs a production run on mpiarray and cunumpy. + +Run ``python benchmarks/rnet_driver.py VARIANT H_MM[,H_MM...]`` with ``VARIANT`` one +of ``coil`` (the coil alone), ``bug`` (coil and shield), ``elise`` (with the EMS at +4 cm); rows are added to ``benchmarks/rnet_driver.json``. Needs simgeometry with +``Helix``. +""" + +from __future__ import annotations + +import hashlib +import json +import sys +import time +from dataclasses import asdict, dataclass +from pathlib import Path + +import numpy as np +import simgeometry as sg + +from wellbalanced.spatial import ( + EdgeGrid, + Feed3D, + ImpedanceConductor3D, + InductiveField3D, + inductive3d, +) +from wellbalanced.spatial.cutcell import CutCells + +HERE = Path(__file__).resolve().parent +PATH = HERE / "rnet_driver.json" +CACHE = HERE / ".rnet_cache" # the cut cells, by conductor and grid (not committed) +COPPER = 5.8e7 +GEOMETRY = [""] # the driver's parameters, part of the cache's key +# COCG's relative residual: the slow part of the residual is the gauge's (gradients), +# which leaves the curl and the impedance unchanged (the coax of test_feed3d.py: 1e-4 +# gives the resistance to 2e-7 of the direct solve, the driver coil at 4 mm 1e-5 and +# 1e-8 the same to 4 digits) +TOLERANCE = 1e-4 + + +def cached_cut_cells(original): + """``cut_cells`` with the results kept on disk: they take minutes on these grids.""" + + def cut_cells(grid, shapes, **options): + key = hashlib.sha1( + repr( + ( + GEOMETRY[0], + tuple(grid.lower), + tuple(grid.upper), + grid.cells, + [type(shape).__name__ for shape in shapes], + options, + ) + ).encode() + ).hexdigest()[:16] + path = CACHE / f"{key}.npz" + if path.exists(): + data = np.load(path) + return CutCells(data["edge_air"], data["face_air"], data["face_width"]) + start = time.perf_counter() + cells = original(grid, shapes, **options) + print(f" cut cells in {time.perf_counter() - start:.0f} s", flush=True) + CACHE.mkdir(exist_ok=True) + np.savez(path, **asdict(cells)) + return cells + + return cut_cells + + +inductive3d.cut_cells_of = cached_cut_cells(inductive3d.cut_cells_of) + + +@dataclass(frozen=True) +class Driver: + """The driver's geometry [m] and the frequency [Hz]; ``*`` in the module: assumed.""" + + plasma_radius: float = 0.1185 + length: float = 0.17 + shield_wall: float = 0.008 + slits: int = 24 + slit_width: float = 0.002 + slit_margin: float = 0.010 + coil_radius: float = 0.1385 + turns: int = 6 + pitch: float = 0.02 + wire_radius: float = 0.004 + lead_radius: float = 0.155 # the axial return bar's distance from the axis + lead_side: float = 0.010 # the leads' square cross-section + feed_gap: float = 0.010 + ems_distance: float | None = None # from the coil's centreline + ems_wall: float = 0.008 + box_half: float = 0.17 + box_end: float = 0.22 + frequency: float = 1.0e6 + + +def coil_shapes(d: Driver) -> tuple[list, Feed3D]: + """The helix, its two radial leads and the return bar with its gap; the feed.""" + start = 0.5 * (d.length - d.pitch * d.turns) + end = start + d.pitch * d.turns + helix = sg.Helix((0, 0, start), d.coil_radius, d.pitch, d.turns, d.wire_radius) + s = d.lead_side + radial = d.lead_radius + s / 2 - d.coil_radius + centre = d.coil_radius + radial / 2 + leads = [ + sg.Cuboid((centre, 0, z), (radial, s, s)) for z in (start, end) + ] # from the helix's ends (at +x) out to the bar + middle = 0.5 * (start + end) + bar = sg.Cuboid((d.lead_radius, 0, middle), (s, s, end - start + s)) - sg.Cuboid( + (d.lead_radius, 0, middle), (2 * s, 2 * s, d.feed_gap) + ) + reach = d.feed_gap / 2 + s / 2 + feed = Feed3D( + (d.lead_radius, 0, middle - reach), + (d.lead_radius, 0, middle + reach), + name="feed", + ) + return [helix, *leads, bar], feed + + +def shield_shape(d: Driver): + """The slotted tube on the back plate (z = 0), closed slits.""" + inner, outer = d.plasma_radius, d.plasma_radius + d.shield_wall + shape = sg.Tube((0, 0, 0), (0, 0, d.length), inner, outer) + slit = sg.Cuboid( + (0.5 * (inner + outer), 0, d.length / 2), + (3 * d.shield_wall, d.slit_width, d.length - 2 * d.slit_margin), + ) + return shape - sg.PolarRepeat(slit, d.slits) + + +def ems_shapes(d: Driver) -> list: + """The EMS tube around the coil and four axial support bars to the back plate.""" + inner = d.coil_radius + d.ems_distance + tube = sg.Tube((0, 0, 0), (0, 0, d.length), inner, inner + d.ems_wall) + side = d.ems_wall + supports = [ + sg.Cuboid((0, 0, d.length / 2), (side, side, d.length)).translated( + ( + (inner + side / 2) * np.cos(np.pi / 4 + k * np.pi / 2), + (inner + side / 2) * np.sin(np.pi / 4 + k * np.pi / 2), + 0.0, + ) + ) + for k in range(4) + ] + return [tube, *supports] + + +def model(d: Driver, h: float, variant: str) -> InductiveField3D: + """The field of the driver ``variant`` on a grid of spacing ``h``.""" + cells = round(2 * d.box_half / h) + grid = EdgeGrid( + (-d.box_half, -d.box_half, 0.0), + (d.box_half, d.box_half, d.box_end), + (cells, cells, round(d.box_end / h)), + ) + coil, feed = coil_shapes(d) + conductors = [ImpedanceConductor3D(coil, COPPER, "coil")] + if variant in ("bug", "elise"): + conductors.append(ImpedanceConductor3D([shield_shape(d)], COPPER, "shield")) + if variant == "elise": + conductors.append(ImpedanceConductor3D(ems_shapes(d), COPPER, "ems")) + return InductiveField3D( + grid, + frequency=d.frequency, + feeds=[feed], + impedance_conductors=conductors, + cut_cell_samples=16, + ) + + +def run(variant: str, h: float) -> dict: + """One row: R_net and its parts, the inductance, the cost.""" + d = Driver(ems_distance=0.04) if variant == "elise" else Driver() + if variant == "elise": + d = Driver(ems_distance=0.04, box_half=0.20) + GEOMETRY[0] = repr(asdict(d)) + setup = time.perf_counter() + field = model(d, h, variant) + setup = time.perf_counter() - setup + start = time.perf_counter() + solution = field.solve(1.0, solver="iterative", tolerance=TOLERANCE, max_iterations=4000) + seconds = time.perf_counter() - start + parts = {name: 2 * p for name, p in solution.conductor_losses.items()} # Ohm at 1 A + row = { + "variant": variant, + "h": h, + "r_net": solution.resistance, + "resistances": parts, + "inductance": solution.inductance, + "power_balance": solution.loss / solution.source_power.real - 1, + "unknowns": int(field.free_edges.sum()), + "iterations": field.last_iterative.iterations, + "setup_seconds": setup, + "seconds": seconds, + "tolerance": TOLERANCE, + "driver": asdict(d), + } + print( + f"{variant} h={1e3 * h:.2f} mm: R_net {row['r_net']:.4f} Ohm" + f" ({', '.join(f'{k} {v:.4f}' for k, v in parts.items())})," + f" L {1e6 * row['inductance']:.3f} uH, {row['unknowns']} unknowns," + f" {row['iterations']} iterations, {setup:.0f} + {seconds:.0f} s", + flush=True, + ) + return row + + +def main(arguments: list[str]) -> None: + """Run one variant on the given grids and merge its rows into the JSON file.""" + variant = arguments[0] if arguments else "bug" + spacings = [ + float(v) * 1e-3 + for v in (arguments[1] if len(arguments) > 1 else "5").split(",") + ] + new = [run(variant, h) for h in spacings] + old = json.loads(PATH.read_text())["results"] if PATH.exists() else [] + keys = {(r["variant"], r["h"]) for r in new} + kept = [r for r in old if (r["variant"], r["h"]) not in keys] + rows = sorted(kept + new, key=lambda r: (r["variant"], -r["h"])) + PATH.write_text(json.dumps({"results": rows}, indent=1, default=float) + "\n") + + +if __name__ == "__main__": + main(sys.argv[1:]) diff --git a/docs/src/content/docs/development/field-3d.md b/docs/src/content/docs/development/field-3d.md index 88b2cbd..6b9b32b 100644 --- a/docs/src/content/docs/development/field-3d.md +++ b/docs/src/content/docs/development/field-3d.md @@ -243,6 +243,62 @@ Loss against NGSolve's surface impedance on fitted meshes, `tube3d` of - On a quiet machine the iterative solve takes 25–50 iterations and 3–11 s for 2·10^5 to 5·10^5 unknowns (two threads). +## Coils of real wire: the feed + +A coil of copper tube is an `ImpedanceConductor3D` driven by a `Feed3D`: the coil current +flows along grid edges across a gap between its terminals (a delta-gap source), from inside +the metal of one terminal to inside the other, and returns on the metal's surface. The +impedance at the terminals then includes the coil's own Joule loss and its proximity +effect. The terminals must lie in one connected piece of metal on the grid (the nodes +joined by edges that are not unknowns), which needs about two cells across the wire; +otherwise the feed is refused. The current's gradient part (the charge at the terminals) +returns through the metal, so the source is compatible with the gauge (the result changes +by below `1e-4` between `eps = 1e-4` and `1e-6`). + +**The coax (tests).** A coaxial line of copper, a = 4 mm and b = 12 mm, shorted by the +box's wall and fed across a 2 mm gap, against `Z' = i omega mu0 / (2 pi) ln(b / a) + Z_s / +(2 pi) (1 / a + 1 / b)` (the second term is the skin's resistance and as much internal +reactance), from two lengths: + +| h | 2 mm | 1 mm | 0.67 mm | 0.5 mm | +| --- | --- | --- | --- | --- | +| reactance X' | +0.96 % | +0.13 % | +0.10 % | +0.04 % | +| resistance R' | −2.6 % | −3.9 % | −0.5 % | −1.7 % | + +The reactance converges; the resistance is within a few per cent and not monotone: the +surface loss takes each cut face's mean `H`, which varies as `1 / r` across a cell next to +the wire (first order in `h / a`). A ring of wire (R / a = 5) gives the thin-ring +inductance to 1 % (the box) and a resistance 18 % above `R_s R / a`, the ring's own +proximity effect. + +**The iterative solver** converges more slowly with a feed (60 iterations at 2 mm and 150 +at 1 mm for the coax, 600–800 for the driver below, against 25–50 for a current sheet), but +the slow part of the residual is the gauge's: at a residual of `1e-4` the coax's +resistance is within `2e-7` of the direct solve. + +## R_net of the BUG driver + +`benchmarks/rnet_driver.py`: the six-turn coil (radius 0.1385 m, pitch 2 cm, an 8 mm +copper tube, simgeometry's `Helix`) with radial leads and a return bar fed across a 1 cm +gap; the Faraday shield (inner radius 0.1185 m, 0.17 m long, 24 axial slits 2 mm wide +ending 1 cm before its ends, `tube - PolarRepeat(slit, 24)`) on the back plate; a +perfectly conducting housing 3 cm outside the coil and 5 cm beyond the driver's open end; +copper at 1 MHz. The pitch, the tube's size, the number and width of the slits and the +housing are assumptions, not IPP's drawings. The shield's wall is 8 mm, thicker than a +real one, because a wall below about two cells lets the field through (an edge reaches +across it). + +| h | unknowns | coil | shield | R_net | L | iterations | cut cells + solve | +| --- | --- | --- | --- | --- | --- | --- | --- | +| 4 mm | 1.12 M | 0.071 Ω | 0.196 Ω | 0.266 Ω | 4.11 µH | 758 | 443 + 382 s | + +The coil alone (no shield) is 0.070 Ω. The split, 27 % in the coil and 73 % in the +shield, is the measured one (25 % and 75 %, Zielke et al. 2021), and R_net is 44 % of the +measured 0.6 Ω, which also has the matching network and the lines. 4 mm is the coarsest +grid on which the 8 mm tube is one piece of metal. The converged value needs the real +shield (about 2 mm thick, so h ≲ 1 mm there) and a resolved wire (h ≲ a / 4): about 10^8 +edges on a uniform grid, a production run on mpiarray and cunumpy, or local refinement. + ## What it costs - **Yee grid, direct:** the direct solve takes 40–60 s for 65 000 unknowns (24 × 24 × 40 @@ -318,4 +374,10 @@ need multipatch domains there. so they cannot refine the 3D field. - The iterative solver on mpiarray and cunumpy (it is all sparse products and a multigrid on a Cartesian grid). -- The helical coil, the Faraday shield of BUG and ELISE, and R_net. +- Done in the third round: coils of real wire (`Feed3D`), the helical coil and the + slotted Faraday shield (simgeometry's `Helix` and `PolarRepeat`), and R_net of BUG on + reduced grids. +- The converged R_net (the real shield's thickness, a resolved wire), and ELISE's with the + EMS and its supports: production runs. +- The resistance's first-order error where `H` varies across a cut face (a wire): an + extrapolation of `H` to the surface in the loss. diff --git a/src/wellbalanced/spatial/__init__.py b/src/wellbalanced/spatial/__init__.py index 892c507..9ddf6c9 100644 --- a/src/wellbalanced/spatial/__init__.py +++ b/src/wellbalanced/spatial/__init__.py @@ -44,6 +44,7 @@ from wellbalanced.spatial.inductive3d import ( Coil3D, Conductor3D, + Feed3D, ImpedanceConductor3D, InductiveField3D, InductiveSolution3D, @@ -72,6 +73,7 @@ "CoupledState", "DriverCoupling", "EdgeGrid", + "Feed3D", "ImpedanceConductor3D", "Fluid2D", "Fluid2DState", diff --git a/src/wellbalanced/spatial/inductive3d.py b/src/wellbalanced/spatial/inductive3d.py index b6c60bd..dba7db4 100644 --- a/src/wellbalanced/spatial/inductive3d.py +++ b/src/wellbalanced/spatial/inductive3d.py @@ -18,7 +18,11 @@ ``j = C^T (T L*)``, so ``G^T j = 0`` exactly: the coil current is divergence-free on the grid whatever the shape. A solenoid is a cylinder with ``T`` along its axis; a rectangular loop is a box with ``T`` normal to the loop; the sheet sits on the region's -staircased surface. +staircased surface. A coil of real wire (copper tube, with its own losses and proximity +effect) is an `ImpedanceConductor3D` driven by a `Feed3D`: the current along grid +edges across a gap between the coil's terminals (a delta-gap source), which returns +through the metal; the terminals must lie in one connected piece of metal on the grid +(about two cells across the wire). Conductors: `Conductor3D`, a conductivity on the edges inside shapes (the midpoint rule, staircased; or a sub-cell average, ``samples``), and a background @@ -77,6 +81,7 @@ __all__ = [ "Coil3D", "Conductor3D", + "Feed3D", "ImpedanceConductor3D", "InductiveField3D", "InductiveSolution3D", @@ -110,6 +115,33 @@ class Coil3D: name: str = "" +@dataclass(frozen=True) +class Feed3D: + """A current source along grid edges: the terminals of a coil of real wire. + + The coil current flows along the edges of the grid from the node nearest ``start`` + to the node nearest ``end``, axis by axis in ``order``. Both terminals must lie + in the metal of the same `ImpedanceConductor3D` (or of the same Dirichlet + wall): the current enters that conductor at one end of the path, flows on its + surface, and leaves it at the other, so the coil's own Joule loss and its proximity + effect are in the solution. A helical coil with its two leads and a gap between + their ends: the feed bridges the gap, from inside one lead to inside the other + (a delta-gap source; ``V = -i omega a`` summed along the path is the terminal + voltage). The edges of the path inside the metal are not unknowns and drop out. + + Attributes: + start: The first terminal (a point inside the metal) [m]. + end: The second terminal [m]. + order: The axes the path follows, in turn. + name: A label. + """ + + start: tuple[float, float, float] + end: tuple[float, float, float] + order: tuple[int, int, int] = (0, 1, 2) + name: str = "" + + @dataclass(frozen=True) class Conductor3D: """A conductor of uniform conductivity: the union of shapes. @@ -328,6 +360,9 @@ class InductiveField3D: grid: The edge grid (its boundaries are the field's). frequency: The frequency [Hz]. coils: The coils, all carrying the coil current. + feeds: Current sources along grid edges (`Feed3D`), also carrying the + coil current: the terminals of coils of real wire, modelled as + `ImpedanceConductor3D`. conductors: Conductors of uniform conductivity (must not overlap). conductivity: A background conductivity [S/m] as a function of ``(3, n)`` points, e.g. a cold plasma's complex ``sigma``; added outside the @@ -346,6 +381,10 @@ class InductiveField3D: true area otherwise). min_sine: The bound of ``sin(theta)`` in the cut faces' surface widths (`wellbalanced.spatial.cutcell.cut_cells`). + cut_cell_samples: The sampling of the conductors' crossings along each edge + and line of the cut cells (simgeometry's ``yee_fractions``), for shapes + without analytic crossings (Booleans, helices); fewer is faster, and + features thinner than a ``samples``-th of a cell can be missed. """ # --- sources and materials ----------------------------------------------------- @@ -390,6 +429,72 @@ def _coil_source(self) -> np.ndarray: t += weight * coil.turns_per_length * direction[normals] return np.asarray(_curl(grid).T @ (t * grid.dual_lengths)) + def _feed_edges(self, feed: Feed3D) -> tuple[np.ndarray, np.ndarray, list]: + """The edges of a feed's path, their signs and the two terminal nodes.""" + grid = self.grid + nodes = [] + for point in (feed.start, feed.end): + index = np.rint( + (np.asarray(point, dtype=float) - grid.lower) / grid.spacing + ) + index = index.astype(int) + if np.any(index < 0) or np.any(index >= np.asarray(grid.nodes)): + raise ValueError(f"feed {feed.name!r}: a terminal is outside the grid") + nodes.append(index) + if sorted(feed.order) != [0, 1, 2]: + raise ValueError("a feed's order is a permutation of the three axes") + edges, signs = [], [] + current = nodes[0].copy() + for axis in feed.order: + step = 1 if nodes[1][axis] >= current[axis] else -1 + while current[axis] != nodes[1][axis]: + low = current.copy() + if step < 0: + low[axis] -= 1 + edges.append( + int(grid._edge_index(axis, [np.array([i]) for i in low])[0]) + ) + signs.append(float(step)) + current[axis] += step + if not edges: + raise ValueError(f"feed {feed.name!r}: its terminals are the same node") + flat = [int(np.ravel_multi_index(tuple(n), grid.nodes)) for n in nodes] + return np.asarray(edges), np.asarray(signs), flat + + def _feed_source(self) -> np.ndarray: + """``j`` of the feeds for 1 A [A]: the current along each path's edges. + + Raises if a feed's terminals are not in the same connected piece of metal (the + nodes joined by edges that are not unknowns): then its current would have no + return path on the grid. + """ + grid = self.grid + source = np.zeros(grid.n_edges) + if not self.feeds: + return source + import scipy.sparse.csgraph as csgraph + + gradient = abs(self._gradient_matrix()) + metal = sp.diags((~self.free_edges).astype(float)) + graph = (gradient.T @ metal @ gradient).tocsr() + _, piece = csgraph.connected_components(graph, directed=False) + isolated = np.asarray(graph.sum(axis=1)).ravel() == 0 + for feed in self.feeds: + edges, signs, (first, last) = self._feed_edges(feed) + if isolated[first] or isolated[last] or piece[first] != piece[last]: + raise ValueError( + f"feed {feed.name!r}: its terminals must be inside the metal of" + " one conductor (or one Dirichlet wall)" + ) + np.add.at(source, edges, signs) + return source + + def _gradient_matrix(self) -> sp.csr_matrix: + grid = self.grid + if not hasattr(grid, "_gradient_matrix"): + grid._gradient_matrix = grid.gradient() # type: ignore[attr-defined] + return grid._gradient_matrix # type: ignore[attr-defined] + def _edge_conductivity( self, background: PointField | None ) -> tuple[np.ndarray, dict[str, np.ndarray]]: @@ -438,6 +543,7 @@ def __init__( *, frequency: float, coils: Sequence[Coil3D] = (), + feeds: Sequence[Feed3D] = (), conductors: Sequence[Conductor3D] = (), conductivity: PointField | None = None, samples: int = 1, @@ -446,6 +552,7 @@ def __init__( cut_cells: bool = True, min_air_fraction: float = 0.01, min_sine: float = 0.1, + cut_cell_samples: int = 64, ) -> None: """Sample the coils and conductors on the grid and assemble the operator.""" if frequency <= 0.0: @@ -467,7 +574,7 @@ def __init__( raise ValueError("min_air_fraction must be in (0, 1]") self.min_air_fraction = float(min_air_fraction) self._surfaces = [ - cut_cells_of(grid, c.shapes, min_sine=min_sine) + cut_cells_of(grid, c.shapes, min_sine=min_sine, samples=cut_cell_samples) if self.cut_cells else staircased_cells(grid, c.shapes) for c in self.impedance_conductors @@ -486,6 +593,8 @@ def __init__( volume |= edges if np.any(volume & (self.edge_air < 1.0)): raise ValueError("volume conductors must not touch impedance conductors") + self.feeds = tuple(feeds) + self._unit_source = self._unit_source + self._feed_source() self._factorizations: dict[float, Any] = {} self._preconditioners: dict[float, Any] = {} self._order: np.ndarray | None = None diff --git a/src/wellbalanced/spatial/iterative.py b/src/wellbalanced/spatial/iterative.py index 2694f1d..50225cc 100644 --- a/src/wellbalanced/spatial/iterative.py +++ b/src/wellbalanced/spatial/iterative.py @@ -107,9 +107,13 @@ def __init__( lattice, active = coarse_lattice, coarse_active self.levels.append((matrix, 1.0 / diagonal, None)) # unknowns no edge reaches (e.g. a component of the nodal field none of a - # node's edges is along) have empty rows: a unit diagonal keeps them at zero + # node's edges is along) have empty rows: a unit diagonal keeps them at zero; + # a coarse node whose active fine nodes are those of another (two coarse nodes + # reaching one fine node through a gap in a conductor) makes the Galerkin + # operator singular: a relative shift of 1e-10 on the diagonal keeps it regular empty = np.asarray(abs(matrix).sum(axis=1)).ravel() == 0 - self.coarse = spla.splu(sp.csc_matrix(matrix + sp.diags(empty.astype(float)))) + shift = 1e-10 * np.abs(matrix.diagonal()) + empty + self.coarse = spla.splu(sp.csc_matrix(matrix + sp.diags(shift))) def _cycle(self, level: int, rhs: np.ndarray) -> np.ndarray: matrix, inverse, prolongation = self.levels[level] diff --git a/src/wellbalanced/tests/unit/test_cutcell.py b/src/wellbalanced/tests/unit/test_cutcell.py index ae113f8..359b296 100644 --- a/src/wellbalanced/tests/unit/test_cutcell.py +++ b/src/wellbalanced/tests/unit/test_cutcell.py @@ -184,3 +184,72 @@ def test_volume_conductors_must_not_touch_impedance_conductors(): ) with pytest.raises(ValueError, match="min_air_fraction"): InductiveField3D(grid, frequency=1e3, min_air_fraction=0.0) + + +# The slotted cylinder of benchmarks/slotted_cylinder_case.py (step R4's test) at 500 kHz +# (skin depth 0.71 mm in the 6 mm wall), against NGSolve with the same surface impedance +# on curved meshes fitted to the tube and its 4 mm slits, finest mesh and order +# (benchmarks/slotted_cylinder_ngsolve_sibc.py): the losses at 1 A. +NGSOLVE_2D = {"none": 0.5820426537591570, "through": 0.6682538180678069} # W/m +NGSOLVE_3D = { # W, one eighth + "none": 0.00923019992605531, + "closed": 0.013396949945715761, + "through": 0.008593391259681745, +} + + +def slotted_tube(slits, half_length): + """The 30/36 mm tube with four axial 4 mm slits, closed 10 mm before its ends.""" + shape = Tube((0, 0, -half_length), (0, 0, half_length), 0.030, 0.036) + if slits == "none": + return shape + reach = 4 * half_length if slits == "through" else 0.040 + return ( + shape + - Cuboid((0, 0, 0), (1.0, 0.004, 2 * reach)) + - Cuboid((0, 0, 0), (0.004, 1.0, 2 * reach)) + ) + + +@pytest.mark.parametrize("slits", ["none", "through"]) +def test_slotted_tube_loss_matches_ngsolve_in_2d(slits): + """z-invariant: cut cells to 0.1 % at h = 2.5 mm (slits 1.6 cells), staircase 20 % off.""" + shapes = [slotted_tube(slits, 1.0)] + _, cut, scale = quarter(24, shapes, 5e5) + _, stair, _ = quarter(24, shapes, 5e5, cut_cells=False) + assert scale * cut.loss == pytest.approx(NGSOLVE_2D[slits], rel=1e-3) + assert scale * stair.loss / NGSOLVE_2D[slits] - 1 > 0.15 + + +def test_slotted_tube_losses_match_ngsolve_in_3d(): + """The finite tube, one eighth, h = 2.5 mm: within 6 %, and the slits' effect. + + The remainder is at the 3D corners (the tube's and the slits' ends), first order; + closed slits, as in a Faraday shield, raise the loss by 45 % (NGSolve) and the cut + cells give 36 %. + """ + h = 0.0025 + grid = EdgeGrid( + (0, 0, 0), + (BOX, BOX, 0.1), + (24, 24, 40), + ("dirichlet", "dirichlet", ("natural", "dirichlet")), + ) + coil = Coil3D(Cuboid((0, 0, 0), (2 * COIL, 2 * COIL, 0.16)), (0, 0, 1), TURNS) + losses = {} + for slits in NGSOLVE_3D: + model = InductiveField3D( + grid, + frequency=5e5, + coils=[coil], + impedance_conductors=[ + ImpedanceConductor3D([slotted_tube(slits, 0.05)], SIGMA, "tube") + ], + ) + solution = model.solve(1.0, solver="iterative") + losses[slits] = solution.loss + assert solution.loss == pytest.approx(NGSOLVE_3D[slits], rel=0.06), slits + assert h == grid.spacing[0] + effect = losses["closed"] / losses["none"] + assert effect == pytest.approx(NGSOLVE_3D["closed"] / NGSOLVE_3D["none"], rel=0.08) + assert losses["through"] < losses["none"] < losses["closed"] diff --git a/src/wellbalanced/tests/unit/test_feed3d.py b/src/wellbalanced/tests/unit/test_feed3d.py new file mode 100644 index 0000000..509f776 --- /dev/null +++ b/src/wellbalanced/tests/unit/test_feed3d.py @@ -0,0 +1,145 @@ +"""Coils of real wire: a feed (delta-gap source) driving surface-impedance conductors. + +Closed form: a coaxial line of surface-impedance conductors, inner radius ``a`` and outer +``b``, has per metre ``L' = mu0 / (2 pi) ln(b / a)`` outside the metal and the surface +impedance's ``Z_s / (2 pi) (1 / a + 1 / b)`` (``R' = Re``, and as much again in +``omega L'``: the skin's internal inductance). The line is shorted at one end by a +perfectly conducting wall and fed across a gap at the other; the difference of two +lengths removes the ends. +""" + +from __future__ import annotations + +import numpy as np +import pytest +from scipy.constants import mu_0 as MU_0 + +from wellbalanced.spatial import ( + EdgeGrid, + Feed3D, + ImpedanceConductor3D, + InductiveField3D, +) + +simgeometry = pytest.importorskip("simgeometry") +Cuboid, Cylinder, Torus, Tube = ( + simgeometry.Cuboid, + simgeometry.Cylinder, + simgeometry.Torus, + simgeometry.Tube, +) + +COPPER = 5.8e7 +INNER, OUTER, OUTER_WALL, HALF, GAP = 0.004, 0.012, 0.014, 0.015, 0.002 + + +def coax(length, h, frequency=1e6, regularization=1e-6, reverse=False): + """A coax along z from the fed end (z = 0) to the shorted end (z = length).""" + cells = round(2 * HALF / h) + grid = EdgeGrid( + (-HALF, -HALF, 0.0), (HALF, HALF, length), (cells, cells, round(length / h)) + ) + inner = Cylinder((0, 0, GAP), (0, 0, length + 1.0), INNER) + outer = Tube((0, 0, -1.0), (0, 0, length + 1.0), OUTER, OUTER_WALL) + # from the wall at z = 0 (joined to the outer conductor) into the inner one + ends = ((0, 0, 0), (0, 0, GAP + 0.003)) + feed = Feed3D(*(ends[::-1] if reverse else ends), name="feed") + return InductiveField3D( + grid, + frequency=frequency, + feeds=[feed], + impedance_conductors=[ + ImpedanceConductor3D([inner], COPPER, "inner"), + ImpedanceConductor3D([outer], COPPER, "outer"), + ], + regularization=regularization, + ) + + +def test_coax_impedance_per_length_matches_the_closed_form(): + h, frequency = 0.001, 1e6 + short, long = ( + coax(length, h, frequency).solve(1.0, solver="iterative") + for length in (0.008, 0.016) + ) + per_length = (long.impedance - short.impedance) / 0.008 + surface = ImpedanceConductor3D((), COPPER).surface_impedance(frequency) + expected = 2j * np.pi * frequency * MU_0 / (2 * np.pi) * np.log( + OUTER / INNER + ) + surface / (2 * np.pi) * (1 / INNER + 1 / OUTER) + # the coil's power is the conductors' losses, to round-off + assert long.loss == pytest.approx(long.source_power.real, rel=1e-8) + assert set(long.conductor_losses) == {"inner", "outer"} + # inductance to 0.3 %, resistance to 5 % (first order: H is a face average + # where it varies as 1 / r, at a / h = 4) + assert per_length.imag == pytest.approx(expected.imag, rel=3e-3) + assert per_length.real == pytest.approx(expected.real, rel=0.05) + + +def test_feed_is_compatible_and_the_gauge_does_not_matter(): + """The feed's current returns through the metal: no gradient of the source.""" + results = [ + coax(0.008, 0.002, regularization=eps).solve(1.0).impedance + for eps in (1e-4, 1e-6) + ] + assert results[0] == pytest.approx(results[1], rel=1e-4) + # the current the other way round: the same impedance + reverse = coax(0.008, 0.002, reverse=True).solve(1.0).impedance + assert reverse == pytest.approx(results[1], rel=1e-10) + + +def test_ring_of_wire(): + """A ring of copper wire with a feed across a gap: L and R of the right size.""" + radius, wire, half, gap, h = 0.02, 0.004, 0.05, 0.002, 0.004 + cells = round(2 * half / h) + grid = EdgeGrid((-half,) * 3, (half,) * 3, (cells,) * 3) + ring = Torus((0, 0, 0), radius, wire) - Cuboid((radius, 0, 0), (0.02, gap, 0.02)) + feed = Feed3D( + (radius, -gap / 2 - 2.5e-3, 0), (radius, gap / 2 + 2.5e-3, 0), order=(1, 0, 2) + ) + model = InductiveField3D( + grid, + frequency=1e6, + feeds=[feed], + impedance_conductors=[ImpedanceConductor3D([ring], COPPER, "ring")], + ) + solution = model.solve(1.0, solver="iterative") + # thin-ring closed forms (a / R = 0.2): the curvature (+18 % in R) and the box + # (-1 % in L) aside + inductance = MU_0 * radius * (np.log(8 * radius / wire) - 2) + skin = ImpedanceConductor3D((), COPPER).surface_impedance(1e6).real + assert solution.inductance == pytest.approx(inductance, rel=0.03) + assert solution.resistance == pytest.approx(skin * radius / wire, rel=0.25) + assert solution.conductor_losses["ring"] == pytest.approx(solution.loss) + + +def test_feed_errors(): + grid = EdgeGrid((0, 0, 0), (1, 1, 1), (4, 4, 4)) + block = ImpedanceConductor3D([Cuboid((0.5, 0.5, 0.5), (0.6, 0.6, 0.6))], COPPER) + with pytest.raises(ValueError, match="outside the grid"): + InductiveField3D(grid, frequency=1.0, feeds=[Feed3D((0, 0, 0), (2, 0, 0))]) + with pytest.raises(ValueError, match="same node"): + InductiveField3D( + grid, frequency=1.0, feeds=[Feed3D((0.5, 0.5, 0.5), (0.5, 0.5, 0.5))] + ) + with pytest.raises(ValueError, match="permutation"): + InductiveField3D( + grid, frequency=1.0, feeds=[Feed3D((0, 0, 0), (1, 1, 1), order=(0, 0, 1))] + ) + # from the metal block into the air: no return path + fine = EdgeGrid((0, 0, 0), (1, 1, 1), (8, 8, 8)) + with pytest.raises(ValueError, match="inside the metal"): + InductiveField3D( + fine, + frequency=1.0, + feeds=[Feed3D((0.5, 0.5, 0.5), (0.5, 0.5, 0.875))], + impedance_conductors=[block], + ) + # from the block to the Dirichlet wall: two pieces of metal + with pytest.raises(ValueError, match="inside the metal"): + InductiveField3D( + fine, + frequency=1.0, + feeds=[Feed3D((0.5, 0.5, 0.5), (0.5, 0.5, 1.0))], + impedance_conductors=[block], + ) diff --git a/tutorials/field_3d.ipynb b/tutorials/field_3d.ipynb index 9f51fba..11778be 100644 --- a/tutorials/field_3d.ipynb +++ b/tutorials/field_3d.ipynb @@ -88,6 +88,115 @@ "$\\sigma$ (`conductivity=`) add their losses to the solution: `solution.loss`,\n", "`solution.conductor_losses` and `solution.resistance`." ] + }, + { + "cell_type": "markdown", + "id": "7fb27b941602401d91542211134fc71a", + "metadata": {}, + "source": "## Coils of real wire\n\nA coil made of copper tube is an `ImpedanceConductor3D` (its inside drops out, Leontovich's\n$E_t = Z_s H \\times n$ on its surface) driven by a `Feed3D`: the coil current flows along\ngrid edges across a gap between the coil's terminals, from inside the metal of one\nterminal to inside the other, and on the metal's surface everywhere else. The impedance\nthe terminals see then includes the coil's own Joule loss and its proximity effect.\n\nThe check is a coaxial line of copper, inner radius $a$ and outer $b$, shorted at the far\nend and fed across a gap at the near one. Per metre,\n$Z' = i\\omega\\frac{\\mu_0}{2\\pi}\\ln\\frac{b}{a} + \\frac{Z_s}{2\\pi}\\left(\\frac1a + \\frac1b\\right)$;\nthe difference of two lengths removes the ends." + }, + { + "cell_type": "code", + "execution_count": null, + "id": "acae54e37e7d407bbb7b55eff062a284", + "metadata": {}, + "outputs": [], + "source": [ + "from simgeometry import Cuboid, Tube\n", + "\n", + "from wellbalanced.spatial import Feed3D, ImpedanceConductor3D\n", + "\n", + "copper, frequency, a, b, gap = 5.8e7, 1e6, 0.004, 0.012, 0.002\n", + "\n", + "\n", + "def coax(length, h):\n", + " \"\"\"The line from z = 0 (fed) to z = length (shorted by the box's wall).\"\"\"\n", + " n = round(0.03 / h)\n", + " grid = EdgeGrid(\n", + " (-0.015, -0.015, 0), (0.015, 0.015, length), (n, n, round(length / h))\n", + " )\n", + " inner = Cylinder((0, 0, gap), (0, 0, length + 1), a)\n", + " outer = Tube((0, 0, -1), (0, 0, length + 1), b, b + 0.002)\n", + " return InductiveField3D(\n", + " grid,\n", + " frequency=frequency,\n", + " feeds=[Feed3D((0, 0, 0), (0, 0, gap + 0.003))],\n", + " impedance_conductors=[\n", + " ImpedanceConductor3D([inner], copper, \"inner\"),\n", + " ImpedanceConductor3D([outer], copper, \"outer\"),\n", + " ],\n", + " )\n", + "\n", + "\n", + "z_s = ImpedanceConductor3D((), copper).surface_impedance(frequency)\n", + "exact = 2j * np.pi * frequency * mu_0 / (2 * np.pi) * np.log(b / a) + z_s / (\n", + " 2 * np.pi\n", + ") * (1 / a + 1 / b)\n", + "rows = []\n", + "for h in (0.002, 0.001):\n", + " short, long = (\n", + " coax(length, h).solve(1.0, solver=\"iterative\") for length in (0.008, 0.016)\n", + " )\n", + " z = (long.impedance - short.impedance) / 0.008\n", + " rows.append((h, z.real / exact.real - 1, z.imag / exact.imag - 1))\n", + "xr.Dataset(\n", + " {\n", + " \"resistance_error\": (\"h\", [r[1] for r in rows], {\"long_name\": \"error of R'\"}),\n", + " \"reactance_error\": (\"h\", [r[2] for r in rows], {\"long_name\": \"error of X'\"}),\n", + " },\n", + " coords={\"h\": (\"h\", [r[0] for r in rows], {\"units\": \"m\"})},\n", + ").to_dataframe()" + ] + }, + { + "cell_type": "markdown", + "id": "9a63283cbaf04dbcab1f6479b197f3a8", + "metadata": {}, + "source": "The reactance converges; the resistance is within a few per cent and first order: the\nsurface loss takes each cut face's average $H$, which here varies as $1/r$ across the cell.\n\n## Slotted cylinder\n\nThe test that decided the discretisation: a copper-like tube (6 mm wall) with four axial\n4 mm slits in a solenoid, at 500 kHz. In the z-invariant cut the cut cells match NGSolve\n(Nédélec elements on meshes fitted to the slits, the same surface impedance) to below\n0.1 %, while staircased cells are 20 % off however fine the grid; the finite tube in 3D is\nin [3D field discretisation](/wellbalanced/development/field-3d/)." + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8dd0d8092fe74a7c96281538738b07e2", + "metadata": {}, + "outputs": [], + "source": [ + "ngsolve = 0.6682538180678069 # W/m at 1 A, the whole cross-section\n", + "\n", + "\n", + "def slotted(cells, cut):\n", + " \"\"\"A quarter of the cross-section, z-periodic, two cells long.\"\"\"\n", + " h = 0.06 / cells\n", + " grid = EdgeGrid(\n", + " (0, 0, 0),\n", + " (0.06, 0.06, 2 * h),\n", + " (cells, cells, 2),\n", + " (\"dirichlet\",) * 2 + (\"periodic\",),\n", + " )\n", + " tube = (\n", + " Tube((0, 0, -1), (0, 0, 1), 0.030, 0.036)\n", + " - Cuboid((0, 0, 0), (1, 0.004, 4))\n", + " - Cuboid((0, 0, 0), (0.004, 1, 4))\n", + " )\n", + " model = InductiveField3D(\n", + " grid,\n", + " frequency=5e5,\n", + " coils=[Coil3D(Cuboid((0, 0, 0), (0.09, 0.09, 4)), (0, 0, 1), 100.0)],\n", + " impedance_conductors=[ImpedanceConductor3D([tube], 1e6)],\n", + " cut_cells=cut,\n", + " )\n", + " return 4 * model.solve(1.0).loss / (2 * h) / ngsolve - 1\n", + "\n", + "\n", + "cells = [12, 24, 48]\n", + "xr.Dataset(\n", + " {\n", + " \"cut_cells\": (\"cells\", [slotted(n, True) for n in cells]),\n", + " \"staircase\": (\"cells\", [slotted(n, False) for n in cells]),\n", + " },\n", + " coords={\"cells\": cells},\n", + ").to_dataframe()" + ] } ], "metadata": { From c21a0a2f57eed429f3ad5e9afd1ef8dcf5f03a9e Mon Sep 17 00:00:00 2001 From: Max Date: Tue, 3 Nov 2026 21:28:05 +0100 Subject: [PATCH 2/3] Add R_net at 3.4 mm and format the benchmark. --- benchmarks/rnet_driver.json | 36 +++++++++++++++++++ benchmarks/rnet_driver.py | 4 ++- docs/src/content/docs/development/field-3d.md | 9 +++-- 3 files changed, 45 insertions(+), 4 deletions(-) diff --git a/benchmarks/rnet_driver.json b/benchmarks/rnet_driver.json index b7f4a81..7055877 100644 --- a/benchmarks/rnet_driver.json +++ b/benchmarks/rnet_driver.json @@ -36,6 +36,42 @@ "frequency": 1000000.0 } }, + { + "variant": "bug", + "h": 0.0034, + "r_net": 0.2855851383706187, + "resistances": { + "coil": 0.0679426879710528, + "shield": 0.21764246298162943 + }, + "inductance": 4.102877478409794e-06, + "power_balance": 4.405713682764656e-08, + "unknowns": 1830128, + "iterations": 860, + "setup_seconds": 418.28152704099193, + "seconds": 1394.2107854999922, + "tolerance": 0.0001, + "driver": { + "plasma_radius": 0.1185, + "length": 0.17, + "shield_wall": 0.008, + "slits": 24, + "slit_width": 0.002, + "slit_margin": 0.01, + "coil_radius": 0.1385, + "turns": 6, + "pitch": 0.02, + "wire_radius": 0.004, + "lead_radius": 0.155, + "lead_side": 0.01, + "feed_gap": 0.01, + "ems_distance": null, + "ems_wall": 0.008, + "box_half": 0.17, + "box_end": 0.22, + "frequency": 1000000.0 + } + }, { "variant": "coil", "h": 0.004, diff --git a/benchmarks/rnet_driver.py b/benchmarks/rnet_driver.py index b3d334a..3807ad9 100644 --- a/benchmarks/rnet_driver.py +++ b/benchmarks/rnet_driver.py @@ -213,7 +213,9 @@ def run(variant: str, h: float) -> dict: field = model(d, h, variant) setup = time.perf_counter() - setup start = time.perf_counter() - solution = field.solve(1.0, solver="iterative", tolerance=TOLERANCE, max_iterations=4000) + solution = field.solve( + 1.0, solver="iterative", tolerance=TOLERANCE, max_iterations=4000 + ) seconds = time.perf_counter() - start parts = {name: 2 * p for name, p in solution.conductor_losses.items()} # Ohm at 1 A row = { diff --git a/docs/src/content/docs/development/field-3d.md b/docs/src/content/docs/development/field-3d.md index 6b9b32b..3196ef9 100644 --- a/docs/src/content/docs/development/field-3d.md +++ b/docs/src/content/docs/development/field-3d.md @@ -291,10 +291,13 @@ across it). | h | unknowns | coil | shield | R_net | L | iterations | cut cells + solve | | --- | --- | --- | --- | --- | --- | --- | --- | | 4 mm | 1.12 M | 0.071 Ω | 0.196 Ω | 0.266 Ω | 4.11 µH | 758 | 443 + 382 s | +| 3.4 mm | 1.83 M | 0.068 Ω | 0.218 Ω | 0.286 Ω | 4.10 µH | 860 | 418 + 1394 s | -The coil alone (no shield) is 0.070 Ω. The split, 27 % in the coil and 73 % in the -shield, is the measured one (25 % and 75 %, Zielke et al. 2021), and R_net is 44 % of the -measured 0.6 Ω, which also has the matching network and the lines. 4 mm is the coarsest +The coil alone (no shield) is 0.070 Ω at 4 mm. From 4 to 3.4 mm the coil's resistance +falls by 4 % and the shield's rises by 11 % (R_net +7 %): not converged, as expected at +first order with a = h. The split, 24–27 % in the coil and 73–76 % in the shield, is the +measured one (25 % and 75 %, Zielke et al. 2021), and R_net is 44–48 % of the measured +0.6 Ω, which also has the matching network and the lines. 4 mm is the coarsest grid on which the 8 mm tube is one piece of metal. The converged value needs the real shield (about 2 mm thick, so h ≲ 1 mm there) and a resolved wire (h ≲ a / 4): about 10^8 edges on a uniform grid, a production run on mpiarray and cunumpy, or local refinement. From 116f5ef658a2cfcc2288483ff80814951b4c2a1e Mon Sep 17 00:00:00 2001 From: Max Date: Tue, 3 Nov 2026 21:29:09 +0100 Subject: [PATCH 3/3] Sort the statements of inductive3d. --- src/wellbalanced/spatial/inductive3d.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/wellbalanced/spatial/inductive3d.py b/src/wellbalanced/spatial/inductive3d.py index dba7db4..f583554 100644 --- a/src/wellbalanced/spatial/inductive3d.py +++ b/src/wellbalanced/spatial/inductive3d.py @@ -461,6 +461,12 @@ def _feed_edges(self, feed: Feed3D) -> tuple[np.ndarray, np.ndarray, list]: flat = [int(np.ravel_multi_index(tuple(n), grid.nodes)) for n in nodes] return np.asarray(edges), np.asarray(signs), flat + def _gradient_matrix(self) -> sp.csr_matrix: + grid = self.grid + if not hasattr(grid, "_gradient_matrix"): + grid._gradient_matrix = grid.gradient() # type: ignore[attr-defined] + return grid._gradient_matrix # type: ignore[attr-defined] + def _feed_source(self) -> np.ndarray: """``j`` of the feeds for 1 A [A]: the current along each path's edges. @@ -489,12 +495,6 @@ def _feed_source(self) -> np.ndarray: np.add.at(source, edges, signs) return source - def _gradient_matrix(self) -> sp.csr_matrix: - grid = self.grid - if not hasattr(grid, "_gradient_matrix"): - grid._gradient_matrix = grid.gradient() # type: ignore[attr-defined] - return grid._gradient_matrix # type: ignore[attr-defined] - def _edge_conductivity( self, background: PointField | None ) -> tuple[np.ndarray, dict[str, np.ndarray]]: