From 3edbfa5bd865c3b681e38d200f1b29c9b990c296 Mon Sep 17 00:00:00 2001 From: vsitaraman Date: Wed, 23 Sep 2026 12:14:05 -0400 Subject: [PATCH] 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 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 --- anasen_fem/README.md | 21 ++- anasen_fem/clean.sh | 4 +- anasen_fem/dotproduct.py | 148 +++++++++++++++ anasen_fem/paraview_dotproduct.py | 279 ++++++++++++++++++++++++++++ anasen_fem/paraview_plotter.py | 170 ++++++++++++++++- anasen_fem/run.py | 121 ++++++++++-- anasen_fem/scalars.dat.names | 4 +- anasen_fem/scalars_weight.dat | 1 + anasen_fem/scalars_weight.dat.names | 8 + anasen_fem/wires2d.sif | 7 + anasen_fem/wires2d_weight.sif | 132 +++++++++++++ anasen_fem/wires_gmsh2d_bc.py | 69 +++++-- 12 files changed, 920 insertions(+), 44 deletions(-) create mode 100644 anasen_fem/dotproduct.py create mode 100755 anasen_fem/paraview_dotproduct.py create mode 100644 anasen_fem/scalars_weight.dat create mode 100644 anasen_fem/scalars_weight.dat.names create mode 100644 anasen_fem/wires2d_weight.sif diff --git a/anasen_fem/README.md b/anasen_fem/README.md index e124e33..7c42725 100755 --- a/anasen_fem/README.md +++ b/anasen_fem/README.md @@ -23,5 +23,24 @@ - 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` - `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 (`_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 + - Sweep the weighting solve over all 24 cathodes rather than one at a time diff --git a/anasen_fem/clean.sh b/anasen_fem/clean.sh index f0c8051..a3feef1 100755 --- a/anasen_fem/clean.sh +++ b/anasen_fem/clean.sh @@ -1,2 +1,2 @@ -rm wires2d.msh -rm wires2d/* +rm -f wires2d.msh +rm -rf wires2d/* diff --git a/anasen_fem/dotproduct.py b/anasen_fem/dotproduct.py new file mode 100644 index 0000000..358d7d7 --- /dev/null +++ b/anasen_fem/dotproduct.py @@ -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() diff --git a/anasen_fem/paraview_dotproduct.py b/anasen_fem/paraview_dotproduct.py new file mode 100755 index 0000000..603a3c7 --- /dev/null +++ b/anasen_fem/paraview_dotproduct.py @@ -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) + diff --git a/anasen_fem/paraview_plotter.py b/anasen_fem/paraview_plotter.py index 43439f9..12866d0 100755 --- a/anasen_fem/paraview_plotter.py +++ b/anasen_fem/paraview_plotter.py @@ -4,10 +4,55 @@ import numpy as np import sys 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.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.ViewSize = [2000,2000] @@ -20,7 +65,7 @@ renderView.MultiSamples = 8 # 0 disables it, 4-8 is usually sufficient ResetCamera() 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 #colorbar contour_display_potentialLUT = GetColorTransferFunction('potential', contour_display, separate=True) @@ -62,23 +107,30 @@ Render() SaveScreenshot("contour_output.png") -#make glyps for field lines -contour_display.LineWidth = 1.0 # Increase this for thicker lines -contour_display.RenderLinesAsTubes = 0 # Makes lines look smoother at high res +# Everything from here on is a quadrant view: denser levels, heavier lines. +contour_filter.Isosurfaces = contour_levels(N_CONTOURS_QUARTER, CONTOUR_BIAS_QUARTER) +contour_display.LineWidth = LINE_WIDTH_QUARTER +contour_display.RenderLinesAsTubes = 1 # 1. Get the active view view = GetActiveView() # 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 -view.CameraParallelScale = 0.015 +view.CameraParallelScale = ZOOM_SCALE # 3. Position the Camera (0.5m away is fine) view.CameraPosition = [zoom_center[0], zoom_center[1], 0.5] view.CameraFocalPoint = zoom_center 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.Scalars = ['POINTS', 'potential'] # pot_threshold.ThresholdMethod = 'Above Upper Threshold' @@ -95,7 +147,7 @@ glyph.ScaleArray = ['POINTS', 'No scale array'] glyph.ScaleFactor = 0.001 glyph.GlyphMode = 'Every Nth Point' -glyph.Stride = 32 +glyph.Stride = GLYPH_STRIDE # --- 3. Display the Glyphs --- 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 ColorBy(glyph_display, ('POINTS', 'potential')) 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 # the existing 'potential' scalar bar. @@ -116,3 +168,101 @@ glyph_display.SetScalarBarVisibility(renderView, False) # --- 4. Final Render --- Render() 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") diff --git a/anasen_fem/run.py b/anasen_fem/run.py index 5a0e15f..bf20198 100755 --- a/anasen_fem/run.py +++ b/anasen_fem/run.py @@ -1,28 +1,113 @@ -import code import os +import sys -# val=-174.3 -val=0 -count=10 -while val<174.3+0.1: +# Cathode wire (0..23) that wires2d_weight.sif drives to 1 V. Passed to the +# mesher, which puts that one wire in its own physical group, tag 40. +SELECTED_CATHODE = 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) - os.system("python3 wires_gmsh2d_bc.py "+str(val)) - os.system("ElmerGrid 14 2 wires2d.msh -2d") - os.system("ElmerSolver wires2d.sif") - os.system("./paraview_plotter.py") + os.system("mkdir -p wires2d/mesh_files wires2d/sif_files wires2d/vtu_files png") + + # Each stage archives its own output before the next one starts, so a + # 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("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 Field_output.png png/Field_ouput_z_%02d_%1.4f.png"%(count,val)) - val=val+17.43 - count = count + 1 - break + + # (a bare "break" here runs only the first z, as the original did) # os.system("tar -cvzf wiress2d/mesh.tar.gz wires2d/mesh_files") # os.system("rm -rf wires2d/mesh_files/*") # os.system("tar -cvzf wires2d/sif.tar.gz wires2d/sif_files") # os.system("rm -rf wires2d/sif_files/*") # os.system("tar -cvzf wires2d/vtu.tar.gz wires2d/vtu_files") -# os.system("rm -rf wires2d/vtu_files/*") \ No newline at end of file +# os.system("rm -rf wires2d/vtu_files/*") diff --git a/anasen_fem/scalars.dat.names b/anasen_fem/scalars.dat.names index 0d58be4..73d4898 100755 --- a/anasen_fem/scalars.dat.names +++ b/anasen_fem/scalars.dat.names @@ -1,8 +1,8 @@ Metadata for SaveScalars file: ./scalars.dat Elmer version: 26.2 -Elmer compilation date: 2026-05-14 +Elmer compilation date: 2026-06-10 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: 1: res: potential difference diff --git a/anasen_fem/scalars_weight.dat b/anasen_fem/scalars_weight.dat new file mode 100644 index 0000000..9776b61 --- /dev/null +++ b/anasen_fem/scalars_weight.dat @@ -0,0 +1 @@ + 1.000000000000E+000 diff --git a/anasen_fem/scalars_weight.dat.names b/anasen_fem/scalars_weight.dat.names new file mode 100644 index 0000000..2540132 --- /dev/null +++ b/anasen_fem/scalars_weight.dat.names @@ -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 diff --git a/anasen_fem/wires2d.sif b/anasen_fem/wires2d.sif index 229d0b4..11a8483 100755 --- a/anasen_fem/wires2d.sif +++ b/anasen_fem/wires2d.sif @@ -118,3 +118,10 @@ 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 = 0 +End diff --git a/anasen_fem/wires2d_weight.sif b/anasen_fem/wires2d_weight.sif new file mode 100644 index 0000000..92afaf1 --- /dev/null +++ b/anasen_fem/wires2d_weight.sif @@ -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 diff --git a/anasen_fem/wires_gmsh2d_bc.py b/anasen_fem/wires_gmsh2d_bc.py index e3cd584..f8c3ccb 100755 --- a/anasen_fem/wires_gmsh2d_bc.py +++ b/anasen_fem/wires_gmsh2d_bc.py @@ -16,16 +16,24 @@ gmsh.option.setNumber("General.NumThreads", 10) # gmsh.option.setNumber("Mesh.MeshSizeMax", 10.0) gmsh.option.setNumber("Geometry.Tolerance", 4e-2) # 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 # z_loc = -174.3 if len(sys.argv) < 2: - print("Usage: python3 wires_gmsh2d_bc.py ") + print("Usage: python3 wires_gmsh2d_bc.py [selected_cathode_index]") quit() 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 k = 2 * np.pi / 24.0 kg = k/2.0 @@ -193,18 +201,44 @@ def get_surfs(disks): 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 [] gwire_surfs = get_surfs(guard_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 [] i2wire_surfs = get_surfs(ic2_wires) if include_ic_wires else [] - 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") gmsh.model.mesh.field.setNumbers(f1, "CurvesList", all_active_wire_surfs) @@ -220,6 +254,10 @@ gmsh.model.mesh.field.setAsBackgroundMesh(f2) # --- 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 if include_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, cwire_surfs, tag=30, name="cathode_wires") -# Gas Volume (2D) -gmsh.model.addPhysicalGroup(2, [anasen_barrel], tag=13, name="gas") +# The one cathode the weighting potential is solved for, kept out of tag 30 so +# 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.model.mesh.refine() -# gmsh.model.mesh.refine() +gmsh.option.setNumber("Mesh.Algorithm", 5) + +gmsh.model.mesh.generate(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.fltk.run() gmsh.finalize()