diff --git a/cherab/openadas/install.py b/cherab/openadas/install.py index 8ed02f6..12ed67c 100644 --- a/cherab/openadas/install.py +++ b/cherab/openadas/install.py @@ -28,6 +28,12 @@ def install_files(configuration, download=False, repository_path=None, adas_path=None): for adf in configuration: + if adf.lower() == 'adf11scd': + for args in configuration[adf]: + install_adf11scd(*args, download=download, repository_path=repository_path, adas_path=adas_path) + if adf.lower() == 'adf11acd': + for args in configuration[adf]: + install_adf11acd(*args, download=download, repository_path=repository_path, adas_path=adas_path) if adf.lower() == 'adf12': for args in configuration[adf]: install_adf12(*args, download=download, repository_path=repository_path, adas_path=adas_path) @@ -45,6 +51,48 @@ def install_files(configuration, download=False, repository_path=None, adas_path install_adf22bme(*args, download=download, repository_path=repository_path, adas_path=adas_path) +def install_adf11scd(element, file_path, download=False, repository_path=None, adas_path=None): + """ + Adds the ionisation rate defined in an ADF11 file to the repository. + + :param element: The element described by the rate file. + :param file_path: Path relative to ADAS root. + :param download: Attempt to download file if not present (Default=True). + :param repository_path: Path to the repository in which to install the rates (optional). + :param adas_path: Path to ADAS files repository (optional). + """ + + print('Installing {}...'.format(file_path)) + path = _locate_adas_file(file_path, download, adas_path) + if not path: + raise ValueError('Could not locate the specified ADAS file.') + + # decode file and write out rates + rate = parse_adf11(element, path) + repository.update_ionisation_rates(rate, repository_path) + + +def install_adf11acd(element, file_path, download=False, repository_path=None, adas_path=None): + """ + Adds the recombination rate defined in an ADF11 file to the repository. + + :param element: The element described by the rate file. + :param file_path: Path relative to ADAS root. + :param download: Attempt to download file if not present (Default=True). + :param repository_path: Path to the repository in which to install the rates (optional). + :param adas_path: Path to ADAS files repository (optional). + """ + + print('Installing {}...'.format(file_path)) + path = _locate_adas_file(file_path, download, adas_path) + if not path: + raise ValueError('Could not locate the specified ADAS file.') + + # decode file and write out rates + rate = parse_adf11(element, path) + repository.update_recombination_rates(rate, repository_path) + + # todo: move print calls to logging def install_adf12(donor_ion, donor_metastable, receiver_ion, receiver_ionisation, file_path, download=False, repository_path=None, adas_path=None): """ diff --git a/cherab/openadas/openadas.py b/cherab/openadas/openadas.py index 0b880b4..ad7c73b 100644 --- a/cherab/openadas/openadas.py +++ b/cherab/openadas/openadas.py @@ -56,7 +56,37 @@ def wavelength(self, ion, ionisation, transition): ion = ion.element return repository.get_wavelength(ion, ionisation, transition) - def beam_cx_rate(self, donor_ion, receiver_ion, receiver_ionisation, transition): + def ionisation_rate(self, ion, ionisation): + + if isinstance(ion, Isotope): + ion = ion.element + + try: + data = repository.get_ionisation_rate(ion, ionisation) + + except RuntimeError: + if self._missing_rates_return_null: + return NullIonisationRate() + raise + + return IonisationRate(data, extrapolate=self._permit_extrapolation) + + def recombination_rate(self, ion, ionisation): + + if isinstance(ion, Isotope): + ion = ion.element + + try: + data = repository.get_recombination_rate(ion, ionisation) + + except RuntimeError: + if self._missing_rates_return_null: + return NullRecombinationRate() + raise + + return RecombinationRate(data, extrapolate=self._permit_extrapolation) + + def beam_cx_pec(self, donor_ion, receiver_ion, receiver_ionisation, transition): """ :param donor_ion: @@ -80,13 +110,13 @@ def beam_cx_rate(self, donor_ion, receiver_ion, receiver_ionisation, transition) except RuntimeError: if self._missing_rates_return_null: - return [NullBeamCXRate()] + return [NullBeamCXPEC()] raise # load and interpolate the relevant transition data from each file rates = [] for donor_metastable, rate_data in data: - rates.append(BeamCXRate(donor_metastable, wavelength, rate_data, extrapolate=self._permit_extrapolation)) + rates.append(BeamCXPEC(donor_metastable, wavelength, rate_data, extrapolate=self._permit_extrapolation)) return rates def beam_stopping_rate(self, beam_ion, plasma_ion, ionisation): @@ -146,7 +176,7 @@ def beam_population_rate(self, beam_ion, metastable, plasma_ion, ionisation): # load and interpolate data return BeamPopulationRate(data, extrapolate=self._permit_extrapolation) - def beam_emission_rate(self, beam_ion, plasma_ion, ionisation, transition): + def beam_emission_pec(self, beam_ion, plasma_ion, ionisation, transition): """ :param beam_ion: @@ -170,13 +200,13 @@ def beam_emission_rate(self, beam_ion, plasma_ion, ionisation, transition): except RuntimeError: if self._missing_rates_return_null: - return NullBeamEmissionRate() + return NullBeamEmissionPEC() raise # load and interpolate data - return BeamEmissionRate(data, wavelength, extrapolate=self._permit_extrapolation) + return BeamEmissionPEC(data, wavelength, extrapolate=self._permit_extrapolation) - def impact_excitation_rate(self, ion, ionisation, transition): + def impact_excitation_pec(self, ion, ionisation, transition): """ :param ion: @@ -194,12 +224,12 @@ def impact_excitation_rate(self, ion, ionisation, transition): except RuntimeError: if self._missing_rates_return_null: - return NullImpactExcitationRate() + return NullImpactExcitationPEC() raise - return ImpactExcitationRate(wavelength, data, extrapolate=self._permit_extrapolation) + return ImpactExcitationPEC(wavelength, data, extrapolate=self._permit_extrapolation) - def recombination_rate(self, ion, ionisation, transition): + def recombination_pec(self, ion, ionisation, transition): """ :param ion: @@ -217,10 +247,10 @@ def recombination_rate(self, ion, ionisation, transition): except (FileNotFoundError, KeyError): if self._missing_rates_return_null: - return NullRecombinationRate() + return NullRecombinationPEC() raise - return RecombinationRate(wavelength, data, extrapolate=self._permit_extrapolation) + return RecombinationPEC(wavelength, data, extrapolate=self._permit_extrapolation) # def stage_resolved_line_ radiation_rate(self, ion, ionisation): # diff --git a/cherab/openadas/parse/__init__.py b/cherab/openadas/parse/__init__.py index 002aebf..c7b3c3f 100644 --- a/cherab/openadas/parse/__init__.py +++ b/cherab/openadas/parse/__init__.py @@ -1,3 +1,5 @@ + +from .adf11 import parse_adf11 from .adf12 import parse_adf12 from .adf15 import parse_adf15 from .adf21 import parse_adf21 diff --git a/cherab/openadas/parse/adf11.py b/cherab/openadas/parse/adf11.py index 1457d42..76f4265 100644 --- a/cherab/openadas/parse/adf11.py +++ b/cherab/openadas/parse/adf11.py @@ -16,10 +16,11 @@ # See the Licence for the specific language governing permissions and limitations # under the Licence. -import numpy as np import re +import numpy as np -# todo: to be implemented in a future release +from cherab.core.atomic import Element +from cherab.core.utility import RecursiveDict, Cm3ToM3, PerCm3ToPerM3 # def parse_adf11acd(ion, ionisation, adf_file_path): @@ -58,6 +59,110 @@ # pass +def parse_adf11(element, adf_file_path): + """ + Reads contents of open adas adf11 files + + :param element: Element described by ADF file. + :param adf_file_path: Path to ADF11 file from ADAS root. + :return: temperature, density, rates as numpy array + """ + + if not isinstance(element, Element): + raise TypeError('The element must be an Element object.') + + with open(adf_file_path, "r") as source_file: + + lines = source_file.readlines() # read file contents by lines + tmp = re.split("\s{2,}", lines[0].strip()) # split into relevant variables + # exctract variables + z_nuclear = int(tmp[0]) + n_densities = int(tmp[1]) + n_temperatures = int(tmp[2]) + z_min = int(tmp[3]) + z_max = int(tmp[4]) + element_name = tmp[5].strip('/').lower() + projectname = tmp[6] + + if element.atomic_number != z_nuclear or element.name != element_name: + raise ValueError("The requested element '{}' does not match the element description in the" + "specified ADF11 file, '{}'.".format(element.name, element_name)) + + # check if it is a resolved file + if re.match("\s*[0-9]+", lines[3]): # is it unresolved? + startsearch = 2 + else: + startsearch = 4 # skip vectors with info about resolved states + + # get temperature and density vectors + for i in range(startsearch, len(lines)): + if re.match("^\s*C{0}-{2,}", lines[i]): + tmp = re.sub("\n*\s+", "\t", + "".join(lines[startsearch:i]).strip()) # replace unwanted chars + tmp = np.fromstring(tmp, sep="\t", dtype=float) # put into nunpy array + densities = tmp[:n_densities] # read density values + densities = 10**densities # convert from log10(cm^-3) to cm^-3 + densities = PerCm3ToPerM3.to(densities) # convert units from cm^-3 to m^-3 + temperatures = tmp[n_densities:] # read temperature values + temperatures = 10**temperatures # convert from log10(eV) to eV + startsearch = i + break + + # process rate data + rates = RecursiveDict() + + # get beginnig and end of requested rates data block and add it to xarray + blockrates_start = None + blockrates_stop = None + for i in range(startsearch, len(lines)): + + if re.match("^\s*C*-{2,}", lines[i]): # is it a rates block header? + + # is it a first data block found? + if not blockrates_start is None: + blockrates_stop = i # end of the requested block + + rates_table = re.sub("\n*\s+", "\t", + "".join(lines[ + blockrates_start:blockrates_stop]).strip()) # replace unwanted chars + rates_table = np.fromstring(rates_table, sep="\t", + dtype=float).reshape((n_temperatures, + n_densities)) # transform into an array + + # convert units from cm^-3 to m^-3 + rates_table = 10**rates_table + rates_table = Cm3ToM3.to(rates_table) + + rates[element][ion_charge]['ne'] = densities + rates[element][ion_charge]['te'] = temperatures + rates[element][ion_charge]['rates'] = np.swapaxes(rates_table, 0, 1) + + print() + print("density", densities.shape) + print("temperatures", temperatures.shape) + print("rates", rates_table.shape) + print() + + # if end of data block beak the loop or reassign start of data block for next stage + if re.match("^\s*C{1}-{2,}", lines[i]) or re.match("^\s*C{0,1}-{2,}", lines[i]) and \ + re.match("^\s*C\n", lines[i + 1]): + break + + z1_pos = re.search("Z1\s*=*\s*[0-9]+\s*", lines[i]).group() # get Z1 part + ion_charge = int(re.sub("Z1[\s*=]", "", z1_pos)) # remove Z1 to avoid getting 1 later + if not re.search("IGRD\s*=*\s*[0-9]+\s*", lines[i]) is None: # get the IGRD part + igrd_pos = re.search("IGRD\s*=*\s*[0-9]+\s*", lines[i]).group() # get the IGRD part + else: + igrd_pos = "No spec" + if not re.search("IPRT\s*=*\s*[0-9]+\s*", lines[i]) is None: + iptr_pos = re.search("IPRT\s*=*\s*[0-9]+\s*", lines[i]).group() # get the IPRT part + else: + iptr_pos = "No spec" + blockrates_start = i + 1 # if block start not known, check if we are at the right position + + return rates + + # def parse_adf11xxx(file_path, ion, ionisation): # # adf11_fh = open(file_path, 'r') diff --git a/cherab/openadas/rates/__init__.pxd b/cherab/openadas/rates/__init__.pxd index d6cd0fa..7cd28f9 100644 --- a/cherab/openadas/rates/__init__.pxd +++ b/cherab/openadas/rates/__init__.pxd @@ -19,3 +19,4 @@ from cherab.openadas.rates.beam cimport * from cherab.openadas.rates.cx cimport * from cherab.openadas.rates.pec cimport * +from cherab.openadas.rates.atomic cimport * diff --git a/cherab/openadas/rates/__init__.py b/cherab/openadas/rates/__init__.py index c2afe1b..09a2fce 100644 --- a/cherab/openadas/rates/__init__.py +++ b/cherab/openadas/rates/__init__.py @@ -18,4 +18,5 @@ from .beam import * from .cx import * -from .pec import * \ No newline at end of file +from .pec import * +from .atomic import * diff --git a/cherab/openadas/rates/atomic.pxd b/cherab/openadas/rates/atomic.pxd new file mode 100644 index 0000000..8767ad7 --- /dev/null +++ b/cherab/openadas/rates/atomic.pxd @@ -0,0 +1,47 @@ +# Copyright 2016-2018 Euratom +# Copyright 2016-2018 United Kingdom Atomic Energy Authority +# Copyright 2016-2018 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + +from cherab.core cimport IonisationRate as CoreIonisationRate +from cherab.core cimport RecombinationRate as CoreRecombinationRate +from cherab.core.math cimport Interpolate2DCubic + + +cdef class IonisationRate(CoreIonisationRate): + + cdef: + readonly dict raw_data + readonly double wavelength + readonly tuple density_range, temperature_range + Interpolate2DCubic _rate + + +cdef class NullImpactExcitationRate(CoreIonisationRate): + pass + + +cdef class RecombinationRate(CoreRecombinationRate): + + cdef: + readonly dict raw_data + readonly double wavelength + readonly tuple density_range, temperature_range + Interpolate2DCubic _rate + + +cdef class NullRecombinationRate(CoreRecombinationRate): + pass diff --git a/cherab/openadas/rates/atomic.pyx b/cherab/openadas/rates/atomic.pyx new file mode 100644 index 0000000..bbdf33c --- /dev/null +++ b/cherab/openadas/rates/atomic.pyx @@ -0,0 +1,99 @@ + +# Copyright 2016-2018 Euratom +# Copyright 2016-2018 United Kingdom Atomic Energy Authority +# Copyright 2016-2018 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + + +cdef class IonisationRate(CoreIonisationRate): + + def __init__(self, dict data, extrapolate=False): + """ + :param data: Dictionary containing rate data. + :param extrapolate: Enable extrapolation (default=False). + """ + + self.raw_data = data + + # unpack + ne = data['ne'] + te = data['te'] + rate = data['rate'] + + # store limits of data + self.density_range = ne.min(), ne.max() + self.temperature_range = te.min(), te.max() + + # interpolate rate + self._rate = Interpolate2DCubic( + ne, te, rate, extrapolate=extrapolate, extrapolation_type="quadratic" + ) + + cpdef double evaluate(self, double density, double temperature) except? -1e999: + + # prevent -ve values (possible if extrapolation enabled) + return max(0, self._rate.evaluate(density, temperature)) + + +cdef class NullIonisationRate(CoreIonisationRate): + """ + A PEC rate that always returns zero. + Needed for use cases where the required atomic data is missing. + """ + + cpdef double evaluate(self, double density, double temperature) except? -1e999: + return 0.0 + + +cdef class RecombinationRate(CoreRecombinationRate): + + def __init__(self, dict data, extrapolate=False): + """ + :param data: Dictionary containing rate data. + :param extrapolate: Enable extrapolation (default=False). + """ + + self.raw_data = data + + # unpack + ne = data['ne'] + te = data['te'] + rate = data['rate'] + + # store limits of data + self.density_range = ne.min(), ne.max() + self.temperature_range = te.min(), te.max() + + # interpolate rate + self._rate = Interpolate2DCubic( + ne, te, rate, extrapolate=extrapolate, extrapolation_type="quadratic" + ) + + + cpdef double evaluate(self, double density, double temperature) except? -1e999: + + # prevent -ve values (possible if extrapolation enabled) + return max(0, self._rate.evaluate(density, temperature)) + + +cdef class NullRecombinationRate(CoreRecombinationRate): + """ + A PEC rate that always returns zero. + Needed for use cases where the required atomic data is missing. + """ + + cpdef double evaluate(self, double density, double temperature) except? -1e999: + return 0.0 diff --git a/cherab/openadas/rates/beam.pxd b/cherab/openadas/rates/beam.pxd index 8e9bece..ba27aaa 100644 --- a/cherab/openadas/rates/beam.pxd +++ b/cherab/openadas/rates/beam.pxd @@ -18,7 +18,7 @@ from cherab.core cimport BeamStoppingRate as CoreBeamStoppingRate from cherab.core cimport BeamPopulationRate as CoreBeamPopulationRate -from cherab.core cimport BeamEmissionRate as CoreBeamEmissionRate +from cherab.core cimport BeamEmissionPEC as CoreBeamEmissionPEC from cherab.core.math cimport Function1D, Function2D @@ -52,7 +52,7 @@ cdef class NullBeamPopulationRate(CoreBeamPopulationRate): pass -cdef class BeamEmissionRate(CoreBeamEmissionRate): +cdef class BeamEmissionPEC(CoreBeamEmissionPEC): cdef readonly: dict raw_data @@ -63,5 +63,5 @@ cdef class BeamEmissionRate(CoreBeamEmissionRate): tuple temperature_range -cdef class NullBeamEmissionRate(CoreBeamEmissionRate): +cdef class NullBeamEmissionPEC(CoreBeamEmissionPEC): pass diff --git a/cherab/openadas/rates/beam.pyx b/cherab/openadas/rates/beam.pyx index 154a16f..8f7f0df 100644 --- a/cherab/openadas/rates/beam.pyx +++ b/cherab/openadas/rates/beam.pyx @@ -151,7 +151,7 @@ cdef class NullBeamPopulationRate(CoreBeamPopulationRate): return 0.0 -cdef class BeamEmissionRate(CoreBeamEmissionRate): +cdef class BeamEmissionPEC(CoreBeamEmissionPEC): """ The beam emission coefficient interpolation class. @@ -206,7 +206,7 @@ cdef class BeamEmissionRate(CoreBeamEmissionRate): return val -cdef class NullBeamEmissionRate(CoreBeamEmissionRate): +cdef class NullBeamEmissionPEC(CoreBeamEmissionPEC): """ A beam rate that always returns zero. Needed for use cases where the required atomic data is missing. diff --git a/cherab/openadas/rates/cx.pxd b/cherab/openadas/rates/cx.pxd index a2680d6..47fdf75 100644 --- a/cherab/openadas/rates/cx.pxd +++ b/cherab/openadas/rates/cx.pxd @@ -16,11 +16,11 @@ # See the Licence for the specific language governing permissions and limitations # under the Licence. -from cherab.core cimport BeamCXRate as CoreBeamCXRate +from cherab.core cimport BeamCXPEC as CoreBeamCXPEC from cherab.core.math cimport Function1D, Function2D -cdef class BeamCXRate(CoreBeamCXRate): +cdef class BeamCXPEC(CoreBeamCXPEC): cdef readonly: dict raw_data @@ -34,5 +34,5 @@ cdef class BeamCXRate(CoreBeamCXRate): readonly tuple b_field_range -cdef class NullBeamCXRate(CoreBeamCXRate): +cdef class NullBeamCXPEC(CoreBeamCXPEC): pass diff --git a/cherab/openadas/rates/cx.pyx b/cherab/openadas/rates/cx.pyx index 3714cee..ece84cf 100644 --- a/cherab/openadas/rates/cx.pyx +++ b/cherab/openadas/rates/cx.pyx @@ -22,7 +22,7 @@ cimport cython from cherab.core.math cimport Interpolate1DCubic -cdef class BeamCXRate(CoreBeamCXRate): +cdef class BeamCXPEC(CoreBeamCXPEC): """ The effective cx rate interpolation class. @@ -106,7 +106,7 @@ cdef class BeamCXRate(CoreBeamCXRate): return rate -cdef class NullBeamCXRate(CoreBeamCXRate): +cdef class NullBeamCXPEC(CoreBeamCXPEC): """ A beam CX rate that always returns zero. Needed for use cases where the required atomic data is missing. diff --git a/cherab/openadas/rates/pec.pxd b/cherab/openadas/rates/pec.pxd index 05b6d59..435b172 100644 --- a/cherab/openadas/rates/pec.pxd +++ b/cherab/openadas/rates/pec.pxd @@ -16,13 +16,13 @@ # See the Licence for the specific language governing permissions and limitations # under the Licence. -from cherab.core cimport ImpactExcitationRate as CoreImpactExcitationRate -from cherab.core cimport RecombinationRate as CoreRecombinationRate -from cherab.core cimport ThermalCXRate as CoreThermalCXRate +from cherab.core cimport ImpactExcitationPEC as CoreImpactExcitationPEC +from cherab.core cimport RecombinationPEC as CoreRecombinationPEC +from cherab.core cimport ThermalCXPEC as CoreThermalCXPEC from cherab.core.math cimport Interpolate2DCubic -cdef class ImpactExcitationRate(CoreImpactExcitationRate): +cdef class ImpactExcitationPEC(CoreImpactExcitationPEC): cdef: readonly dict raw_data @@ -31,11 +31,11 @@ cdef class ImpactExcitationRate(CoreImpactExcitationRate): Interpolate2DCubic _rate -cdef class NullImpactExcitationRate(CoreImpactExcitationRate): +cdef class NullImpactExcitationPEC(CoreImpactExcitationPEC): pass -cdef class RecombinationRate(CoreRecombinationRate): +cdef class RecombinationPEC(CoreRecombinationPEC): cdef: readonly dict raw_data @@ -44,7 +44,7 @@ cdef class RecombinationRate(CoreRecombinationRate): Interpolate2DCubic _rate -cdef class NullRecombinationRate(CoreRecombinationRate): +cdef class NullRecombinationPEC(CoreRecombinationPEC): pass diff --git a/cherab/openadas/rates/pec.pyx b/cherab/openadas/rates/pec.pyx index 0f8c396..c2f3bc4 100644 --- a/cherab/openadas/rates/pec.pyx +++ b/cherab/openadas/rates/pec.pyx @@ -23,7 +23,7 @@ import matplotlib.pyplot as plt from cherab.core.utility.conversion import PhotonToJ -cdef class ImpactExcitationRate(CoreImpactExcitationRate): +cdef class ImpactExcitationPEC(CoreImpactExcitationPEC): def __init__(self, double wavelength, dict data, extrapolate=False): """ @@ -58,7 +58,7 @@ cdef class ImpactExcitationRate(CoreImpactExcitationRate): return max(0, self._rate.evaluate(density, temperature)) -cdef class NullImpactExcitationRate(CoreImpactExcitationRate): +cdef class NullImpactExcitationPEC(CoreImpactExcitationPEC): """ A PEC rate that always returns zero. Needed for use cases where the required atomic data is missing. @@ -68,7 +68,7 @@ cdef class NullImpactExcitationRate(CoreImpactExcitationRate): return 0.0 -cdef class RecombinationRate(CoreRecombinationRate): +cdef class RecombinationPEC(CoreRecombinationPEC): def __init__(self, double wavelength, dict data, extrapolate=False): """ @@ -104,7 +104,7 @@ cdef class RecombinationRate(CoreRecombinationRate): return max(0, self._rate.evaluate(density, temperature)) -cdef class NullRecombinationRate(CoreRecombinationRate): +cdef class NullRecombinationPEC(CoreRecombinationPEC): """ A PEC rate that always returns zero. Needed for use cases where the required atomic data is missing. diff --git a/cherab/openadas/repository/__init__.py b/cherab/openadas/repository/__init__.py index af84b36..e589658 100644 --- a/cherab/openadas/repository/__init__.py +++ b/cherab/openadas/repository/__init__.py @@ -18,6 +18,7 @@ from .beam import * from .pec import * +from .atomic import * from .wavelength import * from .utility import DEFAULT_REPOSITORY_PATH from .create import populate diff --git a/cherab/openadas/repository/atomic.py b/cherab/openadas/repository/atomic.py new file mode 100644 index 0000000..be9b6a9 --- /dev/null +++ b/cherab/openadas/repository/atomic.py @@ -0,0 +1,197 @@ +# Copyright 2016-2018 Euratom +# Copyright 2016-2018 United Kingdom Atomic Energy Authority +# Copyright 2016-2018 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas +# +# Licensed under the EUPL, Version 1.1 or – as soon they will be approved by the +# European Commission - subsequent versions of the EUPL (the "Licence"); +# You may not use this work except in compliance with the Licence. +# You may obtain a copy of the Licence at: +# +# https://joinup.ec.europa.eu/software/page/eupl5 +# +# Unless required by applicable law or agreed to in writing, software distributed +# under the Licence is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR +# CONDITIONS OF ANY KIND, either express or implied. +# +# See the Licence for the specific language governing permissions and limitations +# under the Licence. + + +import os +import json +import numpy as np + +from cherab.core.atomic import Element +from cherab.core.utility import RecursiveDict +from .utility import DEFAULT_REPOSITORY_PATH, valid_ionisation + + +def add_ionisation_rate(species, ionisation, rate, repository_path=None): + """ + Adds a single ionisation rate to the repository. + + If adding multiple rates, consider using the update_ionisation_rates() + function instead. The update function avoids repeatedly opening and closing + the rate files. + + :param repository_path: + :return: + """ + + update_ionisation_rates({ + species: { + ionisation: rate + } + }, repository_path) + + +def update_ionisation_rates(rates, repository_path=None): + """ + Ionisation rate file structure + + /ionisation/.json + + File contains multiple rates, indexed by ion ionisation. + """ + + repository_path = repository_path or DEFAULT_REPOSITORY_PATH + + for species, rate_data in rates.items(): + + # sanitise and validate arguments + if not isinstance(species, Element): + raise TypeError('The species must be an Element object.') + + path = os.path.join(repository_path, 'ionisation/{}.json'.format(species.symbol.lower())) + + _update_and_write_adf11(species, rate_data, path) + + +def add_recombination_rate(species, ionisation, rate, repository_path=None): + """ + Adds a single recombination rate to the repository. + + If adding multiple rates, consider using the update_recombination_rates() + function instead. The update function avoids repeatedly opening and closing + the rate files. + + :param repository_path: + :return: + """ + + update_recombination_rates({ + species: { + ionisation: rate + } + }, repository_path) + + +def update_recombination_rates(rates, repository_path=None): + """ + Ionisation rate file structure + + /recombination/.json + + File contains multiple rates, indexed by ion ionisation. + """ + + repository_path = repository_path or DEFAULT_REPOSITORY_PATH + + for species, rate_data in rates.items(): + + # sanitise and validate arguments + if not isinstance(species, Element): + raise TypeError('The species must be an Element object.') + + path = os.path.join(repository_path, 'recombination/{}.json'.format(species.symbol.lower())) + + _update_and_write_adf11(species, rate_data, path) + + +def _update_and_write_adf11(species, rate_data, path): + + # read in any existing rates + try: + with open(path, 'r') as f: + content = RecursiveDict.from_dict(json.load(f)) + except FileNotFoundError: + content = RecursiveDict() + + for ionisation, rates in rate_data.items(): + + if not valid_ionisation(species, ionisation): + raise ValueError('Ionisation level is larger than the number of protons in the specified species.') + + # sanitise and validate rate data + te = np.array(rates['te'], np.float64) + ne = np.array(rates['ne'], np.float64) + rate_table = np.array(rates['rates'], np.float64) + + if ne.ndim != 1: + raise ValueError('Density array must be a 1D array.') + + if te.ndim != 1: + raise ValueError('Temperature array must be a 1D array.') + + if (ne.shape[0], te.shape[0]) != rate_table.shape: + raise ValueError('Electron temperature, density and rate data arrays have inconsistent sizes.') + + # update file content with new rate + content[ionisation] = { + 'te': te.tolist(), + 'ne': ne.tolist(), + 'rate': rate_table.tolist(), + } + + # create directory structure if missing + directory = os.path.dirname(path) + if not os.path.isdir(directory): + os.makedirs(directory) + + # write new data + with open(path, 'w') as f: + json.dump(content, f, indent=2, sort_keys=True) + + +def get_ionisation_rate(element, ionisation, repository_path=None): + + repository_path = repository_path or DEFAULT_REPOSITORY_PATH + + path = os.path.join(repository_path, 'ionisation/{}.json'.format(element.symbol.lower())) + try: + with open(path, 'r') as f: + content = json.load(f) + d = content[str(ionisation)] + except (FileNotFoundError, KeyError): + print() + print(path) + raise RuntimeError('Requested ionisation rate (element={}, ionisation={})' + ' is not available.'.format(element.symbol, ionisation)) + + # convert to numpy arrays + d['ne'] = np.array(d['ne'], np.float64) + d['te'] = np.array(d['te'], np.float64) + d['rate'] = np.array(d['rate'], np.float64) + + return d + + +def get_recombination_rate(element, ionisation, repository_path=None): + + repository_path = repository_path or DEFAULT_REPOSITORY_PATH + + path = os.path.join(repository_path, 'recombination/{}.json'.format(element.symbol.lower())) + try: + with open(path, 'r') as f: + content = json.load(f) + d = content[str(ionisation)] + except (FileNotFoundError, KeyError): + raise RuntimeError('Requested recombination rate (element={}, ionisation={})' + ' is not available.'.format(element.symbol, ionisation)) + + # convert to numpy arrays + d['ne'] = np.array(d['ne'], np.float64) + d['te'] = np.array(d['te'], np.float64) + d['rate'] = np.array(d['rate'], np.float64) + + return d diff --git a/cherab/openadas/repository/create.py b/cherab/openadas/repository/create.py index f42a31d..9c1d1a8 100644 --- a/cherab/openadas/repository/create.py +++ b/cherab/openadas/repository/create.py @@ -37,6 +37,14 @@ def populate(download=True, repository_path=None, adas_path=None): # install a common selection of open adas files rates = { + 'adf11scd': ( + (hydrogen, 'adf11/scd12/scd12_h.dat'), + (carbon, 'adf11/scd96/scd96_c.dat'), + ), + 'adf11acd': ( + (hydrogen, 'adf11/acd12/acd12_h.dat'), + (carbon, 'adf11/acd96/acd96_c.dat') + ), # 'adf11plt': ( # (hydrogen, "adf11/plt12/plt12_h.dat"), # (helium, "adf11/plt96/plt96_he.dat"),