.gitignore

@@ -32,5 +32,4 @@ Armory/CorrelateQQQ.h
 QQQStage2.C
 anasen_fem/scalars.dat.names
 myenv/
-eloss_calculations/
+eloss_calculations/
-anasen_fem/He96_CO2_4_260Torr.gas

anasen_fem/garfield_sim.py

@@ -4,7 +4,7 @@ import sys
 
 # 1. FIX: Manually load the Garfield library if it's not in the ROOT namespace
 # Update this path to your actual installation location
-garfield_lib_path = "/home/vs19g/garfieldpp/install/lib/libGarfield.so"
+garfield_lib_path = "/home/vsitaraman/garfieldpp/install/lib/libGarfield.so"
 
 if os.path.exists(garfield_lib_path):
     ROOT.gSystem.Load(garfield_lib_path)
@@ -44,13 +44,12 @@ else:
 # --- 3. FIELD MAP SETUP ---
 fm = ROOT.Garfield.ComponentElmer()
 
-
+# Update these filenames to match your Elmer SIF output exactly
-fm.Initialise("wires2d/mesh.header", 
+# Assuming ElmerGrid was run on 'wires2d' directory
+fm.Initialise("wires2d/mesh.nodes", 
               "wires2d/mesh.elements", 
-              "wires2d/mesh.nodes", 
+              "wires2d/mesh.boundary", 
-              "wires2d/dielectrics.dat", # Dielectrics (leave as empty string)
+              "wires2d/elfield_anasen.result", "mm")
-              "wires2d/elstatics.result", 
-              "mm")
 
 # Set the medium (Body 13 from your Gmsh script)
 fm.SetMedium(0, gas)
@@ -74,33 +73,8 @@ x0, y0, z0, t0 = 35.0, 0.0, 0.0, 0.0
 print(f"Simulating heavy ion drift from r={x0}...")
 drift.DriftIon(x0, y0, z0, t0)
 
-# Create a file to store the heavy ion track
-with open("heavy_ion_track.csv", "w") as f:
-    f.write("x,y,z,t\n")
-    
-    # After running drift.DriftIon(x0, y0, z0, t0):
-    n_points = drift.GetNumberOfDriftLinePoints()
-    for i in range(n_points):
-        xi, yi, zi, ti = ROOT.double(0), ROOT.double(0), ROOT.double(0), ROOT.double(0)
-        drift.GetDriftLinePoint(i, xi, yi, zi, ti)
-        f.write(f"{xi},{yi},{zi},{ti}\n")
-
 print(f"Simulating electron avalanche from r={x0}...")
 # AvalancheElectron(x, y, z, t, energy, dx, dy, dz)
 aval.AvalancheElectron(x0, y0, z0, t0, 0.1, 0.0, 0.0, 0.0)
 
-with open("avalanche_endpoints.csv", "w") as f:
-    f.write("x,y,z,t\n")
-    
-    # After aval.AvalancheElectron(...)
-    n_endpoints = aval.GetNumberOfEndpoints()
-    for i in range(n_endpoints):
-        # Get start and end points of each electron in the avalanche
-        x1, y1, z1, t1, e1 = ROOT.double(0), ROOT.double(0), ROOT.double(0), ROOT.double(0), ROOT.double(0)
-        x2, y2, z2, t2, e2, status = ROOT.double(0), ROOT.double(0), ROOT.double(0), ROOT.double(0), ROOT.double(0), ROOT.int(0)
-        
-        aval.GetEndpoint(i, x1, y1, z1, t1, e1, x2, y2, z2, t2, e2, status)
-        # We save the endpoint (x2, y2, z2) where the electron was collected or attached
-        f.write(f"{x2},{y2},{z2},{t2}\n")
-
 print("Simulation complete.")

anasen_fem/paraview_plotter.py

@@ -1,4 +1,4 @@
-#!/home/vs19g/ParaView-6.1.0-MPI-Linux-Python3.12-x86_64/bin/pvbatch
+#!/home/vsitaraman/ParaView-6.1.0-RC1-MPI-Linux-Python3.12-x86_64/bin/pvbatch
 import numpy as np
 import sys
 from paraview.simple import *
@@ -35,14 +35,15 @@ contour_display.SetScalarBarVisibility(renderView, True)
 view = GetActiveView()
 
 # 2. Define your desired coordinate ranges (x_min, x_max, y_min, y_max, z_min, z_max)
-x_min, x_max = -0.05, 0.05
+# Example: Look at a box from -10 to 10 in all dimensions
-y_min, y_max = -0.05, 0.05
+x_min, x_max = -50.0, 50.0
-z_min, z_max = -0.05, 0.05
+y_min, y_max = -50.0, 50.0
+z_min, z_max = -50.0, 50.0
 
 # 3. Calculate Center, Position, and Parallel Scale
 center = [(x_min + x_max) / 2.0, (y_min + y_max) / 2.0, (z_min + z_max) / 2.0]
 # Position the camera far away along Z to look at the center
-position = [center[0], center[1], 1.0] 
+position = [center[0], center[1], z_min - 30.0] 
 # Parallel scale defines how much of the scene is visible. 
 # It is usually half the height of the viewed area.
 view.CameraParallelScale = max((x_max - x_min), (y_max - y_min))/1.6
@@ -68,13 +69,15 @@ contour_display.RenderLinesAsTubes = 0    # Makes lines look smoother at high re
 # 1. Get the active view
 view = GetActiveView()
 
-# 1. Set the Focal Point to the middle of the quadrant in metres
+# 4. Apply settings
-zoom_center = [-0.025, 0.025, 0.0] 
+# 1. Set the Focal Point to the middle of the quadrant
+zoom_center = [-25, 25, 0.0] 
 
-# 2. Tighten the Parallel Scale 
+# 2. Tighten the Parallel Scale
-view.CameraParallelScale = 0.015 
+view.CameraParallelScale = 15 
 
-# 3. Position the Camera (0.5m away is fine)
+# 3. Position the Camera
+# Keep it 0.5m away looking "down" at the new center
 view.CameraPosition = [zoom_center[0], zoom_center[1], 0.5]
 view.CameraFocalPoint = zoom_center
 view.CameraViewUp = [0.0, 1.0, 0.0]
@@ -91,8 +94,9 @@ glyph = Glyph(Input=contour_filter, GlyphType='Arrow') #
 # Orientation Array: Use the 'electric field' vector from Elmer
 glyph.OrientationArray = ['POINTS', 'electric field']
 glyph.ScaleArray = ['POINTS', 'No scale array']
-glyph.ScaleFactor = 0.001  
+glyph.ScaleFactor = 1
 
+# Sampling: Every nth point (Stride 16)
 glyph.GlyphMode = 'Every Nth Point'
 glyph.Stride = 24
 

anasen_fem/run.py

@@ -7,10 +7,9 @@ count=11
 while val<178.3+0.1:
     print(val)
     os.system("python3 wires_gmsh2d_bc.py "+str(val))
-    os.system("ElmerGrid 14 2 wires2d.msh -2d")
+    os.system("ElmerGrid 14 2 wires2d.msh")
     os.system("ElmerSolver wires2d.sif")
     os.system("./paraview_plotter.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))
     os.system("cp wires2d/elfield_anasen_t0001.vtu wires2d/vtu_files/elfield_anasen_%02d_%1.4f.vtu"%(count,val))

anasen_fem/wires2d.sif

@@ -10,7 +10,6 @@ Simulation
   Steady State Max Iterations = 1
   Output File = "elstatics.result"
   Post File = "elstatics.ep"
-  Coordinate Scaling = 0.001  ! Converts mm from Gmsh to meters for Elmer
 End
 
 Constants
@@ -86,7 +85,7 @@ Solver 4
 End
 
 Boundary Condition 1
-  Target Boundaries = 1 
+  Target Boundaries = 1
   Potential = -200
   Calculate Electric Force = True
 End

anasen_fem/wires_gmsh2d_bc.py

@@ -188,17 +188,17 @@ 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
 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)
 
-f2 = gmsh.model.mesh.field.add("Threshold")
+# f2 = gmsh.model.mesh.field.add("Threshold")
-gmsh.model.mesh.field.setNumber(f2, "InField", f1)
+# gmsh.model.mesh.field.setNumber(f2, "InField", f1)
-gmsh.model.mesh.field.setNumber(f2, "SizeMin", 0.05)   # Fine mesh near wires
+# gmsh.model.mesh.field.setNumber(f2, "SizeMin", 0.1)   # Fine mesh near wires
-gmsh.model.mesh.field.setNumber(f2, "SizeMax", 5.0)  # Large mesh in empty space
+# gmsh.model.mesh.field.setNumber(f2, "SizeMax", 10.0)  # Large mesh in empty space
-gmsh.model.mesh.field.setNumber(f2, "DistMin", 0.5)   # Apply SizeMin within 1mm
+# gmsh.model.mesh.field.setNumber(f2, "DistMin", 1.0)   # Apply SizeMin within 1mm
-gmsh.model.mesh.field.setNumber(f2, "DistMax", 15.0)  # Transition to SizeMax by 20mm
+# gmsh.model.mesh.field.setNumber(f2, "DistMax", 20.0)  # Transition to SizeMax by 20mm
 
-gmsh.model.mesh.field.setAsBackgroundMesh(f2)
+# gmsh.model.mesh.field.setAsBackgroundMesh(f2)
 
 
 # --- Physical Groups ---
@@ -222,9 +222,8 @@ gmsh.model.addPhysicalGroup(2, [anasen_barrel], tag=13, name="gas")
 gmsh.option.setNumber("Mesh.Algorithm", 6)
 
 gmsh.model.mesh.generate(dim=2)
-# gmsh.model.mesh.refine()
+gmsh.model.mesh.refine()
-# gmsh.model.mesh.refine()
+gmsh.model.mesh.refine()
 gmsh.write("wires2d.msh")
-gmsh.model.mesh.setOrder(1)
 #gmsh.fltk.run()
 gmsh.finalize()