import os
import re
from xgc_analysis.mesh import Mesh
from xgc_analysis.plane import Plane
from xgc_analysis.magnetic_field import MagneticField
from xgc_analysis.species import Species
[docs]
class Simulation:
def __init__(self, directories=["./"], is_stellarator=False, sim_is_axisymmetric=False):
"""
Initializes a Simulation instance for the gyrokinetic code XGC.
Args:
directories (list of str): List of directories to search for simulation data.
Default is ['./'].
is_stellarator (bool): True if the simulation is for a stellarator (default: False, i.e. tokamak).
sim_is_axisymmetric (bool): True if the simulation is axisymmetric (default: False).
The constructor searches the provided directories for all three files (actually directories):
- xgc.mesh.bp
- xgc.equil.bp
- xgc.bfield.bp
The first directory containing all three files is used to set up the Mesh and MagneticField
instances. If none of the directories contains all required files, a RuntimeError is raised.
"""
required_files = ["xgc.mesh.bp", "xgc.equil.bp", "xgc.bfield.bp", "input"]
found_dir = None
for d in directories:
if all(os.path.exists(os.path.join(d, f)) for f in required_files):
found_dir = d
break
if found_dir is None:
raise RuntimeError("None of the provided directories contains all required files: " +
", ".join(required_files))
self.data_directory = found_dir
# Build full file paths.
mesh_file = os.path.join(found_dir, "xgc.mesh.bp")
equil_file = os.path.join(found_dir, "xgc.equil.bp")
bfield_file = os.path.join(found_dir, "xgc.bfield.bp")
input_file = os.path.join(found_dir, "input")
# Parse the input file to get the simulation parameters
self.input_params = self._parse_namelist_file(input_file)
# Set up species list
self.species = self._initialize_species()
# Create the Mesh and MagneticField instances.
# It is assumed that the Mesh class is defined in mesh.py and accepts a filename and a flag for axisymmetry.
self.mesh = Mesh(is_axisymmetric=(not is_stellarator),data_dir=found_dir)
# MagneticField is assumed to be defined in magnetic_field.py.
# Its constructor accepts a plane instance (e.g., the first plane from the mesh) and the equil and bfield file paths.
self.magnetic_field = MagneticField(plane_instance=self.mesh.get_plane(0),data_dir=found_dir)
self.is_stellarator = is_stellarator
self.sim_is_axisymmetric = sim_is_axisymmetric
# def _parse_namelist_file(self,filename):
# namelists = {}
# current_namelist = None
#
# with open(filename, 'r') as f:
# lines = f.readlines()
# for line in lines:
# line = line.split("!", 1)[0].strip() # Remove comment part and trim whitespace
# if not line:
# continue
# if line.startswith("&"):
# current_nml = line[1:].strip().lower()
# namelists[current_nml] = {}
# elif line.startswith("/"):
# current_nml = None
# elif current_nml:
# keyval = line.split('=')
# if len(keyval) == 2:
# key = keyval[0].strip().lower()
# val = keyval[1].strip()
# namelists[current_nml][key] = val
# return namelists
def _parse_namelist_file(self, filepath):
with open(filepath, 'r') as f:
lines = f.readlines()
namelists = {}
current_nml = None
for line in lines:
# Strip comments
line = line.split("!", 1)[0].strip()
if not line:
continue
# Start or end of namelist
if line.startswith("&"):
current_nml = line[1:].strip().lower()
namelists[current_nml] = {}
elif line.startswith("/"):
current_nml = None
elif current_nml:
if "=" in line:
key, val_str = line.split("=", 1)
key = key.strip().lower()
val_str = val_str.strip()
# Split by space (unless quoted string)
raw_values = re.findall(r"'[^']*'|[^ ]+", val_str)
# Convert values
values = [self.convert_fortran_value(v) for v in raw_values]
# Collapse to scalar if length 1
namelists[current_nml][key] = values if len(values) > 1 else values[0]
return namelists
[docs]
def convert_fortran_value(self,val):
val = val.strip()
# Boolean
if val.lower() in [".true.", "true", "t"]:
return True
if val.lower() in [".false.", "false", "f"]:
return False
# String
if val.startswith("'") and val.endswith("'"):
return val.strip("'")
# Fortran float with D exponent
try:
return float(val.replace('D', 'E'))
except ValueError:
pass
# Default fallback
return val
def _initialize_species(self):
"""
Initialize Species objects based on input parameters.
"""
ptl_param = self.input_params.get("ptl_param", {})
n_species = len(ptl_param.get("ptl_mass_au", []))
return [Species(self, i) for i in range(n_species)]
# ------------------------------------------------------------------------------
# Example usage:
# ------------------------------------------------------------------------------
if __name__ == "__main__":
# Create a Simulation instance searching in the current directory (or a list of directories)
sim = Simulation(directories=["./"], is_stellarator=False, sim_is_axisymmetric=True)
print("Simulation data loaded from directory:", sim.data_directory)
print("Mesh and MagneticField instances have been set up.")