diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000..8c7cc74 --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1,4 @@ +include README.md CHANGELOG.md LICENSE.txt MANIFEST.in setup.py .gitignore +recursive-include cherab *.py *.pyx *.pxd *.json +recursive-include demos *.py *.pyx *.pxd + diff --git a/cherab/openadas/ADF11Xarray.py b/cherab/openadas/ADF11Xarray.py new file mode 100644 index 0000000..8fb99a4 --- /dev/null +++ b/cherab/openadas/ADF11Xarray.py @@ -0,0 +1,442 @@ +import numpy as np +import re +from scipy import interpolate +from math import ceil +from scipy.optimize import nnls +import xarray as xr +from os.path import expanduser +from scipy import interpolate +import scipy.constants as const +import matplotlib.pyplot as plt +from sympy import Matrix +import xarray as xr +from scipy.linalg import solve + +def rates_adf11(file, ionicstate, rate_ne, rate_te, igrd=None, iprt=None): + # TODO:proper handling and naming for resolved cases + # TODO: degree of interpolation, now set to 2 + # TODO: replace current interpolation by interpolation used in cherab + """ + Reads contents of open adas adf11 files + + :param file: path to the adf11 file + :param ionicstate: recombining ion charge to get the rates for (i.e. for ionization it will be Z+1) + :param rate_ne: electron density to calculate rate for + :param rate_te: electron temperature to calculate rate for + :param igrd: something for resolved states, to be filled in later + :param iprt: something for resolved states, to be filled in later + :return: temperature, density, rates as numpy array + """ + with open(file, "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]) + density_n = int(tmp[1]) + temperature_n = int(tmp[2]) + z_min = int(tmp[3]) + z_max = int(tmp[4]) + element = tmp[5] + projectname = tmp[6] + + # 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("C{0}-+[^-^\n]", 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 + density = tmp[:density_n] # read density values + temperature = tmp[density_n:] # read temperature values + startsearch = i + break + + # get beginnig and end of requested rates data block + blockrates_start = None + blockrates_stop = None + + for i in range(startsearch, len(lines)): + if re.match("C{0}-+[^-^\n]", lines[i]): # is it a rates block header? + if not blockrates_start: # if block start not known, check if we are at the right position + z1_pos = re.search("Z1\s*=*\s*[0-9]+\s*", lines[i]).group() # get Z1 part + z1_pos = re.sub("Z1", "", z1_pos) # remove Z1 to avoid getting 1 later + igrd_pos = re.search("IGRD\s*=*\s*[0-9]+\s*", lines[i]).group() # get the IGRD part + iptr_pos = re.search("IPRT\s*=*\s*[0-9]+\s*", lines[i]).group() # get the IPRT part + # is it the requested block? + if not igrd and not iprt or\ + int(re.search("[0-9]+", z1_pos).group()) == ionicstate and \ + int(re.search("[0-9]+", igrd_pos).group()) == igrd and \ + int(re.search("[0-9]+", iptr_pos).group()) == iprt: + blockrates_start = i+1 + else: + blockrates_stop = i # end of the requested block + break + + # get rates into an 2D array and interpolate rates for requested electron density and temperature + if blockrates_start and blockrates_stop: + 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((temperature_n, + density_n)) # transform into an array + interp = interpolate.interp2d(density, temperature, rates_table, kind='cubic') # interpolate + + rates = np.power(10, interp(np.log10(rate_ne), np.log10(rate_te))) # calculate requested rates + return {'DENS': rate_ne, 'TE': rate_te, 'rate': rates, + "[DENS]": "cm^-3", "[Te]": "eV", "[rate]": "cm^3/s"} + else: + print("wrong charge or indexes") + + +class ADF11(): + + def __init__(self, element, year): + + self.element = element + self.year = year + self.coef_ion = None + self.coef_recom = None + self.coef_cxh= None + #self.coef_linepower = None + + + self.readfiles() + + def readfiles(self): + home = expanduser("~") + #ionisation ...scd, recombination ... acd, cx with H ... ccd + filepath = home+"/Python/adas/adf11/{0}{1}/{0}{1}_{2}.dat" # coef, year, element + + try: + self.coef_ion = self.readdata(filepath.format("scd", self.year, self.element)) + except FileNotFoundError: + pass + try: + self.coef_recom = self.readdata(filepath.format("acd", self.year, self.element)) + except FileNotFoundError: + pass + + try: + self.coef_cxh = self.readdata(filepath.format("ccd", self.year, self.element)) + except FileNotFoundError: + pass + + def readdata(self, filepath): + + with open(filepath, "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]) + density_n = int(tmp[1]) + temperature_n = int(tmp[2]) + z_min = int(tmp[3]) + z_max = int(tmp[4]) + element = tmp[5] + projectname = tmp[6] + + # 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 + density = tmp[:density_n] # read density values + temperature = tmp[density_n:] # read temperature values + startsearch = i + break + + # 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((temperature_n, + density_n)) # transform into an array + rates_table = np.expand_dims(rates_table,0) + tmp = xr.DataArray(rates_table, + dims=["recomioncharge", "Te", "ne"], + coords={"Te": temperature, + "ne": density, + "recomioncharge": [ioncharge]} + ) + + try: + ret = xr.concat((ret, tmp.copy()), dim="recomioncharge") + except UnboundLocalError: + ret = tmp.copy() + + # 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 + ioncharge = 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 ret + + def interpvalues(self, data, charge,te,ne): + """ + interpolates adf11 coefficients + :param data: ionization, recombination, cx, ..., data + :param charge: recombining ion charge. In case of ionzation it is the ionic stage charge +1 + :param te: + :param ne: + :return: + """ + #mesh_te, mesh_ne = np.meshgrid(data.Te.data, data.ne.data) + interp = interpolate.interp2d(data.Te.data, data.ne.data, data.sel(recomioncharge=charge).data.T, kind="linear") + values = interp(np.log10(te), np.log10(ne)) + return np.power(10, values) + + def ionbalance(self,te, ne, nh = np.array([None])): + numstates = self.coef_recom.recomioncharge.shape[0]+1 + fracabun = 0 + + for j in range(te.shape[0]): + matbal = np.zeros((numstates, numstates)) + + matbal[0,0] -= self.interpvalues(self.coef_ion,1, te[j], ne[j]) + matbal[0,1] += self.interpvalues(self.coef_recom,1,te[j],ne[j]) + matbal[-1,-1] -= self.interpvalues(self.coef_recom, self.coef_recom.recomioncharge.max(), te[j], ne[j]) + matbal[-1, -2] += self.interpvalues(self.coef_ion, self.coef_ion.recomioncharge.max(), te[j], ne[j]) + + if not nh[0] == None: + try: + matbal[0, 1] += self.interpvalues(self.coef_cxh, 1, te[j], ne[j]) * nh[j] / ne[j] + matbal[-1, -1] -= self.interpvalues(self.coef_cxh, self.coef_recom.recomioncharge.max(), + te[j], ne[j]) * nh[j] / ne[j] + except AttributeError: + pass + + for i in range(1,numstates-1): + matbal[i,i-1] += self.interpvalues(self.coef_ion,i, te[j], ne[j]) + matbal[i,i] -= (self.interpvalues(self.coef_ion,i+1, te[j], ne[j]) + + self.interpvalues(self.coef_recom,i,te[j],ne[j])) + matbal[i, i+1] += self.interpvalues(self.coef_recom,i+1,te[j],ne[j]) + if not nh[0] == None: + try: + matbal[i, i ] -= self.interpvalues(self.coef_cxh,i, te[j], ne[j])* nh[j]/ne[j] + matbal[i, i+1] += self.interpvalues(self.coef_cxh,i+1, te[j], ne[j])* nh[j]/ne[j] + except AttributeError: + pass + + matbal = np.concatenate((matbal, np.ones((1, matbal.shape[1]))), axis=0) + matbal = np.concatenate((matbal, np.zeros((matbal.shape[0], 1))), axis=1) + solution = np.zeros((matbal.shape[0])) + solution[-1] = 1 + tmp = np.linalg.lstsq(matbal,solution)[0][0:-1] + tmp = np.expand_dims(tmp,axis=1) + tmp = xr.DataArray(tmp, + coords={"charge": np.arange(numstates), + "Te":[te[j]], + "ne":("Te",[ne[j]])}, + dims=["charge","Te"]) + try: + ret = xr.concat((ret,tmp.copy()), dim="Te") + except UnboundLocalError: + ret = tmp.copy() + + return ret + + def null(self, A, eps=1e-20): + u, s, vh = np.linalg.svd(A) + null_space = vh[np.argmin(s),:] + return null_space.T + +def speed(temp,amu): + vp= np.sqrt(2*const.k*temp/(amu*const.atomic_mass)) + return vp + +def meanfreepath(): + home = expanduser("~") + amu = {} + amu["li"] = 6.94 + amu["h"] = 1.08 + amu["he"] = 4.0 + amu["c"] = 12.01 + amu["n"] = 14.07 + + elem = "n" + filepath = home+"/Python/adas/adf11/scd96/scd96_{0}.dat".format(elem) + data = ADF11.readData(filepath) + + if True: + compass = np.loadtxt("livap01_profile.dat", delimiter=",") + #compass =np.loadtxt("livap02_profile.dat", delimiter=",") + + step = np.diff(compass[:,0]*1e-3)[0] + te = compass[:, 2] + ne = compass[:, 1]*1e-6 + + alpha = np.zeros_like(te) + mfp = np.zeros_like(te) + particles_loss = np.zeros_like(te) + particles = np.zeros_like(te) + + tli = 25 + spd = speed(tli+273, amu[elem]) + + i = 0 + alpha[i] = data.interpvalues(1, te[i], ne[i]) + mfp[i] = 1 / (alpha[i] * ne[i] / spd) + particles_loss[i] = 1 - np.exp(-1 * step / mfp[i]) + particles[i] = np.exp(-1 * step / mfp[i]) + + for i in range(1,te.shape[0]): + alpha[i] = data.interpvalues( 1, te[i], ne[i]) + mfp[i] = 1/(alpha[i]*ne[i]/spd) + particles_loss[i] = particles[i-1]*(1 - np.exp(-1 * step / mfp[i])) + particles[i] = particles[i-1]*np.exp(-1 * step / mfp[i]) + + + loss_mean = np.sum(np.multiply(particles_loss,compass[:,0]))/np.sum(particles_loss) + + + fig, ax = plt.subplots() + ax.plot(compass[:,0], particles) + + ax.set_title("COMPASS-U, lithium vapour temperature {0:d} K".format(tli+273)) + ax.set_xlabel("divertor distance d [mm]") + ax.set_ylabel("Fractional Density $\\frac{n(d)}{n(0)}$ [a.u.]") + + fig, ax = plt.subplots() + ax.plot(compass[:, 0], particles_loss) + ax.axvline(loss_mean, color="r") + plt.text(loss_mean+0.1, 0.015,"E(Particle loss)={0:2.3f} mm".format(loss_mean)) + ax.set_title("COMPASS-U, lithium vapour temperature {0:d} K".format(tli + 273)) + ax.set_xlabel("divertor distance d [mm]") + ax.set_ylabel("Fractional Particle Loss $\\frac{dn(d)}{n(0)}$ [a.u.]") + + + figpl, ax = plt.subplots() + twinx = ax.twinx() + ax.plot(compass[:,0],compass[:,1], label ="ne") + twinx.plot(compass[:, 0], compass[:, 2],"r", label="Te") + ax.plot([],[],"r", label="Te") + ax.set_title("COMPASS-U") + ax.set_ylabel("ne [m$^{-3}$]") + twinx.set_ylabel("Te [eV]") + ax.set_xlabel("divertor distance d [mm]") + ax.legend() + + if False: + te = np.linspace(0.1, 5, 100) + ne = np.linspace(1e14, 1e15, 100) + alpha = np.zeros((te.shape[0], ne.shape[0])) + mfp = np.zeros_like(alpha) + spd = speed(700, amu[elem]) + + for i in range(te.shape[0]): + for j in range(ne.shape[0]): + interpolated = interpvalues(data, 1, te[i], ne[j]) + alpha[j,i] = interpolated + mfp[j,i] = 1/(interpolated*ne[j]/spd) + + gig, ax = plt.subplots() + cf2 = ax.contour(te, ne * 1e8, np.log10(mfp), 5, colors = "w",ls="-") + cf1 = ax.contourf(te, ne*1e8, np.log10(mfp), 400) + cbar = gig.colorbar(cf1,ax=ax,extend="max") + cf2.clabel() + cbar.set_label("log10(mfp) [m]") + ax.set_xlabel("Te [eV]") + ax.set_ylabel("ne [m$^{-3}$]") + ax.set_title("Mean Free Path {}".format(elem)) + + #te = np.linspace(1, 40, 100) + nec = np.ones_like(te)*1e13 + #alpha = interpvalues(data,1,te,ne)[0,:] + + asd, ax = plt.subplots() + (data.sel(ne = 13, charge=1)+13).plot.line(ax=ax) + (data.sel(ne=14, charge=1)+14).plot.line(ax=ax) + ax.plot(np.log10(te), np.log10(alpha[0,:])+13,"--") + ax.plot(np.log10(te), np.log10(alpha[-1, :]) + 14, "--") + ax.plot(np.log10(te), np.log10(spd/(mfp[-1,:]*1e14)), "--") + + asd, ax = plt.subplots() + ax.plot(te, np.log10(mfp[0, :]), "-") + ax.plot(te, np.log10(mfp[-1,:]), "--") + +if __name__ == "__main__": + te = np.arange(1,500,0.5) + #te = np.concatenate((te,np.arange(20,100,0.5) )) + #te = np.concatenate((te, np.arange(100, 400, 1))) + + ne = np.ones_like(te) + col = ["b","r","g","c","m","y","k"] + carbon = ADF11("c", 96) + boron = ADF11("b", 89) + helium = ADF11("he", 96) + + balc_h = carbon.ionbalance(te,ne*1e13, ne*1e13*np.exp(-1*te/100)/1000) + balc = carbon.ionbalance(te,ne*1e13) + balb_h = boron.ionbalance(te,ne*1e13, ne*1e13*np.exp(-1*te/100)/1000) + balb = boron.ionbalance(te,ne*1e13) + + ionbal, ax = plt.subplots() + for i in balc.charge.data: + balc.sel(charge = i).plot.line(ax=ax, label="{0}+".format(i), color = col[i]) + balc_h.sel(charge=i).plot.line(ax=ax, ls = "--", color = col[i]) + plt.plot([],[],"k-", label = "No TCX") + plt.plot([], [], "k--", label="TCX") + ax.set_title("Carbon Ionisation Balance") + ax.set_ylabel("Fractional Abundance [a.u.]") + ax.legend() + + + ionbal, ax = plt.subplots() + for i in balb.charge.data: + balb.sel(charge = i).plot.line(ax=ax, label="{0}+".format(i), color = col[i]) + balb_h.sel(charge=i).plot.line(ax=ax,ls = "--", color = col[i]) + plt.plot([],[],"k-", label = "No TCX") + plt.plot([], [], "k--", label="TCX") + ax.set_title("Boron Ionisation Balance") + ax.set_ylabel("Fractional Abundance [a.u.]") + ax.legend() + + ionbal, ax = plt.subplots() + for i in balb.charge.data: + balb.sel(charge = i).plot.line(ax=ax, label="{0}".format(i),ls = "--", color = col[i]) + for i in balc.charge.data: + balc.sel(charge = i).plot.line(ax=ax, label="{0}".format(i), color = col[i]) + + ionbal, ax = plt.subplots() + for i in balb.charge.data: + balb_h.sel(charge = i).plot.line(ax=ax, label="{0}".format(i),ls = "--", color = col[i]) + for i in balc.charge.data: + balc_h.sel(charge = i).plot.line(ax=ax, label="{0}".format(i), color = col[i]) + + + figs, ax = plt.subplots() + ax.plot(te, np.log10(np.exp(-1*te/100)/1000)) \ No newline at end of file diff --git a/cherab/openadas/config.py b/cherab/openadas/config.py index d380409..f7e2399 100644 --- a/cherab/openadas/config.py +++ b/cherab/openadas/config.py @@ -136,7 +136,7 @@ def _check_for_adf_file(open_adas_root_path, relative_adf_file_path, download_pa for element, ionisation, relative_adf_path, download_path in DEFAULT_ADF15_FILES: print(relative_adf_path) _check_for_adf_file(ROOT_DATA_PATH, relative_adf_path, download_path) - add_adf15_to_atomic_data(default, element, ionisation, os.path.join(ROOT_DATA_PATH, relative_adf_path)) + default = add_adf15_to_atomic_data(default, element, ionisation, os.path.join(ROOT_DATA_PATH, relative_adf_path)) # TODO: there are concerns about the accuracy of this data # Data from NIST Atomic Spectra Database - http://www.nist.gov/pml/data/asd.cfm diff --git a/cherab/openadas/library/__init__.py b/cherab/openadas/library/__init__.py new file mode 100644 index 0000000..1550aec --- /dev/null +++ b/cherab/openadas/library/__init__.py @@ -0,0 +1,18 @@ +# Copyright 2014-2017 United Kingdom Atomic Energy Authority +# +# 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 .library import ADF11_PLT_FILES, ADF11_PRB_FILES diff --git a/cherab/openadas/library/adf11.json b/cherab/openadas/library/adf11.json new file mode 100644 index 0000000..1b79325 --- /dev/null +++ b/cherab/openadas/library/adf11.json @@ -0,0 +1,87 @@ +{ + "Description": "Recommended set of ADF11 files for PLT coefficients", + "adf11_plt_files": { + "H": { + "ADAS_Path": "adf11/plt12/plt12_h.dat", + "Download_URL": "http://open.adas.ac.uk/download/adf11/plt12/plt12_h.dat" + }, + "He": { + "ADAS_Path": "adf11/plt96/plt96_he.dat", + "Download_URL": "http://open.adas.ac.uk/download/adf11/plt96/plt96_he.dat" + }, + "Li": { + "ADAS_Path": "adf11/plt96/plt96_li.dat", + "Download_URL": "http://open.adas.ac.uk/download/adf11/plt96/plt96_li.dat" + }, + "Be": { + "ADAS_Path": "adf11/plt96/plt96_be.dat", + "Download_URL": "http://open.adas.ac.uk/download/adf11/plt96/plt96_be.dat" + }, + "C": { + "ADAS_Path": "adf11/plt96/plt96_c.dat", + "Download_URL": "http://open.adas.ac.uk/download/adf11/plt96/plt96_c.dat" + }, + "N": { + "ADAS_Path": "adf11/plt96/plt96_n.dat", + "Download_URL": "http://open.adas.ac.uk/download/adf11/plt96/plt96_n.dat" + }, + "O": { + "ADAS_Path": "adf11/plt96/plt96_o.dat", + "Download_URL": "http://open.adas.ac.uk/download/adf11/plt96/plt96_o.dat" + }, + "Ne": { + "ADAS_Path": "adf11/plt96/plt96_ne.dat", + "Download_URL": "http://open.adas.ac.uk/download/adf11/plt96/plt96_ne.dat" + }, + "Ar": { + "ADAS_Path": "adf11/plt40/plt40_ar.dat", + "Download_URL": "http://open.adas.ac.uk/download/adf11/plt40/plt40_ar.dat" + }, + "Kr": { + "ADAS_Path": "adf11/plt89/plt89_kr.dat", + "Download_URL": "http://open.adas.ac.uk/download/adf11/plt89/plt89_kr.dat" + } + }, + "adf11_prb_files": { + "H": { + "ADAS_Path": "adf11/prb12/prb12_h.dat", + "Download_URL": "http://open.adas.ac.uk/download/adf11/prb12/prb12_h.dat" + }, + "He": { + "ADAS_Path": "adf11/prb96/prb96_he.dat", + "Download_URL": "http://open.adas.ac.uk/download/adf11/prb96/prb96_he.dat" + }, + "Li": { + "ADAS_Path": "adf11/prb96/prb96_li.dat", + "Download_URL": "http://open.adas.ac.uk/download/adf11/prb96/prb96_li.dat" + }, + "Be": { + "ADAS_Path": "adf11/prb96/prb96_be.dat", + "Download_URL": "http://open.adas.ac.uk/download/adf11/prb96/prb96_be.dat" + }, + "C": { + "ADAS_Path": "adf11/prb96/prb96_c.dat", + "Download_URL": "http://open.adas.ac.uk/download/adf11/prb96/prb96_c.dat" + }, + "N": { + "ADAS_Path": "adf11/prb96/prb96_n.dat", + "Download_URL": "http://open.adas.ac.uk/download/adf11/prb96/prb96_n.dat" + }, + "O": { + "ADAS_Path": "adf11/prb96/prb96_o.dat", + "Download_URL": "http://open.adas.ac.uk/download/adf11/prb96/prb96_o.dat" + }, + "Ne": { + "ADAS_Path": "adf11/prb96/prb96_ne.dat", + "Download_URL": "http://open.adas.ac.uk/download/adf11/prb96/prb96_ne.dat" + }, + "Ar": { + "ADAS_Path": "adf11/prb89/prb89_ar.dat", + "Download_URL": "http://open.adas.ac.uk/download/adf11/prb89/prb89_ar.dat" + }, + "Kr": { + "ADAS_Path": "adf11/prb89/prb89_kr.dat", + "Download_URL": "http://open.adas.ac.uk/download/adf11/prb89/prb89_kr.dat" + } + } +} \ No newline at end of file diff --git a/cherab/openadas/library/library.py b/cherab/openadas/library/library.py new file mode 100644 index 0000000..b8f89f0 --- /dev/null +++ b/cherab/openadas/library/library.py @@ -0,0 +1,29 @@ +# Copyright 2014-2017 United Kingdom Atomic Energy Authority +# +# 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 + + +LIBRARY_PATH = os.path.split(__file__)[0] + +ADF11_FILES = json.load(open(os.path.join(LIBRARY_PATH, 'adf11.json'))) +ADF11_PLT_FILES = ADF11_FILES['adf11_plt_files'] +ADF11_PRB_FILES = ADF11_FILES['adf11_prb_files'] + + + diff --git a/cherab/openadas/openadas.py b/cherab/openadas/openadas.py index 349efb6..6c27b8a 100644 --- a/cherab/openadas/openadas.py +++ b/cherab/openadas/openadas.py @@ -15,14 +15,18 @@ # under the Licence. import os +import urllib.request -from cherab.core import AtomicData, Isotope +from cherab.core import AtomicData from cherab.core.atomic.elements import Isotope -from pkg_resources import resource_filename +from cherab.core.utility.recursivedict import RecursiveDict -from cherab.openadas.read import adf12, adf15, adf21, adf22 +from cherab.openadas.library import * +from cherab.openadas.read import adf11, adf12, adf15, adf21, adf22 +from cherab.openadas.read.adf15 import add_adf15_to_atomic_data from . import config from .rates import * +from cherab.openadas.rates.radiated_power import StageResolvedRadiation class OpenADAS(AtomicData): @@ -32,8 +36,7 @@ def __init__(self, data_path=None, config=config.default, permit_extrapolation=F super().__init__() self._data_path = data_path or self._setup_data_path() - # configuration is immutable, changes to ADAS state are not tracked - self._config = config.copy() + self._config = config # if true informs interpolation objects to allow extrapolation beyond the limits of the tabulated data self._permit_extrapolation = permit_extrapolation @@ -208,3 +211,69 @@ def recombination_rate(self, ion, ionisation, transition): data = adf15(os.path.join(self._data_path, filename), block_number) return RecombinationRate(wavelength, data, extrapolate=self._permit_extrapolation) + def stage_resolved_line_radiation_rate(self, ion, ionisation): + + # extract element from isotope + if isinstance(ion, Isotope): + ion = ion.element + + try: + plt_files = ADF11_PLT_FILES[ion.symbol] + except KeyError: + raise ValueError("No ADF11 files set for Ion - {}".format(ion.symbol)) + + absolute_file_path = self._check_for_adf_file(plt_files['ADAS_Path'], plt_files['Download_URL']) + + densities, temperatures, rate_data = adf11(absolute_file_path, ion, ionisation) + + name = 'Stage Resolved Line Radiation - ({}, {})'.format(ion.symbol, ionisation) + return StageResolvedRadiation(ion, ionisation, densities, temperatures, rate_data, + name=name, extrapolate=self._permit_extrapolation) + + def stage_resolved_continuum_radiation_rate(self, ion, ionisation): + + # extract element from isotope + if isinstance(ion, Isotope): + ion = ion.element + + try: + prb_files = ADF11_PRB_FILES[ion.symbol] + except KeyError: + raise ValueError("No ADF11 files set for Ion - {}".format(ion.symbol)) + + absolute_file_path = self._check_for_adf_file(prb_files['ADAS_Path'], prb_files['Download_URL']) + + densities, temperatures, rate_data = adf11(absolute_file_path, ion, ionisation) + + name = 'Stage Resolved Continuum Radiation - ({}, {})'.format(ion.symbol, ionisation) + + return StageResolvedRadiation(ion, ionisation, densities, temperatures, rate_data, + name=name, extrapolate=self._permit_extrapolation) + + def _check_for_adf_file(self, relative_adf_file_path, download_path): + + relative_adf_directory, adf_file_name = os.path.split(relative_adf_file_path) + absolute_adf_directory = os.path.join(self._data_path, relative_adf_directory) + absolute_file_path = os.path.join(absolute_adf_directory, adf_file_name) + + if not os.path.exists(absolute_adf_directory): + os.makedirs(absolute_adf_directory) + + if os.path.isfile(absolute_file_path): + return absolute_file_path + else: + urllib.request.urlretrieve(download_path, absolute_file_path) + + return absolute_file_path + + def add_adf15_file(self, element, ionisation, adf_file_path): + + if not os.path.isfile(adf_file_path): + new_path = os.path.join(os.path.expanduser('~/.cherab/openadas'), adf_file_path) + if not os.path.isfile(new_path): + raise ValueError("Could not find ADF15 file - '{}'".format(adf_file_path)) + adf_file_path = new_path + + atomic_data_dict = RecursiveDict.from_dict(self._config) + atomic_data_dict = add_adf15_to_atomic_data(atomic_data_dict, element, ionisation, adf_file_path) + self._config = atomic_data_dict.freeze() diff --git a/cherab/openadas/rates/pec.pyx b/cherab/openadas/rates/pec.pyx index ef0dba2..5b9615b 100644 --- a/cherab/openadas/rates/pec.pyx +++ b/cherab/openadas/rates/pec.pyx @@ -42,7 +42,14 @@ cdef class ImpactExcitationRate(CoreImpactExcitationRate): ) cpdef double evaluate(self, double density, double temperature) except? -1e999: - return self._pec.evaluate(density, temperature) + + cdef double rate + + rate = self._pec.evaluate(density, temperature) + if rate < 0: + return 0.0 + + return rate # todo: evaluate it the interpolation can be done without the log10 operations? or if this should be ported to the cx rates. @@ -66,7 +73,14 @@ cdef class RecombinationRate(CoreRecombinationRate): cpdef double evaluate(self, double density, double temperature) except? -1e999: - return self._pec.evaluate(density, temperature) + + cdef double rate + + rate = self._pec.evaluate(density, temperature) + if rate < 0: + return 0.0 + + return rate # cdef class ThermalCXRate(CoreThermalCXRate): diff --git a/cherab/openadas/rates/radiated_power.pyx b/cherab/openadas/rates/radiated_power.pyx new file mode 100644 index 0000000..8d3c71f --- /dev/null +++ b/cherab/openadas/rates/radiated_power.pyx @@ -0,0 +1,58 @@ + +# Copyright 2014-2017 United Kingdom Atomic Energy Authority +# +# 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 numpy as np +from numpy cimport ndarray +from cherab.core.math.interpolators.interpolators2d cimport Interpolate2DCubic +from cherab.core.atomic.rates cimport StageResolvedLineRadiation as CoreStageResolvedLineRadiation + + +cdef class StageResolvedRadiation(CoreStageResolvedLineRadiation): + + cdef: + readonly bint extrapolate + readonly tuple density_range, temperature_range + readonly ndarray _electron_density, _electron_temperature, _radiated_power + readonly Interpolate2DCubic _power_func + + def __init__(self, element, ionisation, electron_density, electron_temperature, radiated_power, name='', extrapolate=False): + + super().__init__(element, ionisation, name=name) + + self._electron_density = np.array(electron_density, dtype=np.float64) + self._electron_temperature = np.array(electron_temperature, dtype=np.float64) + self._radiated_power = np.array(radiated_power, dtype=np.float64) + + if not all(j >= 0 for i in self._radiated_power for j in i): + raise ValueError("Can't have negative power values!") + + self.density_range = (self._electron_density.min(), self._electron_density.max()) + self.temperature_range = (self._electron_temperature.min(), self._electron_temperature.max()) + + self.extrapolate = extrapolate + self._power_func = Interpolate2DCubic(self._electron_density, self._electron_temperature, self._radiated_power, + extrapolate=extrapolate, extrapolation_type="nearest") + + cdef double evaluate(self, double electron_density, double electron_temperature) except? -1e999: + cdef double rate + rate = self._power_func.evaluate(electron_density, electron_temperature) + if rate < 0: + print("Warning - negative power encountered!") + print("ne - {:.4G}, te - {:.3G}, rate - {:.4G}".format(electron_density, electron_temperature, rate)) + return rate + # return self._power_func.evaluate(electron_density, electron_temperature) + diff --git a/cherab/openadas/read/__init__.py b/cherab/openadas/read/__init__.py index 6d18291..c267305 100644 --- a/cherab/openadas/read/__init__.py +++ b/cherab/openadas/read/__init__.py @@ -15,6 +15,7 @@ # under the Licence. +from .adf11 import adf11 from .adf12 import adf12 from .adf15 import adf15 from .adf21 import adf21 diff --git a/cherab/openadas/read/adf11.py b/cherab/openadas/read/adf11.py new file mode 100644 index 0000000..310db88 --- /dev/null +++ b/cherab/openadas/read/adf11.py @@ -0,0 +1,76 @@ +# Copyright 2014-2017 United Kingdom Atomic Energy Authority +# +# 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 numpy as np +import re + + +def adf11(file_path, ion, ionisation): + + adf11_fh = open(file_path, 'r') + lines = adf11_fh.readlines() + adf11_fh.close() + + file_header = lines.pop(0).split() + if int(file_header[0]) != ion.atomic_number: + raise ValueError("ADF file does not match required atomic number.") + + num_densities = int(file_header[1]) + num_temperatures = int(file_header[2]) + + lines.pop(0) + parameter_lines = [] + while True: + next_line = lines.pop(0) + if next_line[0:2] == '--': + break + components = next_line.split() + for c in components: + parameter_lines.append(float(c)) + + if len(parameter_lines) != num_densities + num_temperatures: + raise ValueError("ADF11 file could not be parsed correctly.") + + densities = np.array([10**x * 1E4 for x in parameter_lines[0:num_densities]]) + temperatures = np.array([10**x for x in parameter_lines[num_densities:]]) + + found = False + while True: + if next_line[0:2] == '--': + match = re.match('.*Z1= ([0-9]*).*', next_line) + if match and int(match.group(1)) == ionisation+1: + found = True + break + if next_line[0] == 'C': + break + next_line = lines.pop(0) + + if found: + rate_data = [] + while True: + next_line = lines.pop(0) + if next_line[0:2] == '--' or next_line[0] == 'C': + break + for c in next_line.split(): + rate_data.append(10**float(c) * 1E-8) + + rate_data = np.array(rate_data) + rate_data = rate_data.reshape((num_densities, num_temperatures)) + + else: + raise ValueError('Requested ADF11 data could not be found in this file.') + + return densities, temperatures, rate_data diff --git a/cherab/openadas/read/adf15.py b/cherab/openadas/read/adf15.py index db3425f..9d24ceb 100644 --- a/cherab/openadas/read/adf15.py +++ b/cherab/openadas/read/adf15.py @@ -137,6 +137,8 @@ def add_adf15_to_atomic_data(atomic_data_dictionary, element, ionisation, adf_fi atomic_data_dictionary[rate_type][element][ionisation][(upper_level, lower_level)] = (adf_file_path, block_num) atomic_data_dictionary["wavelength"][element][ionisation][(upper_level, lower_level)] = wavelength + return atomic_data_dictionary + # Group lines of file into blocks based on precursor ' 6561.9A 24...' def _group_by_block(source_file, match_string): diff --git a/demos/adf15_plots.py b/demos/adf15_plots.py new file mode 100644 index 0000000..51fcc23 --- /dev/null +++ b/demos/adf15_plots.py @@ -0,0 +1,10 @@ + +from cherab.core.atomic import carbon +from cherab.openadas import OpenADAS + + +adas = OpenADAS() + +myrate = adas.impact_excitation_rate(carbon, 1, ("2s1 2p1 3d1 2D4.5", "2s2 4d1 2D4.5")) + +