Compare commits

..

No commits in common. "1a01f332e622604d9725a91c5f7782a0d66f3cdd" and "b936da724a30a45b606dae56c955c26ec59ac13a" have entirely different histories.

6 changed files with 35 additions and 61 deletions

1
.gitignore vendored
View File

@ -33,4 +33,3 @@ QQQStage2.C
anasen_fem/scalars.dat.names anasen_fem/scalars.dat.names
myenv/ myenv/
eloss_calculations/ eloss_calculations/
anasen_fem/He96_CO2_4_260Torr.gas

View File

@ -4,7 +4,7 @@ import sys
# 1. FIX: Manually load the Garfield library if it's not in the ROOT namespace # 1. FIX: Manually load the Garfield library if it's not in the ROOT namespace
# Update this path to your actual installation location # 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): if os.path.exists(garfield_lib_path):
ROOT.gSystem.Load(garfield_lib_path) ROOT.gSystem.Load(garfield_lib_path)
@ -44,13 +44,12 @@ else:
# --- 3. FIELD MAP SETUP --- # --- 3. FIELD MAP SETUP ---
fm = ROOT.Garfield.ComponentElmer() 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.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) # Set the medium (Body 13 from your Gmsh script)
fm.SetMedium(0, gas) 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}...") print(f"Simulating heavy ion drift from r={x0}...")
drift.DriftIon(x0, y0, z0, t0) 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}...") print(f"Simulating electron avalanche from r={x0}...")
# AvalancheElectron(x, y, z, t, energy, dx, dy, dz) # AvalancheElectron(x, y, z, t, energy, dx, dy, dz)
aval.AvalancheElectron(x0, y0, z0, t0, 0.1, 0.0, 0.0, 0.0) 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.") print("Simulation complete.")

View File

@ -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 numpy as np
import sys import sys
from paraview.simple import * from paraview.simple import *
@ -35,14 +35,15 @@ contour_display.SetScalarBarVisibility(renderView, True)
view = GetActiveView() view = GetActiveView()
# 2. Define your desired coordinate ranges (x_min, x_max, y_min, y_max, z_min, z_max) # 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 # 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] 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 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. # Parallel scale defines how much of the scene is visible.
# It is usually half the height of the viewed area. # It is usually half the height of the viewed area.
view.CameraParallelScale = max((x_max - x_min), (y_max - y_min))/1.6 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 # 1. Get the active view
view = GetActiveView() 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.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]
@ -91,8 +94,9 @@ glyph = Glyph(Input=contour_filter, GlyphType='Arrow') #
# Orientation Array: Use the 'electric field' vector from Elmer # Orientation Array: Use the 'electric field' vector from Elmer
glyph.OrientationArray = ['POINTS', 'electric field'] glyph.OrientationArray = ['POINTS', 'electric field']
glyph.ScaleArray = ['POINTS', 'No scale array'] 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.GlyphMode = 'Every Nth Point'
glyph.Stride = 24 glyph.Stride = 24

View File

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

View File

@ -10,7 +10,6 @@ Simulation
Steady State Max Iterations = 1 Steady State Max Iterations = 1
Output File = "elstatics.result" Output File = "elstatics.result"
Post File = "elstatics.ep" Post File = "elstatics.ep"
Coordinate Scaling = 0.001 ! Converts mm from Gmsh to meters for Elmer
End End
Constants Constants

View File

@ -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 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) 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 --- # --- Physical Groups ---
@ -222,9 +222,8 @@ gmsh.model.addPhysicalGroup(2, [anasen_barrel], tag=13, name="gas")
gmsh.option.setNumber("Mesh.Algorithm", 6) gmsh.option.setNumber("Mesh.Algorithm", 6)
gmsh.model.mesh.generate(dim=2) 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.write("wires2d.msh")
gmsh.model.mesh.setOrder(1)
#gmsh.fltk.run() #gmsh.fltk.run()
gmsh.finalize() gmsh.finalize()