Genex-Eirene coupling

.gitignore

@@ -1,7 +1,22 @@
 #files ignored in git repository
 *~
 new_*
-test_data/new_*
 __pycache__/
 .hypothesis
-*/.hypothesis
+*/.hypothesis
+*.egg-info/
+*local.sh
+job.out
+genex.pid
+genex_output.log
+*.pyc
+params_in.txt
+input_sources_??????.nc
+eirene_sources_??????/
+debug_log.txt
+params_*out.txt
+notes.txt
+tests/genex/integration/input_files/fort.31
+*.out
+*.log
+profile-results/

README.md

@@ -1,17 +1,42 @@
-# Overview of Repo
-This repo contains a python library that can be used to couple gyrokinetic codes to the monte-carlo neutral code EIRENE.
+# Neutral Coupling
 
-Modules used to extract/process eirene data are stored in the common folder
+This repository contains the `neutral_coupling` Python package for coupling
+gyrokinetic plasma codes to the Monte Carlo neutral code EIRENE.
 
-Tests for the modules in common along with the test data are in the tests folder 
+GENE-X users should first activate the environment created by Tor-X, then run:
 
-Modules used to couple eirene to gkeyll are in the gkeyllCoupling folder.
+```bash
+python -m pip install -e ".[genex]"
+```
 
-Documentation with helpful details is in the doc folder
+Tor-X is an external prerequisite and is not installed or modified by this
+package.
 
-# Breakdown of Gkeyll Coupling module
-Within the gkeyllCoupling folder are three subfolders:
+Gkeyll users can install independently in a Python 3.10 or newer environment (further details in README in `neutral_coupling\gkeyll_coupling`):
+
+```bash
+python -m pip install -e ".[gkeyll]"
+```
+
+These editable installs make `neutral_coupling` importable without modifying
+`PYTHONPATH`.
+
+Backend-independent EIRENE modules are in `neutral_coupling/common`.
+
+Tests and their fixtures are in `tests`.
+
+Backend implementations are in `neutral_coupling/genex_coupling` and
+`neutral_coupling/gkeyll_coupling`.
+
+Documentation relating filenames used by this repository for transferring Eirene sources to the source type is in `doc\fort_filenames.tex`
+
+# Gkeyll coupling
+
+The Gkeyll package contains `geometry`, `io`, and `examples` subpackages.
 
     1. gkeyllGeometry: Contains files used to generate mappings to interpolate data between Gkeyll and EIRENE
     2. gkeyllIO: Contains files used to pass plasma data from Gkeyll to EIRENE and neutral/plasma source data from EIRENE to Gkeyll
-    3. gkeyllExamples: Contains some example uses of the files in gkeyllGeometry and gkeyllIO 
+    3. gkeyllExamples: Contains some example uses of the files in gkeyllGeometry and gkeyllIO
+
+
+The LLMs GPT and related Codex were used to help in creating the code used for the GENE-X coupling, primarily in the creation of unit tests.

neutral_coupling/__init__.py

@@ -0,0 +1,3 @@
+"""Tools for coupling gyrokinetic plasma codes to EIRENE."""
+
+__version__ = "0.1.0"

neutral_coupling/common/__init__.py

@@ -0,0 +1 @@
+"""Backend-independent EIRENE and SOLPS utilities."""

common/B2IO.py → neutral_coupling/common/b2_io.py

@@ -1,6 +1,4 @@
-import numpy as np
-import warnings
-import scipy.constants as pyconst
+import numpy as np
 from pathlib import Path
 
 class B2:

common/EireneInputParser.py → neutral_coupling/common/eirene_input_parser.py

@@ -1,6 +1,8 @@
 import re
+import math
 from pathlib import Path
-from extra_fort_schema import PSEUDO_SPECIES
+from .extra_fort_schema import PSEUDO_SPECIES
+import scipy.constants as pyconst
 
 class EireneInputParser:
     """
@@ -32,6 +34,13 @@ def __init__(self, input_file: str | Path):
             "bulk_ions": [],
             "electrons": ["ELECTRONS"]
         }
+        self.masses = {                    # particle class -> [species]
+            "atoms": [],
+            "molecules": [],
+            "test_ions": [],
+            "bulk_ions": [],
+            "electrons": [pyconst.electron_mass]
+        }
 
         # unified mapping: atomic + bulk species -> stratum
         self.volume_recombination = {"particle":{}, "energy":{}} # {species_name: stratum}
@@ -98,6 +107,7 @@ def _parse_species_block(self, start_idx):
         n = int(self.lines[i].strip())
         i += 1
         species = []
+        mass = []
         for _ in range(n):
             # advance until we hit a real species-definition line
             while i < len(self.lines) and not self._is_species_line(self.lines[i]):
@@ -123,10 +133,16 @@ def _parse_species_block(self, start_idx):
 
             stride = ints[9]   # true 10th integer
             species.append(name)
+            # The mass number immediately follows the species name. Reading
+            # it by token position avoids treating digits in names such as
+            # D2, T2+, or D2T2(B) as numeric fields.
+            name_index = tokens.index(name)
+            mass_number = int(tokens[name_index + 1])
+            mass.append(mass_number * pyconst.proton_mass)
             # jump over irrelevant lines
             i += 1 + 2 * stride
 
-        return species
+        return species, mass
 
     def _is_species_line(self, line: str) -> bool:
         tokens = line.split()
@@ -140,11 +156,18 @@ def parse_species(self):
         for i, line in enumerate(self.lines):
             for cls, header in self.PARTICLE_HEADERS.items():
                 if line.strip().startswith(header):
-                    sp = self._parse_species_block(i)
+                    sp, mass = self._parse_species_block(i)
                     if cls == "bulk_ions":
                         # remove (B) species
-                        sp = [s for s in sp if "(B)" not in s]
+                        physical = [
+                            (species, species_mass)
+                            for species, species_mass in zip(sp, mass)
+                            if "(B)" not in species
+                        ]
+                        sp = [species for species, _ in physical]
+                        mass = [species_mass for _, species_mass in physical]
                     self.species[cls].extend(sp)
+                    self.masses[cls].extend(mass)
 
     # ---------------------------
     # VOLUME RECOMBINATION
@@ -248,6 +271,40 @@ def validate_schema(self):
             if not sp_list:
                 self.report.warn(f"No species found for particle class '{cls}'")
 
+        # Every physical species must have exactly one finite, positive mass.
+        pseudo_species = set(PSEUDO_SPECIES.values())
+        for cls, sp_list in self.species.items():
+            physical_species = [sp for sp in sp_list if sp not in pseudo_species]
+            mass_list = self.masses.get(cls)
+
+            if mass_list is None:
+                self.report.error(
+                    f"No masses found for particle class '{cls}'"
+                )
+                continue
+
+            if len(mass_list) != len(physical_species):
+                self.report.error(
+                    f"Mass count mismatch for particle class '{cls}': "
+                    f"{len(physical_species)} physical species, "
+                    f"{len(mass_list)} masses"
+                )
+
+            for index, mass in enumerate(mass_list):
+                species = (
+                    physical_species[index]
+                    if index < len(physical_species)
+                    else f"mass index {index}"
+                )
+                if (
+                    not isinstance(mass, (int, float))
+                    or not math.isfinite(mass)
+                    or mass <= 0
+                ):
+                    self.report.error(
+                        f"Invalid mass for '{species}' in particle class "
+                        f"'{cls}': {mass!r}"
+                    )
 
         # volume recombination consistency
         if "SUM" in self.requested_strata:
@@ -282,6 +339,10 @@ def _print_debug_summary(self):
         for cls, sp in self.species.items():
             print(f"  {cls}: {sp}")
 
+        print("\nMasses by particle class:")
+        for cls, masses in self.masses.items():
+            print(f"  {cls}: {masses}")
+
         print("\nVolume recombination mapping:")
         if self.volume_recombination:
             for sp, st in self.volume_recombination.items():
@@ -319,4 +380,4 @@ def print(self):
         if self.warnings:
             print("\n[EIRENE INPUT VALIDATION WARNINGS]")
             for w in self.warnings:
-                print(" -", w)
+                print(" -", w)

common/eireneIO.py → neutral_coupling/common/eirene_io.py

@@ -1,14 +1,13 @@
-import StrataAssigner as SA
-import extra_fort_schema
-import triangle_mesh as triangles
+from . import strata_assigner as SA
+from . import extra_fort_schema
+from . import triangle_mesh as triangles
 
 import numpy as np
-import glob
 import scipy.constants as pyconst
-import platform
+import sys
 from pathlib import Path
+from copy import deepcopy
 # for debugging
-import pdb
 
 default_coll_to_adjust = {}
 default_coll_to_adjust["testion-plasma"] = True
@@ -22,11 +21,16 @@ def safe_float(s):
 
 class eirene:
     def __init__(self, filepath=None):
-        version = platform.python_version()
-        if not(int(version[0]) >= 3 and int(version[2:4])>=10):
-            print("Python version 3.10 or greater needed for eirene.load_extra_forts()")
+        if sys.version_info < (3, 10):
+            print("Python version 3.10 or greater needed for match/case in load_extra_forts")
             return
-        self.fort44 = {"meta":{}, "neut":{}, "wld":{}, "res":{}}
+        self.fort44 = {
+            "meta": {},
+            "neut": {},
+            "neut_int": {},
+            "wld": {},
+            "res": {},
+        }
         self.fort44_expanded = {}
         self.fort46 = {}
         self.triangle_mesh = triangles.triangle_mesh(filepath)
@@ -388,7 +392,7 @@ def load_extra_forts(self, eirene_path, extension="???",
                          requested_strata=None, expected_species=None,
                          vol_rec_mapping=None, debug=False,
                          coll_to_adjust=default_coll_to_adjust, print_info=False,
-                         species_patt_to_ignore=">"):
+                         species_patt_to_ignore=">", convert_units=False):
         if isinstance(eirene_path, str):
             eirene_path = Path(eirene_path)
         filelist = eirene_path.glob("fort."+extension)
@@ -415,7 +419,7 @@ def load_extra_forts(self, eirene_path, extension="???",
             header_lines = 8
             Ncells = int(lines_list[header_lines-1].split()[0])-1
             if (Ncells != ntria):
-                raise ValueError(f"Number of triangles from mesh unequal to number from {current_file}")
+                raise ValueError(f"Number of triangles from mesh ({ntria}) unequal to number from {current_file} ({Ncells})")
             add_cells = int(lines_list[header_lines-1].split()[4])-1 - Ncells
             start_line = header_lines
             while x*Ncells < len(lines_list):
@@ -442,6 +446,11 @@ def load_extra_forts(self, eirene_path, extension="???",
         self.loaded_sources = loader
         self.sources = loader.sum_over_collisions()
 
+        if convert_units:
+            self.convert_source_dict_to_SI(self.sources)
+            self.full_source_in_SI = deepcopy(self.loaded_sources.sources)
+            self.convert_source_dict_to_SI(self.full_source_in_SI)
+
     def read_ft30(self, filename):
         self.plasma_gmtry = {}
         with open(filename, "r") as f:
@@ -1073,3 +1082,26 @@ def __write_ft31_field(self, fid, fieldname):
                             formatted_number = " "+mantissa+exponent
                         fid.write(formatted_number)
                     fid.write("\n")
+
+    def convert_source_dict_to_SI(self, d):
+        for root_key, subtree in d.items():
+            self._apply(subtree, root_key)
+
+    def _apply(self, d, root_key):
+        for k, v in d.items():
+            if isinstance(v, dict):
+                self._apply(v, root_key)
+            elif isinstance(v, np.ndarray):
+                d[k] = self._convert_source_to_SI(v, root_key)
+
+    def _convert_source_to_SI(self, arr, moment):
+        match moment:
+            case(0 | "0" | "n" | "particle"):
+                return arr*1e6/pyconst.elementary_charge
+            case(1 | "1" | "momentum"):
+                return arr*10/pyconst.elementary_charge
+            case(2 | "2" | "energy"):
+                return arr*1e6
+            case _:
+                raise ValueError(f"Unknown moment '{moment}' in "
+                                 f"_convert_source_to_SI")

common/extra_fort_schema.py → neutral_coupling/common/extra_fort_schema.py

@@ -37,4 +37,4 @@
 PSEUDO_SPECIES = MappingProxyType(PSEUDO_SPECIES_WRITABLE)
 PARTICLE_CLASS_MAP = MappingProxyType(PARTICLE_CLASS_MAP_WRITABLE)
 COLLISION_MAP = MappingProxyType(COLLISION_MAP_WRITABLE)
-MOMENT_MAP = MappingProxyType(MOMENT_MAP_WRITABLE)
+MOMENT_MAP = MappingProxyType(MOMENT_MAP_WRITABLE)