rtang/ANASEN_analysis
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

modified: anasen_fem/README.md added the Shockley-Ramo dot calculation capability adapted fom Sudarsan's branch with a new visusalisation scheme to show the +ve and-ve sides on a logscale

Browse Source 
modified:   anasen_fem/clean.sh
	new file:   anasen_fem/dotproduct.py
	new file:   anasen_fem/paraview_dotproduct.py
	modified:   anasen_fem/paraview_plotter.py
	modified:   anasen_fem/run.py
	modified:   anasen_fem/scalars.dat.names
	new file:   anasen_fem/scalars_weight.dat
	new file:   anasen_fem/scalars_weight.dat.names
	modified:   anasen_fem/wires2d.sif
	new file:   anasen_fem/wires2d_weight.sif
	modified:   anasen_fem/wires_gmsh2d_bc.py
This commit is contained in:
Vignesh Sitaraman 2026-09-23 12:14:05 -04:00
parent aebb9854ec
commit 3edbfa5bd8
12 changed files with 920 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
anasen_fem/README.md
View File

@ -23,5 +23,24 @@
- 2d simulations of fields only. gmsh for meshing, elmer for fem, paraview to plot - 2d simulations of fields only. gmsh for meshing, elmer for fem, paraview to plot
- Before running, open `paraview_plotter.py` to make the bash shebang (#!) point to the location of `pvpython` or `pvbatch` - Before running, open `paraview_plotter.py` to make the bash shebang (#!) point to the location of `pvpython` or `pvbatch`
- `python3 run.py` should run everything in order, and is hopefully all the files are self-documenting - `python3 run.py` should run everything in order, and is hopefully all the files are self-documenting
* v0.0.2, planned TODO * v0.0.2, September 2026
- Adds the Ramo weighting field and the `E . E_w` dot product for one cathode wire.
- `wires_gmsh2d_bc.py` now cuts the wires into the gas disk with `occ.fragment` instead of `mesh.embed`, so each wire boundary is a real edge of the domain. One cathode (optional 2nd argument, default 1) gets its own physical group, **tag 40**; the other 23 stay in tag 30.
- Two solves, same mesh, identical apart from the potentials:
| file | tag 40 | output |
|---|---|---|
| `wires2d.sif` | 0 V, like every other cathode | `wires2d/elfield_anasen_t0001.vtu` |
| `wires2d_weight.sif` | 1 V, all else 0 V | `wires2d/elfield_weight_t0001.vtu` |
- `dotproduct.py` reads both and writes `wires2d/dotproduct.vtu` with `DotProduct = E . E_w` (V/m^2). It aborts if the node sets differ or if either field is all zeros — both are silent failures otherwise.
- `paraview_dotproduct.py` renders it. `DOT_MIN`/`DOT_MAX` set the colour range and decide whether the picture shows anything.
- `paraview_plotter.py` gains two views: the equipotentials over one quadrant, and streamlines over the same quadrant. Streamlines are drift paths up to diffusion, so they show which wire collects charge from where; `SEED_POTENTIAL` and `SEED_STRIDE` control them.
- Six PNGs per z-locus, archived as before (`<Stem>_z_<count>_<z>[_quarter].png`). `Field_ouput` keeps its typo so the new files sort with the dozen already in `png/`.
- `STAGE` in `run.py`: `1` = mesh + physical solve, `2` = weighting solve + dot product reusing stage 1's mesh, `0` = both. Stage 2 never re-meshes, so it applies to whichever z stage 1 ran last.
- **Do not pass `-autoclean` to ElmerGrid.** It renumbers the physical groups (13/10/20/30/40 become 1/4/5/6/7), so every `Target Bodies`/`Target Boundaries` in the sifs matches nothing, Elmer solves nothing, and the potential comes out identically zero with no error.
- `mesh.recombine()` and `mesh.refine()` are both off: recombination ran over half an hour on the barrel surface without finishing, and refine took the mesh to ~16M nodes.
- Note on the solver stack: both sifs are the same file bar the potentials, so whatever `FluxSolver` does to `Electric Field` it does equally to both and the dot product stays consistent.
* v0.0.3, planned TODO
- Garfield to take Elmer results and perform charge-transport - Garfield to take Elmer results and perform charge-transport
- Sweep the weighting solve over all 24 cathodes rather than one at a time

4
anasen_fem/clean.sh
View File

@ -1,2 +1,2 @@
rm wires2d.msh rm -f wires2d.msh
rm wires2d/* rm -rf wires2d/*

148
anasen_fem/dotproduct.py Normal file
View File

@ -0,0 +1,148 @@
#!/usr/bin/env python3
"""Node-by-node dot product of the physical field with a Ramo weighting field.
Reads the two Elmer VTUs from wires2d.sif and wires2d_weight.sif and writes one
carrying DotProduct = E . E_w. E is in V/m and E_w in 1/m (its electrode is
driven to 1 V), so DotProduct is in V/m^2.
By Shockley-Ramo a charge q drifting at v induces i = -q v . E_w, and for
v = mu E that is -q mu (E . E_w). So this map is the induced-current density up
to the constant -q*mu; --scale folds it in. The signs of q and mu differ
between the electron and ion components, so they are left to the caller.
python3 dotproduct.py
python3 dotproduct.py --scale -1.0e-4
"""
import argparse
import sys
import numpy as np
try:
import meshio
except ImportError:
sys.exit("meshio is required: pip install meshio")
def find_array(point_data, wanted, label):
"""Look up a point-data array, ignoring case and spacing.
Elmer lowercases variable names into the VTU and the spelling has varied
between versions, so match loosely and list what is there on failure.
"""
def norm(s):
return "".join(s.lower().split())
target = norm(wanted)
for key in point_data:
if norm(key) == target:
return point_data[key]
sys.exit(
"%s: no point-data array matching %r.\nAvailable arrays: %s"
% (label, wanted, sorted(point_data))
)
def as_vectors(arr, label):
"""Return an (N, ncomp) float array, tolerating flat or 2-component data."""
a = np.asarray(arr, dtype=float)
if a.ndim == 1:
sys.exit("%s: expected a vector field, got a scalar array" % label)
if a.ndim != 2:
sys.exit("%s: expected a 2D array, got shape %s" % (label, a.shape))
return a
def main():
ap = argparse.ArgumentParser(description=__doc__,
formatter_class=argparse.RawDescriptionHelpFormatter)
ap.add_argument("--phys", default="wires2d/elfield_anasen_t0001.vtu",
help="VTU from the physical solve (wires2d.sif)")
ap.add_argument("--weight", default="wires2d/elfield_weight_t0001.vtu",
help="VTU from the weighting solve (wires2d_weight.sif)")
ap.add_argument("--out", default="wires2d/dotproduct.vtu",
help="output VTU")
ap.add_argument("--field", default="electric field",
help="name of the electric-field array in both inputs")
ap.add_argument("--scale", type=float, default=None,
help="constant (e.g. -q*mu) multiplying the dot product; "
"writes an extra ScaledDotProduct array")
ap.add_argument("--rtol", type=float, default=0.0,
help="relative tolerance for the node-position check")
ap.add_argument("--atol", type=float, default=1e-12,
help="absolute tolerance (m) for the node-position check")
args = ap.parse_args()
phys = meshio.read(args.phys)
weight = meshio.read(args.weight)
# Taken node by node, so both solves must have produced the same node list.
# If they did not, every value below is meaningless and nothing else checks.
if phys.points.shape != weight.points.shape:
sys.exit(
"Node counts differ: %s has %d nodes, %s has %d. Both solves must "
"read the same mesh (same wires2d/ directory)."
% (args.phys, len(phys.points), args.weight, len(weight.points))
)
if not np.allclose(phys.points, weight.points, rtol=args.rtol, atol=args.atol):
worst = np.max(np.linalg.norm(phys.points - weight.points, axis=1))
sys.exit(
"Node positions do not match (largest discrepancy %.3e m). The two "
"solves were run on different meshes." % worst
)
E = as_vectors(find_array(phys.point_data, args.field, args.phys), args.phys)
Ew = as_vectors(find_array(weight.point_data, args.field, args.weight), args.weight)
# A zero field means the solve never ran: Laplace with no Dirichlet
# condition applied has the trivial solution and Elmer reports no error.
# Without this the chain completes and writes zeros, and the only symptom
# is a blank plot.
for label, arr, path in (
("physical", E, args.phys),
("weighting", Ew, args.weight),
):
if not np.any(arr):
sys.exit(
"the %s field in %s is identically zero.\n"
"Elmer returned the trivial solution, which means no boundary "
"condition was applied. Compare the Target Bodies and Target "
"Boundaries numbers in the .sif against wires2d/mesh.names -- "
"ElmerGrid's -autoclean renumbers them." % (label, path)
)
# Elmer may pad 2D vectors to three components; compare only the
# components both arrays actually carry.
ncomp = min(E.shape[1], Ew.shape[1])
if E.shape[1] != Ew.shape[1]:
print("note: component counts differ (%d vs %d), using the first %d"
% (E.shape[1], Ew.shape[1], ncomp))
dot = np.einsum("ij,ij->i", E[:, :ncomp], Ew[:, :ncomp])
out = meshio.Mesh(
points=phys.points,
cells=phys.cells,
point_data=dict(phys.point_data),
cell_data=dict(phys.cell_data),
)
out.point_data["DotProduct"] = dot
out.point_data["WeightingField"] = Ew
if args.scale is not None:
out.point_data["ScaledDotProduct"] = args.scale * dot
meshio.write(args.out, out)
finite = dot[np.isfinite(dot)]
print("nodes : %d" % len(dot))
print("components used : %d" % ncomp)
print("DotProduct min : %.6e V/m^2" % finite.min())
print("DotProduct max : %.6e V/m^2" % finite.max())
print("DotProduct mean : %.6e V/m^2" % finite.mean())
if len(finite) != len(dot):
print("WARNING: %d non-finite nodes" % (len(dot) - len(finite)))
print("wrote %s" % args.out)
if __name__ == "__main__":
main()

279
anasen_fem/paraview_dotproduct.py Executable file
View File

@ -0,0 +1,279 @@
#!/home/vsitaraman/ParaView-6.1.0-RC1-MPI-Linux-Python3.12-x86_64/bin/pvbatch
# Point the shebang above at your own pvbatch/pvpython, as for paraview_plotter.py.
#
# Renders the E . E_w map written by dotproduct.py.
#
# By Shockley-Ramo the dot product describes the current induced on ONE
# electrode: the cathode that wires2d_weight.sif held at 1 V. Everywhere else
# the number is just the tail of that electrode's weighting field, so the view
# is centred on that wire and nothing else. Its position is read out of the
# weighting solve rather than recomputed from the geometry -- the driven
# electrode is wherever the weighting potential reaches 1.
#
# Usage (from anasen_fem/):
# ./paraview_dotproduct.py [dotproduct.vtu] [output.png] [weight.vtu]
#
# Re-plotting an archived z needs its archived weighting solve as well, since
# that is where the driven cathode's position is read from:
# ./paraview_dotproduct.py wires2d/vtu_files/dotproduct_10_0.0000.vtu \
# png/DotProduct_output_z_10_0.0000.png \
# wires2d/vtu_files/elfield_weight_10_0.0000.vtu
import sys
import numpy as np
from paraview.simple import *
INPUT_VTU = sys.argv[1] if len(sys.argv) > 1 else "wires2d/dotproduct.vtu"
OUTPUT_IMAGE = sys.argv[2] if len(sys.argv) > 2 else "DotProduct_output.png"
# Weighting solve, used only to locate the driven cathode. Must be the one that
# belongs to this dot product -- the wire moves with z.
WEIGHT_VTU = sys.argv[3] if len(sys.argv) > 3 else "wires2d/elfield_weight_t0001.vtu"
# Half-height of the view in metres. Cathodes sit every 15 degrees at a radius
# of about 36.4 mm, so the neighbours either side are 9.5 mm away; 14 mm takes
# them both in with a margin.
# Image size. The structure round the driven cathode is much wider than it is
# tall, so a square frame spends most of its area on empty gas.
VIEW_SIZE = [2000, 1200]
# CameraParallelScale, which is the half-HEIGHT of the view in metres. The
# width follows from the aspect: half-height * (VIEW_SIZE[0] / VIEW_SIZE[1]).
# So widening the frame alone does not crop anything, it just shows more
# sideways; the scale has to come down with it. 0.0084 at 2000x1200 keeps the
# same 28 mm width as the old square 0.014 and trims the height to 16.8 mm.
ZOOM_SCALE = 0.0084
# The driven cathode swings round the axis as z changes, so with a fixed "up"
# every z comes out rotated by a different amount and the shots cannot be laid
# side by side. Aligning up with the outward radial direction through the
# cathode puts the detector axis below and the barrel wall above in every
# frame, whatever z. Set False for plain screen-aligned axes.
ALIGN_RADIAL = True
# Shift of the view centre along the outward radial, in metres. Up in the frame
# is the outward radial, so a NEGATIVE value walks the camera in towards the
# detector axis and the cathode rides higher in the picture. Positive drops it.
# Purely framing; it does not move the wire.
CENTER_SHIFT = -0.0036
# Second, wider shot over the same quadrant paraview_plotter.py frames, so the
# dot product can be laid beside the field and contour views of that z. Centre
# is at x < 0, y > 0; negate the y for the lower-left quadrant instead.
# Colour range for DotProduct. Both ends negative: E and E_w point in broadly
# opposite senses over most of the drift volume.
# Most of the painted patch lies between -2e7 and -2e6. Against a range
# reaching -3.1e8 that is the top few percent of the scale, so it all came out
# at the white end of Greens. Narrowing the range spreads those values over the
# green; anything stronger simply clamps to the darkest end, which only affects
# the last fraction of a millimetre at the wire surface.
DOT_MIN, DOT_MAX = -8.0e7, -3.0e6
# Opacity of the painted map. Anything at least as strong as OPAQUE_BELOW is
# drawn solid; from there it ramps down to transparent at FADE_TO, so the grey
# shows through with a soft round edge rather than the ragged one a hard cut
# leaves.
#
# The plateau matters more than either endpoint. Ramping all the way from
# DOT_MIN (-3.1e8) meant the opacity at -5e7 was only about 0.13, so the whole
# patch read as washed out however the fade point was set.
OPAQUE_BELOW = -2.0e7
FADE_TO = -2.0e6
# Coarse cut, only to keep the far field out of the render. Must be nearer zero
# than FADE_TO so its hard edge lands where the surface is already fully
# transparent and cannot be seen.
PAINT_CUT = -5.0e5
# "signed-log" maps sign and magnitude together through a signed logarithm and
# a diverging colour map; "linear" is the earlier one-sided Greens scale, which
# the constants above still drive.
COLOR_MODE = "signed-log"
# X0 sets where the mapping turns from linear to logarithmic, and so how much
# of the weak far field survives at all: at 1e6 the data map to -2.51 .. 1.62,
# against -3.51 .. 2.61 at 1e5. Raise it to paint less.
SYMLOG_X0 = 1.0e6
SYMLOG_LIMIT = 2.5 # symmetric colour limit; match it to the mapped range
SYMLOG_OPAQUE = 0.6 # |value| at which the map is fully opaque
SYMLOG_FADE = 0.3 # |value| below which it is transparent
# Shape of the two opacity ramps, as opposed to where they sit.
#
# A piecewise function's midpoint describes the segment running RIGHTWARD from
# each node, and the two ramps run in opposite senses -- falling on the
# negative side, rising on the positive. So a single midpoint below 0.5 applied
# to both makes the negative lobe fade out sooner while the positive lobe goes
# opaque sooner, rendering the two polarities at different effective
# thresholds. Mirroring it keeps them symmetric: below 0.5 both sides go opaque
# sooner, so more is painted, equally on each.
#
# Sharpness is the independent dial for how abrupt the edge is: 0 smooth, 1 a
# hard step. 0.5 is already close to a step.
SYMLOG_MIDPOINT = 0.3
SYMLOG_SHARPNESS = 0.0
DIVERGING_PRESET = "PRGn" # purple / white / green; zero sits at white
BACKGROUND = [0.35, 0.34, 0.33] # grey
# Equipotential lines over the whole frame, painted gas or not. Off by default:
# the map is the subject here, and the lines are already in the contour plots.
SHOW_CONTOURS = False
CONTOUR_BY = "potential"
CONTOUR_LEVELS = list(np.arange(0, 660, 660 / 40.0))
def threshold_between(source, array, low, high):
"""Threshold filter, tolerating the pre-5.10 property names."""
t = Threshold(Input=source, Scalars=["POINTS", array])
try:
t.LowerThreshold = low
t.UpperThreshold = high
t.ThresholdMethod = "Between"
except Exception:
t.ThresholdRange = [low, high]
return t
# --- locate the driven cathode: the weighting potential is 1 on it ---
weight_reader = XMLUnstructuredGridReader(FileName=[WEIGHT_VTU])
driven = threshold_between(weight_reader, "potential", 0.99, 1.01)
driven.UpdatePipeline()
info = driven.GetDataInformation()
if info.GetNumberOfPoints() == 0:
sys.exit(
"no nodes with weighting potential ~1 in %s -- the weighting solve did\n"
"not drive any electrode, so there is nothing to centre on." % WEIGHT_VTU
)
b = info.GetBounds()
CATHODE = [0.5 * (b[0] + b[1]), 0.5 * (b[2] + b[3]), 0.0]
print("driven cathode at x = %.5f m, y = %.5f m" % (CATHODE[0], CATHODE[1]))
radius = (CATHODE[0] ** 2 + CATHODE[1] ** 2) ** 0.5
outward = [CATHODE[0] / radius, CATHODE[1] / radius, 0.0] if radius > 0 else [0.0, 1.0, 0.0]
# Frame centre, offset from the wire along that radial.
CENTER = [
CATHODE[0] + CENTER_SHIFT * outward[0],
CATHODE[1] + CENTER_SHIFT * outward[1],
0.0,
]
Delete(driven)
Delete(weight_reader)
# --- the dot product map ---
reader = XMLUnstructuredGridReader(FileName=[INPUT_VTU])
renderView = GetActiveViewOrCreate("RenderView")
renderView.ViewSize = VIEW_SIZE
renderView.OrientationAxesVisibility = 1
renderView.UseColorPaletteForBackground = 0
renderView.Background = BACKGROUND
renderView.MultiSamples = 8
if COLOR_MODE == "signed-log":
# E . E_w spans about 3.5 decades negative and 2.6 positive, and 38% of the
# nodes are positive -- charge drifting there induces current of the
# opposite polarity on the cathode. A linear scale over one sign shows
# neither the decades nor the sign, so map through a signed log:
#
# sign(d) * log10(1 + |d| / SYMLOG_X0)
#
# sign() is not in the calculator namespace, hence d/(|d|+1); the +1 is
# negligible beside the 1e5 scale and the expression tends to 0 anyway
# where d does.
calc = PythonCalculator(Input=reader)
calc.ArrayName = "SignedLogDot"
calc.ArrayAssociation = "Point Data"
calc.Expression = (
"(DotProduct/(abs(DotProduct)+1.0))*log10(1.0+abs(DotProduct)/%g)"
% SYMLOG_X0
)
painted = calc
array = "SignedLogDot"
colour_range = (-SYMLOG_LIMIT, SYMLOG_LIMIT) # symmetric, so zero is neutral
preset = DIVERGING_PRESET
# Transparent either side of zero, opaque past SYMLOG_OPAQUE, both signs.
_m, _s = SYMLOG_MIDPOINT, SYMLOG_SHARPNESS
opacity_points = [
-SYMLOG_LIMIT, 1.0, 0.5, _s,
-SYMLOG_OPAQUE, 1.0, 1.0 - _m, _s, # mirrored: this ramp falls
-SYMLOG_FADE, 0.0, 0.5, _s,
SYMLOG_FADE, 0.0, _m, _s, # this one rises
SYMLOG_OPAQUE, 1.0, 0.5, _s,
SYMLOG_LIMIT, 1.0, 0.5, _s,
]
bar_title = "sign(E.Ew) log10(1+|E.Ew|/%g)" % SYMLOG_X0
else:
painted = threshold_between(reader, "DotProduct", -1.0e30, PAINT_CUT)
array = "DotProduct"
colour_range = (DOT_MIN, DOT_MAX)
preset = "Greens"
opacity_points = [
DOT_MIN, 1.0, 0.5, 0.0,
OPAQUE_BELOW, 1.0, 0.5, 0.0,
FADE_TO, 0.0, 0.5, 0.0,
]
bar_title = "E . E_w (V/m^2)"
surface_display = Show(painted, renderView)
surface_display.Representation = "Surface"
ColorBy(surface_display, ("POINTS", array))
dotLUT = GetColorTransferFunction(array)
dotLUT.ApplyPreset(preset, True)
dotLUT.RescaleTransferFunction(colour_range[0], colour_range[1])
# Points go (value, opacity, midpoint, sharpness) in ascending value.
dotPWF = GetOpacityTransferFunction(array)
dotPWF.Points = opacity_points
dotLUT.ScalarOpacityFunction = dotPWF
dotLUT.EnableOpacityMapping = 1
surface_display.SetScalarBarVisibility(renderView, True)
try:
GetScalarBar(dotLUT, renderView).Title = bar_title
except Exception as exc:
print("note: could not set the scalar bar title (%s)" % exc)
# --- contour lines over the whole field, including the unpainted gas ---
if SHOW_CONTOURS:
contour_filter = Contour(Input=reader, ContourBy=["POINTS", CONTOUR_BY])
contour_filter.Isosurfaces = CONTOUR_LEVELS
contour_display = Show(contour_filter, renderView)
contour_display.LineWidth = 1.0
contour_display.DiffuseColor = [1.0, 1.0, 1.0]
# ColorBy(rep, None) is rejected by ParaView 6.1 -- it forwards the
# association as the string 'NONE' and raises -- so set the property.
contour_display.ColorArrayName = [None, ""]
contour_display.SetScalarBarVisibility(renderView, False)
# --- camera: parallel projection, centred on the driven cathode ---
# ParaView resets the camera on a view's FIRST render, throwing away anything
# set beforehand, so burn that render here and set the camera after it.
# Measured: scale 0.008 set before the first Render comes back as 0.021; set
# again afterwards it holds.
Render()
view = GetActiveView()
view.CameraParallelScale = ZOOM_SCALE
view.CenterOfRotation = CENTER
view.CameraPosition = [CENTER[0], CENTER[1], 1.0]
view.CameraFocalPoint = CENTER
if ALIGN_RADIAL:
up = outward
print("view rotated %.1f deg so the outward radial is up"
% np.degrees(np.arctan2(up[0], up[1])))
else:
up = [0.0, 1.0, 0.0]
view.CameraViewUp = up
view.CameraParallelProjection = 1
Render()
SaveScreenshot(OUTPUT_IMAGE)
print("Saved %s" % OUTPUT_IMAGE)

170
anasen_fem/paraview_plotter.py
View File

@ -4,10 +4,55 @@ import numpy as np
import sys import sys
from paraview.simple import * from paraview.simple import *
reader = XMLUnstructuredGridReader(FileName=["wires2d/elfield_anasen_t0001.vtu"]) # Optional argument so an archived VTU can be re-plotted without re-solving:
# ./paraview_plotter.py wires2d/vtu_files/elfield_anasen_10_0.0000.vtu
INPUT_VTU = sys.argv[1] if len(sys.argv) > 1 else "wires2d/elfield_anasen_t0001.vtu"
# --- contour settings ---
# The whole-PC view keeps the original sparse, evenly spaced levels; at that
# scale more lines just fill in solid. The quadrant gets the denser, biased
# set, since that is where the structure between wires is actually readable.
POT_MIN, POT_MAX = 0.0, 660.0
N_CONTOURS_FULL, CONTOUR_BIAS_FULL = 40, 1.0
N_CONTOURS_QUARTER, CONTOUR_BIAS_QUARTER = 80, 1.0
# Line widths. Thin lines wash out their own colour, so the quadrant views draw
# heavier than the full-PC one; tubes keep them smooth at 2000 px.
LINE_WIDTH_FULL = 3.0
LINE_WIDTH_QUARTER = 2.5
STREAM_LINE_WIDTH = 2.5
# Arrows are glyphed onto the contour points, so doubling the quadrant contour
# density doubled the arrow count too. This offsets that.
GLYPH_STRIDE = 96
# --- quadrant view ---
# Centre sits at radius 35.4 mm, near the outer edge of the wire band, which
# pushes the wires to one side and leaves the far corner empty. Pulled in to
# 31.8 mm so the band sits more centrally.
ZOOM_CENTER = [-0.0225, 0.0225, 0.0]
ZOOM_SCALE = 0.022 # half-height of the view in metres; larger = more in frame
def contour_levels(n, bias):
"""Isosurface values.
phi ~ ln(r) around a wire, so evenly spaced levels crowd onto the wire
surfaces and thin out in the gas. Raising a symmetric parameter to `bias`
pulls levels towards mid-range, where the saddles between wires are,
without adding more of them. 1.0 is plain linear spacing; much above 1.5
and the middle levels start landing on top of each other.
"""
u = np.linspace(-1.0, 1.0, n)
u = np.sign(u) * np.abs(u) ** bias
mid = 0.5 * (POT_MIN + POT_MAX)
return list(mid + 0.5 * (POT_MAX - POT_MIN) * u)
reader = XMLUnstructuredGridReader(FileName=[INPUT_VTU])
contour_filter = Contour(Input=reader,ContourBy = 'potential') contour_filter = Contour(Input=reader,ContourBy = 'potential')
contour_filter.Isosurfaces = [i for i in np.arange(0,660,660/40.)] contour_filter.Isosurfaces = contour_levels(N_CONTOURS_FULL, CONTOUR_BIAS_FULL)
renderView = GetActiveViewOrCreate('RenderView') renderView = GetActiveViewOrCreate('RenderView')
renderView.ViewSize = [2000,2000] renderView.ViewSize = [2000,2000]
@ -20,7 +65,7 @@ renderView.MultiSamples = 8 # 0 disables it, 4-8 is usually sufficient
ResetCamera() ResetCamera()
contour_display = Show(contour_filter, renderView) contour_display = Show(contour_filter, renderView)
contour_display.LineWidth = 3.0 # Increase this for thicker lines contour_display.LineWidth = LINE_WIDTH_FULL
contour_display.RenderLinesAsTubes = 1 # Makes lines look smoother at high res contour_display.RenderLinesAsTubes = 1 # Makes lines look smoother at high res
#colorbar #colorbar
contour_display_potentialLUT = GetColorTransferFunction('potential', contour_display, separate=True) contour_display_potentialLUT = GetColorTransferFunction('potential', contour_display, separate=True)
@ -62,23 +107,30 @@ Render()
SaveScreenshot("contour_output.png") SaveScreenshot("contour_output.png")
#make glyps for field lines # Everything from here on is a quadrant view: denser levels, heavier lines.
contour_display.LineWidth = 1.0 # Increase this for thicker lines contour_filter.Isosurfaces = contour_levels(N_CONTOURS_QUARTER, CONTOUR_BIAS_QUARTER)
contour_display.RenderLinesAsTubes = 0 # Makes lines look smoother at high res contour_display.LineWidth = LINE_WIDTH_QUARTER
contour_display.RenderLinesAsTubes = 1
# 1. Get the active view # 1. Get the active view
view = GetActiveView() view = GetActiveView()
# 1. Set the Focal Point to the middle of the quadrant in metres # 1. Set the Focal Point to the middle of the quadrant in metres
zoom_center = [-0.025, 0.025, 0.0] zoom_center = ZOOM_CENTER
# 2. Tighten the Parallel Scale # 2. Tighten the Parallel Scale
view.CameraParallelScale = 0.015 view.CameraParallelScale = ZOOM_SCALE
# 3. Position the Camera (0.5m away is fine) # 3. Position the Camera (0.5m away is fine)
view.CameraPosition = [zoom_center[0], zoom_center[1], 0.5] view.CameraPosition = [zoom_center[0], zoom_center[1], 0.5]
view.CameraFocalPoint = zoom_center view.CameraFocalPoint = zoom_center
view.CameraViewUp = [0.0, 1.0, 0.0] view.CameraViewUp = [0.0, 1.0, 0.0]
# Equipotentials over one quadrant, saved before the arrows go on so this and
# Field_output.png below share a camera and overlay exactly.
Render()
SaveScreenshot("contour_quarter_output.png")
# pot_threshold = Threshold(Input=reader) # pot_threshold = Threshold(Input=reader)
# pot_threshold.Scalars = ['POINTS', 'potential'] # pot_threshold.Scalars = ['POINTS', 'potential']
# pot_threshold.ThresholdMethod = 'Above Upper Threshold' # pot_threshold.ThresholdMethod = 'Above Upper Threshold'
@ -95,7 +147,7 @@ glyph.ScaleArray = ['POINTS', 'No scale array']
glyph.ScaleFactor = 0.001 glyph.ScaleFactor = 0.001
glyph.GlyphMode = 'Every Nth Point' glyph.GlyphMode = 'Every Nth Point'
glyph.Stride = 32 glyph.Stride = GLYPH_STRIDE
# --- 3. Display the Glyphs --- # --- 3. Display the Glyphs ---
glyph_display = Show(glyph, renderView) glyph_display = Show(glyph, renderView)
@ -107,7 +159,7 @@ glyph_display.Representation = 'Surface'
# This is the critical line: Color the arrows by the 'potential' scalar # This is the critical line: Color the arrows by the 'potential' scalar
ColorBy(glyph_display, ('POINTS', 'potential')) ColorBy(glyph_display, ('POINTS', 'potential'))
glyph_display.LookupTable = contour_display_potentialLUT glyph_display.LookupTable = contour_display_potentialLUT
contour_display_potentialLUT.RescaleTransferFunction(0.0, 660.0) contour_display_potentialLUT.RescaleTransferFunction(POT_MIN, POT_MAX)
# Optional: Disable the scalar bar for the arrows to avoid cluttering # Optional: Disable the scalar bar for the arrows to avoid cluttering
# the existing 'potential' scalar bar. # the existing 'potential' scalar bar.
@ -116,3 +168,101 @@ glyph_display.SetScalarBarVisibility(renderView, False)
# --- 4. Final Render --- # --- 4. Final Render ---
Render() Render()
SaveScreenshot("Field_output.png") SaveScreenshot("Field_output.png")
# --- Streamlines over the same quadrant ---
# Field lines rather than arrows: with no B field the drift velocity is
# parallel to E, so these are drift paths up to diffusion, i.e. a map of which
# wire collects charge from where. The glyphs above are fixed length, so they
# show direction only. Quadrant only; over the whole PC it is an unreadable mat.
# Equipotentials the streamlines start from. One level near the anodes seeds
# lines only in a ring around each anode, so with the lines truncated short the
# rest of the frame comes out empty. Spanning the range instead puts seeds
# right across the gas: every equipotential is a closed curve winding between
# the wires, so together they cover the volume.
SEED_POTENTIALS = [660.0, 600.0, 550.0, 500.0, 450.0, 400.0,350.0, 300.0,250.0, 250.0, 200.0, 150.0, 100.0, 50.0]
# Keep 1 seed point in N along the seed contours. Measured point counts per
# level, for working out the total:
#
# 600 V 4970 300 V 29971
# 500 V 8861 200 V 51558
# 400 V 16331 100 V 51645
# 50 V 27259
#
# 250 V and 150 V are not measured; from the trend they are of order 40000 and
# 50000, putting the nine levels near 280000 points. Stride 3 is then ~94000
# seeds, of which roughly an eighth land in the rendered quadrant. That is
# three times the previous seed count, but each line is now 4 mm instead of
# 20 mm, so the integration work is comparable.
SEED_STRIDE = 3
# Arc length each line is allowed to run, in metres. Field lines terminate on
# electrodes, so any line long enough to reach one will converge with every
# other line reaching the same wire -- that convergence is the field, not a
# plotting artefact, and the only way to avoid drawing it is to stop the lines
# short. The anode-cathode gap runs 4.3-6.2 mm depending on z, so 3 mm keeps
# every line clear of the cathodes. Raise it towards 0.005 for longer lines,
# and they will start meeting again at the wires.
# Short lines, and the gaps between them closed by seeding more planes rather
# than by making each line longer. 4 mm is just under the 4.3 mm closest
# anode-cathode approach, so no line quite reaches a wire and they stay
# separate instead of converging.
STREAM_MAX_LENGTH = 0.0005
# StreamTracer counts steps as well as length, and its step sizes are in
# cell-length units rather than metres: 0.2 of a cell by default, and cells are
# 0.05 mm at the wires. 2000 steps therefore buys about 20 mm at best and less
# wherever the mesh is fine, so the default caps the lines before
# STREAM_MAX_LENGTH ever does. Raise it well past what the length needs.
MAX_STEPS = 20000
STREAM_DIRECTION = 'BOTH' # 'BOTH' also draws the half that runs into the anode
STREAM_RANGE = [-200.0, 650.0] # full electrode span, needle to anode
# Clear the quadrant view of everything else so only the field lines show.
Hide(glyph, renderView)
Hide(contour_filter, renderView)
seed_contour = Contour(Input=reader, ContourBy=['POINTS', 'potential'])
seed_contour.Isosurfaces = SEED_POTENTIALS
# StreamTracer's own Point Cloud and Line seeds are 3D and this mesh is a
# single plane at z = 0, so they would mostly land outside it. Seeding from a
# contour of the mesh keeps every seed in the plane by construction.
seed_points = MaskPoints(Input=seed_contour)
seed_points.OnRatio = SEED_STRIDE
seed_points.RandomSampling = 0
stream = StreamTracerWithCustomSource(Input=reader, SeedSource=seed_points)
stream.Vectors = ['POINTS', 'electric field']
# FORWARD follows E, which runs from the anode outwards, so each line leaves
# its seed and heads for a cathode. With BOTH, the backward half of every line
# also climbs to the anode the seed ring encircles, and all of them meet there
# -- that is the joining. BOTH is still the honest full field line if you want
# it; this only changes what gets drawn.
stream.IntegrationDirection = STREAM_DIRECTION
stream.MaximumStreamlineLength = STREAM_MAX_LENGTH
stream.MaximumSteps = MAX_STEPS
# Keep integration in the z = 0 plane. Newer ParaView only, so tolerate it
# being absent.
try:
stream.SurfaceStreamlines = 1
except Exception as exc:
print("note: SurfaceStreamlines unavailable (%s)" % exc)
stream_display = Show(stream, renderView)
stream_display.Representation = 'Surface'
stream_display.LineWidth = STREAM_LINE_WIDTH
stream_display.RenderLinesAsTubes = 1
ColorBy(stream_display, ('POINTS', 'potential'))
streamLUT = GetColorTransferFunction('potential', stream_display, separate=True)
streamLUT.ApplyPreset('Cool to Warm', True)
streamLUT.RescaleTransferFunction(STREAM_RANGE[0], STREAM_RANGE[1])
stream_display.LookupTable = streamLUT
stream_display.SetScalarBarVisibility(renderView, True)
renderView.Background = [0.0, 0.0, 0.0]
Render()
SaveScreenshot("Streamlines_quarter_output.png")

119
anasen_fem/run.py
View File

@ -1,24 +1,109 @@
import code
import os import os
import sys
# val=-174.3 # Cathode wire (0..23) that wires2d_weight.sif drives to 1 V. Passed to the
val=0 # mesher, which puts that one wire in its own physical group, tag 40.
count=10 SELECTED_CATHODE = 1
while val<174.3+0.1:
# 0 = everything in one pass
# 1 = mesh and the physical field only
# 2 = weighting field and dot product only, reusing stage 1's mesh
STAGE = 0
# gmsh comes from apt and lives in the system python3; meshio is pip-installed
# into the venv, which cannot see apt packages. Neither interpreter has both,
# so name them explicitly rather than relying on whichever python3 is on PATH.
MESH_PY = "/usr/bin/python3"
DOT_PY = os.path.expanduser("~/myenv/bin/python3")
def run(cmd):
"""os.system that stops on failure.
Worth the four lines: a solve that quietly did nothing once produced a
field of zeros, and every later step ran happily on it.
"""
print(cmd)
if os.system(cmd) != 0:
sys.exit("failed: " + cmd)
# z-loci to run, in mm.
#
# The wires twist: over z = 0 -> 174.3 the anode lattice turns 60 deg against
# the cathode lattice, which is exactly 4 cathode pitches of 15 deg. So the
# field pattern goes through the same cycle four times across that range, while
# the wire radii drift slowly and monotonically (anodes 32.0 -> 37.0 mm).
#
# A uniform 17.43 step samples the fast cycle 2.75 times per period, which
# aliases it -- consecutive slices land at relative azimuths 7.5, 14.1, 5.7,
# 12.2, 3.5 deg, jumping around rather than sweeping. Sample the two variations
# separately instead:
# - 8 slices across one pitch, which walks the azimuth in even ~2 deg steps
# at essentially fixed radius
# - 5 slices one whole pitch apart, which holds the azimuth near-fixed and
# walks the radius over its full range
# 12 z-loci in total, against 11 for the uniform sweep.
PITCH_Z = 174.3 / 4.0
Z_VALUES = sorted(set(
[round(i * PITCH_Z / 8.0, 4) for i in range(8)]
+ [round(i * PITCH_Z, 4) for i in range(5)]
))
# Running only two loci for now; delete this line for the full sweep.
# 19.794 is the crossover: the relative azimuth wraps through zero there, so an
# anode sits directly radially outside a cathode, and the closest anode-cathode
# approach over the whole range (4.339 mm) is at 19.707. It is deliberately off
# the grid above, whose nearest point is 21.7875.
Z_VALUES = [0.0, 19.794]
# Stage 2 reuses whatever single mesh is sitting in wires2d/, so it cannot span
# several z. Without this it would happily solve the weighting field on one z's
# mesh and archive the result under another z's name -- and dotproduct.py could
# not catch it, since the node sets would agree perfectly.
if STAGE == 2 and len(Z_VALUES) > 1:
sys.exit(
"STAGE 2 applies to the one mesh in wires2d/, which is whichever z\n"
"stage 1 ran last. Set Z_VALUES to that single z."
)
for count, val in enumerate(Z_VALUES, start=10):
print(val) print(val)
os.system("python3 wires_gmsh2d_bc.py "+str(val)) os.system("mkdir -p wires2d/mesh_files wires2d/sif_files wires2d/vtu_files png")
os.system("ElmerGrid 14 2 wires2d.msh -2d")
os.system("ElmerSolver wires2d.sif") # Each stage archives its own output before the next one starts, so a
os.system("./paraview_plotter.py") # failure in stage 2 does not throw away what stage 1 already produced.
if STAGE in (0, 1):
run("%s wires_gmsh2d_bc.py %s %d" % (MESH_PY, val, SELECTED_CATHODE))
run("ElmerGrid 14 2 wires2d.msh -2d")
run("ElmerSolver wires2d.sif")
run("./paraview_plotter.py")
os.system("cp wires2d.msh wires2d/mesh_files/wires2d%02d_%1.4f.msh"%(count,val))
os.system("cp wires2d.sif wires2d/sif_files/wires2d_%02d_%1.4f.sif"%(count,val))
os.system("cp wires2d/elfield_anasen_t0001.vtu wires2d/vtu_files/elfield_anasen_%02d_%1.4f.vtu"%(count,val))
os.system("cp contour_output.png png/Contour_output_z_%02d_%1.4f.png"%(count,val))
os.system("cp contour_quarter_output.png png/Contour_output_z_%02d_%1.4f_quarter.png"%(count,val))
os.system("cp Field_output.png png/Field_ouput_z_%02d_%1.4f.png"%(count,val))
os.system("cp Streamlines_quarter_output.png png/Streamlines_output_z_%02d_%1.4f_quarter.png"%(count,val))
# Stage 2 does not re-mesh: it reuses wires2d/ as stage 1 left it, which is
# what makes the node-by-node dot product valid.
if STAGE in (0, 2):
run("ElmerSolver wires2d_weight.sif")
run("%s dotproduct.py" % DOT_PY)
run("./paraview_dotproduct.py")
os.system("cp wires2d_weight.sif wires2d/sif_files/wires2d_weight_%02d_%1.4f.sif"%(count,val))
os.system("cp wires2d/elfield_weight_t0001.vtu wires2d/vtu_files/elfield_weight_%02d_%1.4f.vtu"%(count,val))
os.system("cp wires2d/dotproduct.vtu wires2d/vtu_files/dotproduct_%02d_%1.4f.vtu"%(count,val))
os.system("cp DotProduct_output.png png/DotProduct_output_z_%02d_%1.4f.png"%(count,val))
os.system("cp DotProduct_quarter_output.png png/DotProduct_output_z_%02d_%1.4f_quarter.png"%(count,val))
# os.system("python3 garfield_sim.py") # os.system("python3 garfield_sim.py")
os.system("cp wires2d.msh wires2d/mesh_files/wires2d%02d_%1.4f.msh"%(count,val))
os.system("cp wires2d.sif wires2d/sif_files/wires2d_%02d_%1.4f.sif"%(count,val)) # (a bare "break" here runs only the first z, as the original did)
os.system("cp wires2d/elfield_anasen_t0001.vtu wires2d/vtu_files/elfield_anasen_%02d_%1.4f.vtu"%(count,val))
os.system("cp contour_output.png png/Contour_output_z_%02d_%1.4f.png"%(count,val))
os.system("cp Field_output.png png/Field_ouput_z_%02d_%1.4f.png"%(count,val))
val=val+17.43
count = count + 1
break
# os.system("tar -cvzf wiress2d/mesh.tar.gz wires2d/mesh_files") # os.system("tar -cvzf wiress2d/mesh.tar.gz wires2d/mesh_files")
# os.system("rm -rf wires2d/mesh_files/*") # os.system("rm -rf wires2d/mesh_files/*")

4
anasen_fem/scalars.dat.names
View File

@ -1,8 +1,8 @@
Metadata for SaveScalars file: ./scalars.dat Metadata for SaveScalars file: ./scalars.dat
Elmer version: 26.2 Elmer version: 26.2
Elmer compilation date: 2026-05-14 Elmer compilation date: 2026-06-10
Solver input file: wires2d.sif Solver input file: wires2d.sif
File started at: 2026/05/19 22:16:54 File started at: 2026/09/22 20:59:24
Variables in columns of matrix: Variables in columns of matrix:
1: res: potential difference 1: res: potential difference

1
anasen_fem/scalars_weight.dat Normal file
View File

@ -0,0 +1 @@
1.000000000000E+000

8
anasen_fem/scalars_weight.dat.names Normal file
View File

@ -0,0 +1,8 @@
Metadata for SaveScalars file: ./scalars_weight.dat
Elmer version: 26.2
Elmer compilation date: 2026-06-10
Solver input file: wires2d_weight.sif
File started at: 2026/09/22 21:21:43
Variables in columns of matrix:
1: res: potential difference

7
anasen_fem/wires2d.sif
View File

@ -118,3 +118,10 @@ Boundary Condition 6
Target Boundaries = 30 Target Boundaries = 30
Potential = 0 Potential = 0
End End
! The one cathode singled out in wires_gmsh2d_bc.py. Grounded here, exactly
! like the other 23; it is only separate so wires2d_weight.sif can drive it.
Boundary Condition 7
Target Boundaries = 40
Potential = 0
End

132
anasen_fem/wires2d_weight.sif Normal file
View File

@ -0,0 +1,132 @@
Check Keywords Warn
! Ramo weighting potential for the single cathode at tag 40: that electrode
! sits at 1 V and every other one at 0 V. Identical to wires2d.sif apart from
! the potentials and the output names, so both solves discretise the same way
! and E . E_w is meaningful node by node.
Header
Mesh DB "." "wires2d"
End
Simulation
Coordinate System = Cartesian 2D
Simulation Type = Steady State
Steady State Max Iterations = 1
Output File = "elstatics_weight.result"
Post File = "elstatics_weight.ep"
Coordinate Scaling = 0.001 ! Converts mm from Gmsh to meters for Elmer
End
Constants
Permittivity Of Vacuum = 8.8542e-12
End
Body 1
Target Bodies(1) = 13
Equation = 1
Material = 1
End
Equation 1
Active Solvers(2) = 1 2
End
Material 1
Relative Permittivity = 1
End
Solver 1
Equation = Electrostatics
Procedure = "StatElecSolve" "StatElecSolver"
Variable = Potential
Variable DOFs = 1
Calculate Electric Field = True
Calculate Electric Flux = False
Linear System Solver = Iterative
Linear System Iterative Method = CG
Linear System Max Iterations = 5000
Linear System Convergence Tolerance = 1.0e-8
Linear System Preconditioning = ILU1
Calculate Vectors = Logical True
End
Solver 2
Equation = "Electric Field"
Procedure = "FluxSolver" "FluxSolver"
! Calculate from the potential
Target Variable = "Potential"
! Name of the output vector field in VTU
Flux Variable = String "Electric Field"
! Use 2D components (x, y)
Flux Coefficient = String "Permittivity"
Calculate Vectors = Logical True
End
Solver 3
Equation = Result Output
Procedure = "ResultOutputSolve" "ResultOutputSolver"
Output File Name = elfield_weight ! Sets prefix for output files
Output Format = Vtu
! Optional: Select specific variables to save
Scalar Field 1 = Potential
Vector Field 1 = Electric Field
End
Solver 4
Exec Solver = After All
Equation = SaveScalars
Procedure = "SaveData" "SaveScalars"
Filename = "scalars_weight.dat"
End
Boundary Condition 1
Target Boundaries = 1
Potential = 0
Calculate Electric Force = True
End
Boundary Condition 2
Target Boundaries = 2
Potential = 0
Calculate Electric Force = True
End
Boundary Condition 3
Target Boundaries = 3
Potential = 0
Calculate Electric Force = True
End
Boundary Condition 4
Target Boundaries = 10
Potential = 0
End
Boundary Condition 5
Target Boundaries = 20
Potential = 0
Calculate Electric Force = True
End
Boundary Condition 6
Target Boundaries = 30
Potential = 0
End
! The one cathode singled out in wires_gmsh2d_bc.py. Grounded here, exactly
! like the other 23; it is only separate so wires2d_weight.sif can drive it.
Boundary Condition 7
Target Boundaries = 40
Potential = 1
End

69
anasen_fem/wires_gmsh2d_bc.py
View File

@ -16,16 +16,24 @@ gmsh.option.setNumber("General.NumThreads", 10)
# gmsh.option.setNumber("Mesh.MeshSizeMax", 10.0) # gmsh.option.setNumber("Mesh.MeshSizeMax", 10.0)
gmsh.option.setNumber("Geometry.Tolerance", 4e-2) gmsh.option.setNumber("Geometry.Tolerance", 4e-2)
# gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) # gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
# Only needed alongside mesh.recombine(), which is off. See the note there.
# gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
# gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.option.setNumber("Mesh.SecondOrderIncomplete", 1)
lc = 0.04 lc = 0.04
# z_loc = -174.3 # z_loc = -174.3
if len(sys.argv) < 2: if len(sys.argv) < 2:
print("Usage: python3 wires_gmsh2d_bc.py <z_locus in mm>") print("Usage: python3 wires_gmsh2d_bc.py <z_locus in mm> [selected_cathode_index]")
quit() quit()
z_loc = float(sys.argv[1]) z_loc = float(sys.argv[1])
# Cathode wire (0..23) given its own physical group, tag 40, so that
# wires2d_weight.sif can hold it at 1 V while everything else stays grounded.
selected_cathode = int(sys.argv[2]) if len(sys.argv) > 2 else 1
wireShift = 4.0 wireShift = 4.0
k = 2 * np.pi / 24.0 k = 2 * np.pi / 24.0
kg = k/2.0 kg = k/2.0
@ -193,18 +201,44 @@ def get_surfs(disks):
return surfs return surfs
if not 0 <= selected_cathode < len(cathode_wires):
print("selected_cathode must be 0..%d" % (len(cathode_wires) - 1))
quit()
cathode_wires_other = [
disk for i, disk in enumerate(cathode_wires) if i != selected_cathode
]
# Cut the wire disks into the gas disk, replacing the old mesh.embed() call.
# embed only makes the mesher honour the curves; fragment splits the barrel
# into conforming surfaces, which is what lets one wire carry its own BC.
gmsh.model.occ.synchronize()
all_wire_disks = needle + guard_wires + cathode_wires + anode_wires
if include_ic_wires:
all_wire_disks += ic1_wires + ic2_wires
gmsh.model.occ.fragment([(2, anasen_barrel)], [(2, d) for d in all_wire_disks])
gmsh.model.occ.synchronize()
# fragment re-derives the bounding curves, so re-extract them.
needle_surfs = get_surfs(needle) if include_needle else [] needle_surfs = get_surfs(needle) if include_needle else []
gwire_surfs = get_surfs(guard_wires) gwire_surfs = get_surfs(guard_wires)
awire_surfs = get_surfs(anode_wires) awire_surfs = get_surfs(anode_wires)
cwire_surfs = get_surfs(cathode_wires) cwire_surfs = get_surfs(cathode_wires_other)
cwire_sel_surfs = get_surfs([cathode_wires[selected_cathode]])
i1wire_surfs = get_surfs(ic1_wires) if include_ic_wires else [] i1wire_surfs = get_surfs(ic1_wires) if include_ic_wires else []
i2wire_surfs = get_surfs(ic2_wires) if include_ic_wires else [] i2wire_surfs = get_surfs(ic2_wires) if include_ic_wires else []
all_active_wire_surfs = ( all_active_wire_surfs = (
needle_surfs + gwire_surfs + awire_surfs + cwire_surfs + i1wire_surfs + i2wire_surfs needle_surfs
+ gwire_surfs
+ awire_surfs
+ cwire_surfs
+ cwire_sel_surfs
+ i1wire_surfs
+ i2wire_surfs
) )
gmsh.model.mesh.embed(1, all_active_wire_surfs, 2, anasen_barrel)
f1 = gmsh.model.mesh.field.add("Distance") f1 = gmsh.model.mesh.field.add("Distance")
gmsh.model.mesh.field.setNumbers(f1, "CurvesList", all_active_wire_surfs) gmsh.model.mesh.field.setNumbers(f1, "CurvesList", all_active_wire_surfs)
@ -220,6 +254,10 @@ gmsh.model.mesh.field.setAsBackgroundMesh(f2)
# --- Physical Groups --- # --- Physical Groups ---
# These tags are the Target Bodies / Target Boundaries numbers in the .sif
# files. Do not renumber them without changing both sifs -- and note that
# ElmerGrid's -autoclean renumbers them for you, which silently breaks every
# boundary condition.
# Needle # Needle
if include_needle: if include_needle:
gmsh.model.addPhysicalGroup(1, needle_surfs, tag=1, name="hot_needle") gmsh.model.addPhysicalGroup(1, needle_surfs, tag=1, name="hot_needle")
@ -234,15 +272,24 @@ gmsh.model.addPhysicalGroup(1, gwire_surfs, tag=10, name="guard_wires")
gmsh.model.addPhysicalGroup(1, awire_surfs, tag=20, name="anode_wires") gmsh.model.addPhysicalGroup(1, awire_surfs, tag=20, name="anode_wires")
gmsh.model.addPhysicalGroup(1, cwire_surfs, tag=30, name="cathode_wires") gmsh.model.addPhysicalGroup(1, cwire_surfs, tag=30, name="cathode_wires")
# Gas Volume (2D) # The one cathode the weighting potential is solved for, kept out of tag 30 so
gmsh.model.addPhysicalGroup(2, [anasen_barrel], tag=13, name="gas") # wires2d_weight.sif can drive it independently.
gmsh.model.addPhysicalGroup(1, cwire_sel_surfs, tag=40, name="cathode_wire_selected")
gmsh.option.setNumber("Mesh.Algorithm", 6) # Gas Volume (2D). fragment split the barrel into many surfaces, so collect
# them all; naming the original disk alone would solve one sliver of the domain.
all_surfaces_2d = [s[1] for s in gmsh.model.getEntities(dim=2)]
gmsh.model.addPhysicalGroup(2, all_surfaces_2d, tag=13, name="gas")
gmsh.model.mesh.generate(dim=2) gmsh.option.setNumber("Mesh.Algorithm", 5)
# gmsh.model.mesh.refine()
# gmsh.model.mesh.refine() gmsh.model.mesh.generate(2)
gmsh.model.mesh.setOrder(2) gmsh.model.mesh.setOrder(2)
# recombine() ran for over half an hour on the barrel surface without
# finishing, and refine() took the mesh to ~16M nodes. Neither is needed: the
# Threshold field above already gives 0.05 mm elements on 0.254 mm wires.
# gmsh.model.mesh.recombine()
# gmsh.model.mesh.refine()
gmsh.write("wires2d.msh") gmsh.write("wires2d.msh")
# gmsh.fltk.run() # gmsh.fltk.run()
gmsh.finalize() gmsh.finalize()