Source code for exo_k.chemistry

# -*- coding: utf-8 -*-
"""
@author: jeremy leconte
"""
import numpy as np
from scipy.interpolate import RectBivariateSpline
import pandas as pd
from .util.interp import bilinear_interpolation_array, interp_ind_weights

[docs] class EquChemTable(object): """Class to load and interpolate chemistry data. """ def __init__(self, filename=None, remove_zeros=False): """Initializes chemical composition tables from various init files Parameters ---------- filename : str Name of the input file """ self.filename=filename self.tgrid=None self.pgrid=None self.logpgrid=None self.p_unit=None self.Nt=None self.Np=None self.molecules=None self.vol_mix_ratio={} self.logx_interp={} if self.filename.lower().endswith('.in'): self.read_composition_in(filename=self.filename, remove_zeros=remove_zeros) elif self.filename.lower().endswith('.dat'): self.read_composition_dat(filename=self.filename, remove_zeros=remove_zeros) self.setup_interpolation()
[docs] def read_composition_in(self, filename=None, remove_zeros=False, skiprows=7): """Initializes chemical composition tables from composition.in files Parameters ---------- filename : str Name of the input file skiprows : int, optional Number of lines to skip in the file """ data = np.loadtxt(filename,skiprows=skiprows,unpack=True) self.Nt=np.where(data[2,1:]==data[2,0])[0][0]+1 self.tgrid=data[2,0:self.Nt] #print('T grid:',self.tgrid) self.Np=data[0].size//self.Nt self.pgrid=data[1,::self.Nt]/0.986923267*1.e5 self.p_unit='Pa' #print('P grid:',self.pgrid) #print('Np, Nt:',self.Np, self.Nt) self.logpgrid=np.log10(self.pgrid) with open(filename, 'r') as file: last_line=file.readline() while True: new_line=file.readline() if new_line.split()[0]=='z': break last_line=new_line molecules=last_line.split() #print(molecules) for ii, mol in enumerate(molecules): self.vol_mix_ratio[mol]=data[ii+3].reshape((self.Np,self.Nt)) if remove_zeros: self.remove_zeros()
[docs] def read_composition_dat(self, filename=None, remove_zeros=False, skiprows=1): """Initializes chemical composition tables from composition.in files Parameters ---------- filename : str Name of the input file skiprows : int, optional Number of lines to skip in the file """ with open(filename, 'r') as file: molecules=file.readline().split()[2:] molecules[0]=molecules[0].replace('[mol/mol]','') data = np.loadtxt(filename,skiprows=skiprows,unpack=True) self.pgrid=np.sort(np.array(list(set(data[0])))) self.logpgrid=np.log10(self.pgrid) self.tgrid=np.sort(np.array(list(set(data[1])))) self.Np=self.pgrid.size self.Nt=self.tgrid.size for ii, mol in enumerate(molecules): self.vol_mix_ratio[mol]=data[ii+2].reshape((self.Np,self.Nt))[:,-1::-1] if remove_zeros: self.remove_zeros()
[docs] def remove_zeros(self,deltalog_min_value=30.): """Finds zeros in the chem data and set them to (10.^-deltalog_min_value) times the minimum positive value in the table. This is to be able to work in logspace. """ for mol,x_ratio in self.vol_mix_ratio.items(): mask = np.zeros((self.Np,self.Nt),dtype=bool) mask[np.nonzero(x_ratio)] = True minvalue=np.amin(x_ratio[mask]) self.vol_mix_ratio[mol][~mask]=minvalue/(10.**deltalog_min_value)
[docs] def setup_interpolation(self): """Creates interpolating functions to be called later on. """ for mol,x_ratio in self.vol_mix_ratio.items(): self.logx_interp[mol]=RectBivariateSpline(self.logpgrid,self.tgrid,np.log(x_ratio))
def __getitem__(self, molecule): """Overrides getitem so that EquChemTable['molecule'] directly accesses the database for that molecule. """ if molecule not in self.vol_mix_ratio: raise KeyError('The requested molecule is not available.') return self.vol_mix_ratio[molecule]
[docs] def vmr(self, logP, T, mol): """Interpolates a single molecule on a logP-T value """ return np.exp(self.logx_interp[mol](logP,T,grid=False)[0])
[docs] def interpolate_vmr(self, logp_array=None,t_array=None, mols=None,grid=False): """Interpolates all molecules in mols on a logP-T grid """ res={} for mol in mols: res[mol]=np.exp(self.logx_interp[mol](logp_array, t_array, grid=grid)) return res
# def vmr_logPT(self, logPgrid, tgrid, mols, is_sorted=False): # """Interpolates all molecules in mols on a logP-T grid # """ # res={} # #lP,T=np.meshgrid(logPgrid,tgrid) # for mol in mols: # #res[mol]=10**self.logx_interp[mol](logPgrid.reshape((self.Np,1)),tgrid) # res[mol]=np.exp(self.logx_interp[mol](logPgrid, tgrid, \ # assume_sorted=is_sorted).transpose()) # return res
[docs] class InterpolationChemistry(object): """Chemistry model interpolating Volume Mixing Ratios based on a datafile provided by V. Parmentier. """ def __init__(self, filename=None, t_min=260): """ Parameters ---------- filename: str The path to the datafile. t_min: float minimal temperature to consider. Cannot be below 260K because not enough data in the file. """ self.df = pd.read_csv(filename,sep=r"\s{1,}", engine='python') self.df = self.df[self.df["T"] >= t_min] self.tgrid = np.unique(self.df["T"]) self.logpgrid = np.log10(np.exp(np.unique(self.df["P"])))+2. #self.logpgrid = np.unique(self.df["P"]) self.data = {} for gas in self.df.columns: if gas in ['T','P']: continue self.data[gas] = np.empty(self.logpgrid.shape + self.tgrid.shape) for i_t in range(self.tgrid.size): for i_p in range(self.logpgrid.size): if gas=='e-': pass self.data[gas][i_p, i_t] = self.df[gas].iloc[i_t * len(self.logpgrid) + i_p]
[docs] def compute_vmr(self, molecule = None, logp_array = None, t_array = None): """ Parameters ---------- molecule: str Name of molecule to deal with logp_array: array, np.ndarray Array of log10 pressures (Pa) t_array: array, np.ndarray Array of temperatures (K) """ if (molecule is None) or (logp_array is None) or (t_array is None): raise RuntimeError('Missing keyword argument.') logp_array=np.array(logp_array, dtype=float) t_array=np.array(t_array, dtype=float) i_p,weight_p = interp_ind_weights(logp_array, self.logpgrid) i_t,weight_t = interp_ind_weights(t_array, self.tgrid) p1t1=self.data[molecule][i_p, i_t ] p0t1=self.data[molecule][i_p-1,i_t ] p1t0=self.data[molecule][i_p, i_t-1] p0t0=self.data[molecule][i_p-1,i_t-1] return np.exp(bilinear_interpolation_array(p0t0, p1t0, p0t1, p1t1, weight_p, weight_t))
[docs] def compute_vmr_grid(self, molecule = None, logp_array = None, t_array = None): if (molecule is None) or (logp_array is None) or (t_array is None): raise RuntimeError('Missing keyword argument.') i_p,weight_p = interp_ind_weights(logp_array, self.logpgrid) i_t,weight_t = interp_ind_weights(t_array, self.tgrid) # now put everything on 2D grids. tmp=np.meshgrid(i_p, i_t, indexing='ij') shape=tmp[0].shape i_p_grid=tmp[0].ravel() i_t_grid=tmp[1].ravel() tmp=np.meshgrid(weight_p, weight_t, indexing='ij') weight_p_grid=tmp[0].ravel() weight_t_grid=tmp[1].ravel() p1t1=self.data[molecule][i_p_grid, i_t_grid ] p0t1=self.data[molecule][i_p_grid-1,i_t_grid ] p1t0=self.data[molecule][i_p_grid, i_t_grid-1] p0t0=self.data[molecule][i_p_grid-1,i_t_grid-1] return np.reshape(np.exp(bilinear_interpolation_array(p0t0, p1t0, p0t1, p1t1, weight_p_grid, weight_t_grid)), shape)