diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json new file mode 100644 index 0000000..6b7a787 --- /dev/null +++ b/.devcontainer/devcontainer.json @@ -0,0 +1,25 @@ +// For format details, see https://aka.ms/devcontainer.json. For config options, see the +// README at: https://github.com/devcontainers/templates/tree/main/src/python +{ + "name": "Python 3", + // Or use a Dockerfile or Docker Compose file. More info: https://containers.dev/guide/dockerfile + "image": "mcr.microsoft.com/devcontainers/python:2-3.9-bullseye", + "features": { + "ghcr.io/devcontainers/features/python:1": { + "installTools": true, + "version": "3.9" + } + }, + + // Use 'forwardPorts' to make a list of ports inside the container available locally. + // "forwardPorts": [], + + // Use 'postCreateCommand' to run commands after the container is created. + "postCreateCommand": "pip3 install --user -r requirements.txt" + + // Configure tool-specific properties. + // "customizations": {}, + + // Uncomment to connect as root instead. More info: https://aka.ms/dev-containers-non-root. + // "remoteUser": "root" +} diff --git a/.github/workflows/renate-od.yml b/.github/workflows/renate-od.yml index c2c8ffe..87e9c1d 100644 --- a/.github/workflows/renate-od.yml +++ b/.github/workflows/renate-od.yml @@ -5,15 +5,13 @@ on: jobs: build: - - runs-on: ubuntu-20.04 - + runs-on: ubuntu-24.04 steps: - uses: actions/checkout@v4 - - name: Set up Python 3.6.4 + - name: Set up Python 3.9.12 uses: actions/setup-python@v5 with: - python-version: '3.6.4' + python-version: '3.9.12' - name: Install dependencies run: | python -m pip install --upgrade pip @@ -23,11 +21,11 @@ jobs: python -m unittest -v crm_solver.odetest.OdeTest python -m unittest -v crm_solver.crmsystemtest.CrmRegressionTest python -m unittest -v crm_solver.crmsystemtest.CrmAcceptanceTest - python -m unittest -v crm_solver.atomic_dbtest.AtomicDBTest python -m unittest -v crm_solver.beamlettest.BeamletTest - python -m unittest -v crm_solver.atomic_dbtest.RenateDBTest - python -m unittest -v crm_solver.neutral_dbtest.NeutralDBTest python -m unittest -v crm_solver.coefficientmatrixtest.CoefficientMatrixTest + python -m unittest -v atomic.atomic_dbtest.AtomicDBTest + python -m unittest -v atomic.atomic_dbtest.RenateDBTest + python -m unittest -v atomic.neutral_dbtest.NeutralDBTest python -m unittest -v utility.accessdatatest.AccessDataTest python -m unittest -v utility.getdatatest.GetDataTest python -m unittest -v utility.putdatatest.PutDataTest diff --git a/atomic/aladdin.py b/atomic/aladdin.py index 76f32dd..cdbd3e2 100644 --- a/atomic/aladdin.py +++ b/atomic/aladdin.py @@ -1,12 +1,95 @@ -# import requests +import pycollisiondb +from pycollisiondb.pycollisiondb import PyCollision +from scipy.interpolate import interp1d +class FitFunction(object): + def __init__(self, data): + self.coefficients = self.assign_from_fit(data, 'coeffs') + self.elo = self.assign_from_fit(data, 'Elo') + self.ehi = self.assign_from_fit(data, 'Ehi') + self.name = self.assign_from_fit(data, 'func') + self.precision = self.assign_from_fit(data, 'fit_unc_perc') + + @staticmethod + def assign_from_fit(data, tag): + try: + return data[tag] + except KeyError: + return None + +class AladdinCrossSection(object): + def __init__(self, data): + if not isinstance(data, pycollisiondb.pycoll_ds.PyCollDataSet): + raise TypeError('The expected type from aladdin API is PyCollision') + self.energy = data.x + self.cross_section = data.y + self.len = data.n + self.reaction = self.assign_data(data, tag='reaction') + self.type = self.assign_data(data, tag='process_types') + self.uncertainty = float(self.assign_from_json(data, tag='unc_perc')) + #self.threshold = float(self.assign_from_json(data, tag='threshold')) + self.function = self.assign_from_json(data, tag='data_from_fit') + if self.function: + self.fit = FitFunction(data.metadata['json_data']['fit']) + + @staticmethod + def assign_from_json(data, tag): + try: + return data.metadata['json_data'][tag] + except KeyError: + return None + + @staticmethod + def assign_data(data, tag): + try: + return data.metadata[tag] + except KeyError: + return None class AladdinData(object): - def __init__(self): - pass + def __init__(self, url_source = 'aladdin'): + if isinstance(url_source, str): + if url_source == 'aladdin': + self.api_url = 'https://db-amdis.org/aladdin2/api/' + elif url_source == 'collisiondb': + self.api_url = 'https://db-amdis.org/collisiondb/api/' + else: + raise ValueError('The requested data source is not supported: ' + url_source) + else: + raise TypeError('The required data type for url specification is .') + + @staticmethod + def _set_sign(particle): + if particle.charge >= 1: + return '+' + else: + return '-' + + def rod_to_pycoll(self, transition): + reactants = str(transition.projectile) + ' ' + transition.from_level + ' + ' + \ + str(transition.target) + self._set_sign(transition.target) + + if transition.name == 'ex' or transition.name == 'de-ex': + products = str(transition.projectile) + ' ' + transition.to_level + ' + ' + \ + str(transition.target) + self._set_sign(transition.target) + elif transition.name == 'ion': + products = str(transition.projectile) + '+ + ' + \ + str(transition.target) + self._set_sign(transition.target) + ' + ' + 'e-' + elif transition.name == 'cx': + products = str(transition.projectile) + '+ + ' + \ + str(transition.target) + else: + products = str(transition.projectile) + '+ + ' + return reactants + ' -> ' + products - def fetch_data(self, transition): - pass + def get_data_from_iaea(self, pycoll_reaction): + pycoll_raw = PyCollision.get_datasets(query = {'reaction_text':pycoll_reaction, 'data_type':'cross section'}, + API_URL=self.api_url) + return AladdinCrossSection(pycoll_raw.datasets[list(pycoll_raw.datasets.keys())[0]]) - def check_data(self, transition): - pass + def get_cross_section(self, transition, energy_grid): + reaction = self.rod_to_pycoll(transition) + aladdin_data = self.get_data_from_iaea(pycoll_reaction=reaction) + cross_function = interp1d(aladdin_data.energy, aladdin_data.cross_section, + kind='cubic', fill_value='extrapolate') + return cross_function(energy_grid) diff --git a/atomic/cross_section.py b/atomic/cross_section.py index 8943815..ee79058 100644 --- a/atomic/cross_section.py +++ b/atomic/cross_section.py @@ -6,13 +6,13 @@ class CrossSection: - CROSSEC_SOURCES = {'aladdin': [], 'internal': ['hydrogenic']} + CROSSEC_SOURCES = {'aladdin': [], 'collisiondb':[], 'internal': ['hydrogenic']} def __new__(cls, source, projectile=None): if not isinstance(source, str): raise TypeError('A string is expected to decide which data source to pursue.\n'+CrossSection.GetAvailableData()) - if source == 'aladdin': - return AladdinData() + if source == 'aladdin' or source == 'collisiondb': + return AladdinData(url_source=source) elif source == 'internal': if not isinstance(projectile, str): raise TypeError('The data source requires the projectile parameter to be a valid string.\n'+CrossSection.GetAvailableData()) diff --git a/atomic/internal_db.py b/atomic/internal_db.py index e9f6f15..94d9f6c 100644 --- a/atomic/internal_db.py +++ b/atomic/internal_db.py @@ -33,8 +33,20 @@ def __projectile_parameters(self): self.projectile_particle = Particle(label='T', mass_number=3, atomic_number=1) self.atomic_dict = {'1': 0, '2': 1, '3': 2, '4': 3, '5': 4, '6': 5, None: None} self.atomic_levels = 6 - self.inv_atomic_dict = {index: name for name, index in self.atomic_dict.items()} - self.mass = self.projectile_particle.mass + elif self.projectile == 'Li': + self.projectile_type = 'alkali' + self.projectile_particle = Particle(label='Li', charge=0, mass_number=7, atomic_number=3) + self.atomic_dict = {'2s': 0, '2p': 1, '3s': 2, '3p': 3, '3d': 4, '4s': 5, '4p': 6, '4d': 7, '4f': 8} + self.atomic_levels = 9 + elif self.projectile == 'Na': + self.projectile_type = 'alkali' + self.projectile_particle = Particle(label='Na', charge=0, mass_number=23, atomic_number=11) + self.atomic_dict = {'3s': 0, '3p': 1, '3d': 2, '4s': 3, '4p': 4, '4d': 5, '4f': 6, '5s': 7} + self.atomic_levels = 8 + else: + raise ValueError('The requested projectile is not yet supported.') + self.inv_atomic_dict = {index: name for name, index in self.atomic_dict.items()} + self.mass = self.projectile_particle.mass def _get_projectile_velocity(self): self.velocity = uc.calculate_velocity_from_energy(uc.convert_keV_to_eV(float(self.energy)), self.mass) @@ -43,6 +55,10 @@ def _get_projectile_velocity(self): def set_default_atomic_levels(self): if self.projectile in ['H', 'D', 'T']: return '3', '2', '1', '3n-->2n' + elif self.projectile == 'Li': + return '2p', '2s', '2s', '2p-->2s' + elif self.projectile == 'Na': + return '3p', '3s', '3s', '3p-->3s' def get_rate_interpolator(self, reaction_type, target, from_level, to_level=None): #print(reaction_type, target, from_level, to_level) diff --git a/crm_solver/beamlettest.py b/crm_solver/beamlettest.py index d9b4212..6ec9e56 100644 --- a/crm_solver/beamlettest.py +++ b/crm_solver/beamlettest.py @@ -8,7 +8,7 @@ class BeamletTest(unittest.TestCase): EXPECTED_ATTR = ['param', 'components', 'profiles', 'coefficient_matrix', 'atomic_db', 'initial_condition'] - EXPECTED_INITIAL_CONDITION = [4832583106.4753895, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + EXPECTED_INITIAL_CONDITION = [4832583046.75342, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] EXPECTED_PARAM_ATTR = ['beamlet_source', 'beamlet_energy', 'beamlet_species', 'beamlet_current'] EXPECTED_COMPONENTS_KEYS = ['q', 'Z', 'A'] EXPECTED_COMPONENTS_SPECIES = ['electron', 'ion1', 'ion2'] @@ -42,8 +42,8 @@ def test_initial_conditions(self): for element in range(len(self.beamlet.initial_condition)): self.assertIsInstance(self.beamlet.initial_condition[element], float, msg='Expected type for initial' ' conditions is float.') - self.assertEqual(self.beamlet.initial_condition[element], self.EXPECTED_INITIAL_CONDITION[element], - msg='Computed Init conditions do not match expected init conditions.') + self.assertAlmostEqual(self.beamlet.initial_condition[element], self.EXPECTED_INITIAL_CONDITION[element], + delta=1e-1, msg='Computed Init conditions do not match expected init conditions.') def test_param_xml(self): self.assertIsInstance(self.beamlet.param, etree._ElementTree, diff --git a/requirements.txt b/requirements.txt index d06f94c..11d0e26 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,9 +1,10 @@ -scipy==1.0.0 -numpy==1.14.0 -pandas==0.22.0 -h5py==2.9.0 -matplotlib==2.1.2 -lxml==4.6.3 -tables==3.4.2 -paramiko==2.4.2 -scp==0.13.2 +scipy==1.8.0 +numpy==1.23.1 +pandas==1.4.2 +h5py==3.6.0 +matplotlib==3.6.0 +lxml==4.8.0 +tables==3.6.1 +paramiko==2.8.1 +scp==0.14.4 +pycollisiondb==0.1.3 diff --git a/utility/getdata.py b/utility/getdata.py index 45249bc..bb24106 100644 --- a/utility/getdata.py +++ b/utility/getdata.py @@ -87,12 +87,8 @@ def read_h5_to_array(self): print('Data could NOT be read to array from HD5 file: ' + self.access_path + '. Key is missing!') raise ValueError try: - with h5py.File(self.access_path, 'r') as hdf5_id: - hdf5_group = hdf5_id - for key in self.data_key: - hdf5_group = hdf5_group[key] - self.data = hdf5_group.value - hdf5_id.close() + with h5py.File(self.access_path, 'r') as file: + self.data = numpy.array(file[self.data_key[0]]) print("Data read to array from HD5 file: " + self.access_path + " with key: " + str(self.data_key)) except ValueError: print("Data could NOT be read to array from HD5 file: " + self.access_path +