From ab7c15c5c30aa14513d88b79d6204906d4b94d0d Mon Sep 17 00:00:00 2001 From: Mario Romero Date: Wed, 12 Jun 2019 15:12:47 +0200 Subject: [PATCH 01/20] Mis clases para el gas --- phases/Example.py | 57 +++++++++++++++ phases/GasGenerics.py | 166 ++++++++++++++++++++++++++++++++++++++++++ phases/GasSpecies.py | 115 +++++++++++++++++++++++++++++ 3 files changed, 338 insertions(+) create mode 100644 phases/Example.py create mode 100644 phases/GasGenerics.py create mode 100644 phases/GasSpecies.py diff --git a/phases/Example.py b/phases/Example.py new file mode 100644 index 0000000..9299714 --- /dev/null +++ b/phases/Example.py @@ -0,0 +1,57 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +''' +An example + +Mario Romero June 2019 +''' + +import numpy as np +import astropy.units as u +import astropy.constants as cte +import GasGenerics as gas +import GasSpecies as species + +''' +Each species has its own class defined in 'GasSpecies.py' +''' + +#(1) Let's declare auxiliary variables + +uV = u.cm*u.cm*u.cm + +M0_a = 0.*u.solMass +T0_a = 100.*u.K +P0 = ((1.e4*u.K/uV)*cte.k_B).to(u.erg/uV) + +M0_i = 1000.*u.solMass +T0_i = 1000.*u.K + +dt = 1.*u.yr + +#(2) Let's declare hydrogen objects +HI = species.Neutral_Hydrogen() +HI.init(P=P0,T=T0_a,M=M0_a) #Your initial conditions goes here + +HII = species.Ionized_Hydrogen() +HII.init(P=P0,T=T0_i,M=M0_i) + +#(3) Now evolve + +t = 0.*u.yr +tf = 1000.*u.yr + + +while(t < tf): + print(t,HI.number_density(),HI.mass()) + print(t,HII.number_density(),HII.mass()) + print("\n") + + HI.update_derivatives([HI,HII]) + HII.update_derivatives([HI,HII]) + + HI.evolve(dt) + HII.evolve(dt) + t += dt + + #break diff --git a/phases/GasGenerics.py b/phases/GasGenerics.py new file mode 100644 index 0000000..c9523b2 --- /dev/null +++ b/phases/GasGenerics.py @@ -0,0 +1,166 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +''' +Ideal gas classes + +Mario Romero May 2019 +''' + +import numpy as np +import astropy.units as u +import astropy.constants as cte + + +''' +PARENT CLASS +''' +class Gas(MultiphaseMedium): +#class Gas(object): + + #--------------------- + #DEFAULT SETTINGS + #--------------------- + + def default_settings(self): + self._set() + return{ + 'current_mass':0.0*u.solMass, + 'dm_dt': 0.0*(u.solMass/u.yr), + + 'particle_mass':0.0*u.g, + } + + #--------------------- + #COMMON ATRIBUTES/METHODS TO ALL DERIVED CLASSES + #--------------------- + def _set(self): + + uV = (u.cm*u.cm*u.cm) #volume units + no_units = u.m/u.m #dimensionless units + #Internal variables + self._pressure = 0.0*(u.erg/u.s) + self._temperature = 0.0*u.K + self._number_particles = 0.0*no_units + self._thermal_energy = 0.0*u.erg + self._fugacity = 0.0*no_units + self._volume = 0.0*uV + self._gas_mass = 0.0*u.solMass + self._chemical_potential = 0.0*u.erg + self._number_density = 0.0 / uV + self._mass_density = 0.0*u.g/ uV + self._thermal_energy_density = 0.0*u.erg/uV + self._gamma = 0.0*no_units + + def _compute_variables_indepedent_of_state(self): + uV = (u.cm*u.cm*u.cm) #volume units + no_units = u.m/u.m #dimensionless units + + self._gas_mass = (self._number_particles * self._particle_mass).to(u.solMass) + self._number_density = (self._number_particles / self._volume).to(1./uV) + self._mass_density = (self._number_density * self._particle_mass).to(u.g/uV) + self._thermal_energy_density = (self._thermal_energy / self._volume).to(u.erg/uV) + self._chemical_potential = (cte.k_B*self._temperature*np.log(self._fugacity)).to(u.erg) + self._gamma = (1.+(self._pressure/self._thermal_energy_density)).to(no_units) + + #--------------------- + #CONSTRUCTORS + #--------------------- + #Update + def init(self,P=-1,T=-1,M=-1,N=-1): + self._state(P,T,M,N) + + #--------------------- + #OUTPUTS + #--------------------- + def mass(self,units=u.solMass): #Polymorph from Multiphase + return (self._gas_mass).to(units) + def pressure(self,units=u.erg/(u.cm*u.cm*u.cm)): + return (self._pressure).to(units) + def temperature(self,units=u.K): + return (self._temperature).to(units) + def number_particles(self,units=u.m/u.m): + return (self._number_particles).to(units) + def thermal_energy(self,units=u.erg): + return (self._thermal_energy).to(units) + def fugacity(self,units=u.m/u.m): + return (self._fugacity).to(units) + def volume(self,units=(u.cm*u.cm*u.cm)): + return (self._volume).to(uV) + def chemical_potential(self,units=u.erg): + return (self._chemical_potential).to(units) + def number_density(self,units=1./(u.cm*u.cm*u.cm)): + return (self._number_density).to(units) + def mass_density(self,units=u.g/(u.cm*u.cm*u.cm)): + return (self._mass_density).to(units) + def thermal_energy_density(self,units=u.erg/(u.cm*u.cm*u.cm)): + return (self._thermal_energy_density(self)).to(units) + def adiabatic_constant(self,units=u.m/u.m): + return (self._gamma).to(units) + + #--------------------- + #POLYMORPHIC METHODS + #--------------------- + #I.e.: These are the ONLY functions you have to redefine in each subclass + def _state(self,P,T,M,N):#inputs: pressure, temperature, number of particles, thermal energy, mass + raise NameError("Called by parent class. Should not happen!") + + def update_derivatives(self,term): + raise NameError("Called by parent class. Should not happen!") + +''' +CHILD CLASS +''' +#These classes must have '_state()' well defined and not polymorphic +#You can call this class as a generic monoatomic gas instead of a particular species, but you have to define the mass of its particles. +class Monoatomic(Gas): + + def default_settings(self): + self._set() + return{ + 'current_mass':0.0*u.solMass, + 'dm_dt': 0.0*(u.solMass/u.yr), + + 'particle_mass':0.0*u.g, + } + + def _state(self,P,T,M,N):#inputs: pressure, temperature, number of particles, thermal energy, mass + uV = (u.cm*u.cm*u.cm) #volume units + no_units = u.m/u.m #dimensionless units + + #Convert mass, if given + if M != -1: + N = (M.to(u.g))/(self._particle_mass).to(u.g) + + #You MUST give pressure and number of particles (or total mass) + if (P == -1) or (N == -1): + raise NameError("Gas._state() bad initialized! You must give pressure and total mass/number or particles!") + else: + #Initialize thermodynamic variables + self._pressure = P.to(u.erg/uV) + self._number_particles = N.to(no_units) + + #If temperature is given. I will follow a different path + if T != -1: + self._temperature = T.to(u.K) + + #We compute energy first + self._thermal_energy = (1.5*self._number_particles * (cte.k_B) * self._temperature).to(u.erg) + + #Volume... + self._volume = ((2.*self._thermal_energy) / (3.*self._pressure)).to(uV) + + #And fugacity + self._fugacity = (self._pressure * (2.*np.pi*self._particle_mass / (cte.h*cte.h) )**(-1.5) * (cte.k_B*self._temperature)**(-2.5)).to(no_units) + + else: + raise NameError("Gas.state bad initialized! You must give either temperature or internal energy") + + #Compute the other variables + self._compute_variables_indepedent_of_state() + + #--------------------- + #(STILL) POLYMORPHIC METHODS + #--------------------- + #Note that, this time, the method is well defined here as well + def update_derivatives(self,term): + raise NameError("Work in Progress...") diff --git a/phases/GasSpecies.py b/phases/GasSpecies.py new file mode 100644 index 0000000..7e64086 --- /dev/null +++ b/phases/GasSpecies.py @@ -0,0 +1,115 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +''' +Ideal gas classes + +Mario Romero May 2019 +''' + +import numpy as np +import astropy.units as u +import astropy.constants as cte +import GasGenerics as gas + + +''' +This script contains all realistic gas species (e.g. Hydrogen, Helium, etc). They are subclasses of classes from 'GasGenerics.py' + +You only need to redefine here 'update_derivatives' and 'default_settings' +''' + +class Neutral_Hydrogen(gas.Monoatomic): + + _particle_mass = (1.008*cte.u).to(u.g) #1.008*atomic_mass + + def __init__(self): + print("Warning: code written for testing purposes!") + self.default_settings() + + def default_settings(self): + self._set() + return{ + 'current_mass':0.0*u.solMass, + 'dm_dt': 0.0*(u.solMass/u.yr), + } + + def update_derivatives(self,term): + #DISCLAIMER: THIS IS AN EXAMPLE! + #print("Warning: code written for testing purposes!") + uV = u.cm*u.cm*u.cm + + self.dm_dt = 0.0*u.solMass/u.yr + dt = 0.0*u.yr + for element in term: + + add = 0.0*u.solMass/u.yr + + #Recombination (There is ionized hydrogen) + if (isinstance(element,Ionized_Hydrogen)): + + #ov = 4.1e-10*uV/u.s * (element._temperature/u.K)**(-0.8) + #tau = 1./(ov * (element._number_density).to(1./uV)) + tau = 50.*u.yr + + add = (element._gas_mass).to(u.solMass) / (tau).to(u.yr) + + self.dm_dt += add + + def evolve(self,dt): + + #Now evolve (simple Euler scheme) + self._gas_mass += self.dm_dt * (dt).to(u.yr) + + #And get the new state + self._state(self._pressure,self._temperature,self._gas_mass,-1) + + +class Ionized_Hydrogen(gas.Monoatomic): + + _particle_mass = (cte.m_p).to(u.g) #Proton mass + + def __init__(self): + print("Warning: code written for testing purposes!") + self.default_settings() + + def default_settings(self): + self._set() + return{ + 'current_mass':0.0*u.solMass, + 'dm_dt': 0.0*(u.solMass/u.yr), + } + + def update_derivatives(self,term): + #DISCLAIMER: THIS IS AN EXAMPLE! + #print("Warning: code written for testing purposes!") + + uV = u.cm*u.cm*u.cm + + self.dm_dt = 0.0*u.solMass/u.yr + dt = 0.0*u.yr + + + for element in term: + + add = 0.0*u.solMass/u.yr + + #Recombination (There is neutral hydrogen) + if (isinstance(element,Neutral_Hydrogen)): + + #ov = 4.1e-10*uV/u.s * (self._temperature/u.K)**(-0.8) + #tau = 1./(ov * (self._number_density).to(1./uV)) + tau = 50.*u.yr + + add = - (self._gas_mass).to(u.solMass) / ((tau).to(u.yr)) + + self.dm_dt += add + + #print(self.dm_dt) + + def evolve(self,dt): + + #Now evolve (simple Euler scheme) + self._gas_mass += self.dm_dt * (dt).to(u.yr) + + #And get the new state + self._state(self._pressure,self._temperature,self._gas_mass,-1) From 554f570c836e5c7adf507636f32a473f4d97be5b Mon Sep 17 00:00:00 2001 From: Mario Romero Date: Fri, 14 Jun 2019 15:55:58 +0200 Subject: [PATCH 02/20] =?UTF-8?q?[WIP]=20Clase=20para=20gases=20diatomicos?= =?UTF-8?q?=20a=C3=B1adida?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- phases/Example.py | 57 ------ phases/GasGenerics.py | 386 +++++++++++++++++++++++++---------------- phases/GasSpecies.py | 128 ++++---------- phases/HowToUse_gas.py | 63 +++++++ 4 files changed, 338 insertions(+), 296 deletions(-) delete mode 100644 phases/Example.py create mode 100644 phases/HowToUse_gas.py diff --git a/phases/Example.py b/phases/Example.py deleted file mode 100644 index 9299714..0000000 --- a/phases/Example.py +++ /dev/null @@ -1,57 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -''' -An example - -Mario Romero June 2019 -''' - -import numpy as np -import astropy.units as u -import astropy.constants as cte -import GasGenerics as gas -import GasSpecies as species - -''' -Each species has its own class defined in 'GasSpecies.py' -''' - -#(1) Let's declare auxiliary variables - -uV = u.cm*u.cm*u.cm - -M0_a = 0.*u.solMass -T0_a = 100.*u.K -P0 = ((1.e4*u.K/uV)*cte.k_B).to(u.erg/uV) - -M0_i = 1000.*u.solMass -T0_i = 1000.*u.K - -dt = 1.*u.yr - -#(2) Let's declare hydrogen objects -HI = species.Neutral_Hydrogen() -HI.init(P=P0,T=T0_a,M=M0_a) #Your initial conditions goes here - -HII = species.Ionized_Hydrogen() -HII.init(P=P0,T=T0_i,M=M0_i) - -#(3) Now evolve - -t = 0.*u.yr -tf = 1000.*u.yr - - -while(t < tf): - print(t,HI.number_density(),HI.mass()) - print(t,HII.number_density(),HII.mass()) - print("\n") - - HI.update_derivatives([HI,HII]) - HII.update_derivatives([HI,HII]) - - HI.evolve(dt) - HII.evolve(dt) - t += dt - - #break diff --git a/phases/GasGenerics.py b/phases/GasGenerics.py index c9523b2..54d7899 100644 --- a/phases/GasGenerics.py +++ b/phases/GasGenerics.py @@ -1,166 +1,254 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- ''' Ideal gas classes -Mario Romero May 2019 +Mario Romero May 2019 ''' import numpy as np import astropy.units as u import astropy.constants as cte - +import basic ''' -PARENT CLASS +PARENT CLASSES +classes defined here should not be defined as objects in the main code. ''' -class Gas(MultiphaseMedium): +class Gas(basic.MultiphaseMedium): #class Gas(object): - #--------------------- - #DEFAULT SETTINGS - #--------------------- - - def default_settings(self): - self._set() - return{ - 'current_mass':0.0*u.solMass, - 'dm_dt': 0.0*(u.solMass/u.yr), - - 'particle_mass':0.0*u.g, - } - - #--------------------- - #COMMON ATRIBUTES/METHODS TO ALL DERIVED CLASSES - #--------------------- - def _set(self): - - uV = (u.cm*u.cm*u.cm) #volume units - no_units = u.m/u.m #dimensionless units - #Internal variables - self._pressure = 0.0*(u.erg/u.s) - self._temperature = 0.0*u.K - self._number_particles = 0.0*no_units - self._thermal_energy = 0.0*u.erg - self._fugacity = 0.0*no_units - self._volume = 0.0*uV - self._gas_mass = 0.0*u.solMass - self._chemical_potential = 0.0*u.erg - self._number_density = 0.0 / uV - self._mass_density = 0.0*u.g/ uV - self._thermal_energy_density = 0.0*u.erg/uV - self._gamma = 0.0*no_units - - def _compute_variables_indepedent_of_state(self): - uV = (u.cm*u.cm*u.cm) #volume units - no_units = u.m/u.m #dimensionless units - - self._gas_mass = (self._number_particles * self._particle_mass).to(u.solMass) - self._number_density = (self._number_particles / self._volume).to(1./uV) - self._mass_density = (self._number_density * self._particle_mass).to(u.g/uV) - self._thermal_energy_density = (self._thermal_energy / self._volume).to(u.erg/uV) - self._chemical_potential = (cte.k_B*self._temperature*np.log(self._fugacity)).to(u.erg) - self._gamma = (1.+(self._pressure/self._thermal_energy_density)).to(no_units) - - #--------------------- - #CONSTRUCTORS - #--------------------- - #Update - def init(self,P=-1,T=-1,M=-1,N=-1): - self._state(P,T,M,N) - - #--------------------- - #OUTPUTS - #--------------------- - def mass(self,units=u.solMass): #Polymorph from Multiphase - return (self._gas_mass).to(units) - def pressure(self,units=u.erg/(u.cm*u.cm*u.cm)): - return (self._pressure).to(units) - def temperature(self,units=u.K): - return (self._temperature).to(units) - def number_particles(self,units=u.m/u.m): - return (self._number_particles).to(units) - def thermal_energy(self,units=u.erg): - return (self._thermal_energy).to(units) - def fugacity(self,units=u.m/u.m): - return (self._fugacity).to(units) - def volume(self,units=(u.cm*u.cm*u.cm)): - return (self._volume).to(uV) - def chemical_potential(self,units=u.erg): - return (self._chemical_potential).to(units) - def number_density(self,units=1./(u.cm*u.cm*u.cm)): - return (self._number_density).to(units) - def mass_density(self,units=u.g/(u.cm*u.cm*u.cm)): - return (self._mass_density).to(units) - def thermal_energy_density(self,units=u.erg/(u.cm*u.cm*u.cm)): - return (self._thermal_energy_density(self)).to(units) - def adiabatic_constant(self,units=u.m/u.m): - return (self._gamma).to(units) - - #--------------------- - #POLYMORPHIC METHODS - #--------------------- - #I.e.: These are the ONLY functions you have to redefine in each subclass - def _state(self,P,T,M,N):#inputs: pressure, temperature, number of particles, thermal energy, mass - raise NameError("Called by parent class. Should not happen!") - - def update_derivatives(self,term): - raise NameError("Called by parent class. Should not happen!") + #--------------------- + #DEFAULT SETTINGS + #--------------------- + def default_settings(self): + self._set() + return{ + 'current_mass':0.0*u.solMass, + 'dm_dt': 0.0*(u.solMass/u.yr), + } + + #--------------------- + #COMMON ATRIBUTES/METHODS TO ALL DERIVED CLASSES + #--------------------- + def _set(self): + + uV = (u.cm*u.cm*u.cm) #volume units + no_units = u.m/u.m #dimensionless units + #Internal variables + self._pressure = 0.0*(u.erg/u.s) + self._temperature = 0.0*u.K + self._number_particles = 0.0*no_units + self._thermal_energy = 0.0*u.erg + self._fugacity = 0.0*no_units + self._volume = 0.0*uV + self._gas_mass = 0.0*u.solMass + self._chemical_potential = 0.0*u.erg + self._number_density = 0.0 / uV + self._mass_density = 0.0*u.g/ uV + self._thermal_energy_density = 0.0*u.erg/uV + self._gamma = 0.0*no_units + + def _compute_variables_indepedent_of_state(self): + uV = (u.cm*u.cm*u.cm) #volume units + no_units = u.m/u.m #dimensionless units + + self._gas_mass = (self._number_particles * self._particle_mass).to(u.solMass) + self._number_density = (self._number_particles / self._volume).to(1./uV) + self._mass_density = (self._number_density * self._particle_mass).to(u.g/uV) + self._thermal_energy_density = (self._thermal_energy / self._volume).to(u.erg/uV) + self._chemical_potential = (cte.k_B*self._temperature*np.log(self._fugacity)).to(u.erg) + self._gamma = (1.+(self._pressure/self._thermal_energy_density)).to(no_units) + + #--------------------- + #CONSTRUCTORS + #--------------------- + def __init__(self): + print("Warning: '__init__()' code written for testing purposes!") + self.default_settings() + + #--------------------- + #OUTPUTS + #--------------------- + def mass(self,units=u.solMass): #Polymorph from Multiphase + return (self._gas_mass).to(units) + def pressure(self,units=u.erg/(u.cm*u.cm*u.cm)): + return (self._pressure).to(units) + def temperature(self,units=u.K): + return (self._temperature).to(units) + def number_particles(self,units=u.m/u.m): + return (self._number_particles).to(units) + def thermal_energy(self,units=u.erg): + return (self._thermal_energy).to(units) + def fugacity(self,units=u.m/u.m): + return (self._fugacity).to(units) + def volume(self,units=(u.cm*u.cm*u.cm)): + return (self._volume).to(units) + def chemical_potential(self,units=u.erg): + return (self._chemical_potential).to(units) + def number_density(self,units=1./(u.cm*u.cm*u.cm)): + return (self._number_density).to(units) + def mass_density(self,units=u.g/(u.cm*u.cm*u.cm)): + return (self._mass_density).to(units) + def thermal_energy_density(self,units=u.erg/(u.cm*u.cm*u.cm)): + return (self._thermal_energy_density(self)).to(units) + def adiabatic_constant(self,units=u.m/u.m): + return (self._gamma).to(units) + + #--------------------- + #POLYMORPHIC METHODS + #--------------------- + #I.e.: These are the ONLY functions you have to redefine in each subclass + def state(self,P=-1,T=-1,M=-1,N=-1):#inputs: pressure, temperature, number of particles, thermal energy, mass + raise NameError("Called by parent class. Should not happen!") + + def update_derivatives(self,term): + raise NameError("Called by parent class. Should not happen!") ''' -CHILD CLASS +CHILDS OF GAS +These classes must have 'state()' well defined and not polymorphic ''' -#These classes must have '_state()' well defined and not polymorphic -#You can call this class as a generic monoatomic gas instead of a particular species, but you have to define the mass of its particles. class Monoatomic(Gas): - - def default_settings(self): - self._set() - return{ - 'current_mass':0.0*u.solMass, - 'dm_dt': 0.0*(u.solMass/u.yr), - - 'particle_mass':0.0*u.g, - } - - def _state(self,P,T,M,N):#inputs: pressure, temperature, number of particles, thermal energy, mass - uV = (u.cm*u.cm*u.cm) #volume units - no_units = u.m/u.m #dimensionless units - - #Convert mass, if given - if M != -1: - N = (M.to(u.g))/(self._particle_mass).to(u.g) - - #You MUST give pressure and number of particles (or total mass) - if (P == -1) or (N == -1): - raise NameError("Gas._state() bad initialized! You must give pressure and total mass/number or particles!") - else: - #Initialize thermodynamic variables - self._pressure = P.to(u.erg/uV) - self._number_particles = N.to(no_units) - - #If temperature is given. I will follow a different path - if T != -1: - self._temperature = T.to(u.K) - - #We compute energy first - self._thermal_energy = (1.5*self._number_particles * (cte.k_B) * self._temperature).to(u.erg) - - #Volume... - self._volume = ((2.*self._thermal_energy) / (3.*self._pressure)).to(uV) - - #And fugacity - self._fugacity = (self._pressure * (2.*np.pi*self._particle_mass / (cte.h*cte.h) )**(-1.5) * (cte.k_B*self._temperature)**(-2.5)).to(no_units) - - else: - raise NameError("Gas.state bad initialized! You must give either temperature or internal energy") - - #Compute the other variables - self._compute_variables_indepedent_of_state() - - #--------------------- - #(STILL) POLYMORPHIC METHODS - #--------------------- - #Note that, this time, the method is well defined here as well - def update_derivatives(self,term): - raise NameError("Work in Progress...") + + #--------------------- + #STATE METHOD + #--------------------- + def state(self,P=-1,T=-1,M=-1,N=-1):#inputs: pressure, temperature, number of particles, thermal energy, mass + uV = (u.cm*u.cm*u.cm) #volume units + no_units = u.m/u.m #dimensionless units + + #Convert mass, if given + if M != -1: + N = (M.to(u.g))/(self._particle_mass).to(u.g) + + #You MUST give pressure, temperature and number of particles (or total mass) + if (P == -1) or (T == -1) or (N == -1): + raise NameError("Gas._state() bad initialized! You must give pressure, temperature and total mass/number or particles!") + else: + #Initialize thermodynamic variables + self._pressure = P.to(u.erg/uV) + self._number_particles = N.to(no_units) + + self._temperature = T.to(u.K) + + #We compute energy first + self._thermal_energy = (1.5*self._number_particles * (cte.k_B) * self._temperature).to(u.erg) + + #Volume... + self._volume = ((2.*self._thermal_energy) / (3.*self._pressure)).to(uV) + + #And fugacity + self._fugacity = (self._pressure * (2.*np.pi*self._particle_mass / (cte.h*cte.h) )**(-1.5) * (cte.k_B*self._temperature)**(-2.5)).to(no_units) + + #Compute the other variables + self._compute_variables_indepedent_of_state() + +class Diatomic(Gas): + + #--------------------- + #AUXILIARY METHODS + #--------------------- + def _molecule_mass(self): + return (self._atom_mass[0]+self._atom_mass[1]).to(u.g) + + def _moment_of_inertia(self): + #Get reduced mass + + m_r = ( self._atom_mass[0]*self._atom_mass[1] / (self._atom_mass[0]+self._atom_mass[1]) ).to(u.g) + return m_r*self._atom_distance*self._atom_distance + + def _angular_frequency(self): + return (self._wavenumber * cte.c ).to(1./u.s) + + #--------------------- + #PARTITION FUNCTION METHODS + #--------------------- + def _Z_rot(self): + + no_units = u.cm/u.cm + #Compute a parameter to check if quantum effects are important + theta = ( 0.5 * cte.hbar*cte.hbar / (cte.k_B * self._temperature * self._moment_of_inertia()) ).to(no_units) + + #Now do stuff + Z = None + dZ_db = None #Partial derivative of Z with respect to beta = 1./kT + if theta <= 1: #All rotational modes are enabled (=classical) + Z = 1./theta + dZ_db = - cte.k_B * self._temperature / theta + else: #All rotational modes are disabled + Z = 1.0 + dZ_db = 0.0 + + return [Z*no_units , dZ_db*(u.erg)] + ''' + Note that this is a simplification. + The exact method has Z = sum (2l+1)*exp(-theta*l(l+1)) ; for l=0 to infinity. + Be aware that it's computationally expensive. + ''' + + def _Z_vib(self): + + no_units = u.cm/u.cm + #Compute a parameter to check if quantum effects are important + xi = ( cte.hbar * self._angular_frequency() / (cte.k_B * self._temperature) ).to(no_units) + + #Now do stuff + Z = None + dZ_db = None #Partial derivative of Z with respect to beta = 1./kT + if xi > 1: #All vibrational modes are disabled + Z = 1.0 + dZ_db = 0.0 + else: + Z = 1./xi + dZ_db = - cte.k_B * self._temperature / xi + + return [Z*no_units , dZ_db*(u.erg)] + ''' + Again, this is a simplification. + The exact method has Z = sum exp(-n*xi) = 1./(1-exp(-xi)) . + Note that the zero of energies is set as all vibrations at the fundamental state of the harmonic oscillator. + ''' + + #--------------------- + #STATE METHOD + #--------------------- + def state(self,P=-1,T=-1,M=-1,N=-1):#inputs: pressure, temperature, number of particles, thermal energy, mass + + uV = (u.cm*u.cm*u.cm) #volume units + no_units = u.m/u.m #dimensionless units + + #Because you give each atom mass individually, you need to construct the particle mass first. + self._particle_mass = self._molecule_mass() + + #Convert mass, if given + if M != -1: + N = (M.to(u.g))/(self._particle_mass).to(u.g) + + #You MUST give pressure, temperature and number of particles (or total mass) + if (P == -1) or (T == -1) or (N == -1): + raise NameError("Gas._state() bad initialized! You must give pressure, temperature and total mass/number or particles!") + else: + #Initialize thermodynamic variables + self._pressure = P.to(u.erg/uV) + self._number_particles = N.to(no_units) + self._temperature = T.to(u.K) + + #Useful variable + NkT = (self._number_particles * (cte.k_B).to(u.erg/u.K) * self._temperature).to(u.erg) + + #Partition functions (they are arrays, [0] is the Z itself, and [1] the derivative with respect to beta) + Zrot = self._Z_rot() + Zvib = self._Z_vib() + + #We compute volume first (it's the ideal gas relation) + self._volume = ( NkT / self._pressure ).to(uV) + + #Fugacity... + self._fugacity = (self._pressure * (2.*np.pi*self._particle_mass / (cte.h*cte.h) )**(-1.5) * (cte.k_B*self._temperature)**(-2.5) / (Zrot[0]*Zvib[0]) ).to(no_units) + + #Thermal energy... + const_factor = 1.5 - (Zrot[1]/( Zrot[0] * cte.k_B * self._temperature )).to(no_units) - (Zvib[1]/( Zvib[0] * cte.k_B * self._temperature )).to(no_units) + self._thermal_energy = NkT * const_factor + + + #Compute the other variables + self._compute_variables_indepedent_of_state() \ No newline at end of file diff --git a/phases/GasSpecies.py b/phases/GasSpecies.py index 7e64086..3efc6d5 100644 --- a/phases/GasSpecies.py +++ b/phases/GasSpecies.py @@ -3,7 +3,7 @@ ''' Ideal gas classes -Mario Romero May 2019 +Mario Romero May 2019 ''' import numpy as np @@ -13,103 +13,51 @@ ''' +SPECIFIC GASES This script contains all realistic gas species (e.g. Hydrogen, Helium, etc). They are subclasses of classes from 'GasGenerics.py' - -You only need to redefine here 'update_derivatives' and 'default_settings' ''' +#-------------- +#MONOATOMIC SPECIES +#Required data = particle mass +#-------------- class Neutral_Hydrogen(gas.Monoatomic): - _particle_mass = (1.008*cte.u).to(u.g) #1.008*atomic_mass - - def __init__(self): - print("Warning: code written for testing purposes!") - self.default_settings() - - def default_settings(self): - self._set() - return{ - 'current_mass':0.0*u.solMass, - 'dm_dt': 0.0*(u.solMass/u.yr), - } - - def update_derivatives(self,term): - #DISCLAIMER: THIS IS AN EXAMPLE! - #print("Warning: code written for testing purposes!") - uV = u.cm*u.cm*u.cm - - self.dm_dt = 0.0*u.solMass/u.yr - dt = 0.0*u.yr - for element in term: - - add = 0.0*u.solMass/u.yr - - #Recombination (There is ionized hydrogen) - if (isinstance(element,Ionized_Hydrogen)): + _particle_mass = (1.008*cte.u).to(u.g) #1.008*atomic_mass - #ov = 4.1e-10*uV/u.s * (element._temperature/u.K)**(-0.8) - #tau = 1./(ov * (element._number_density).to(1./uV)) - tau = 50.*u.yr + def update_derivatives(self,term): + print("Neutral_Hydrogen: Work in Progress...") - add = (element._gas_mass).to(u.solMass) / (tau).to(u.yr) - - self.dm_dt += add - - def evolve(self,dt): - - #Now evolve (simple Euler scheme) - self._gas_mass += self.dm_dt * (dt).to(u.yr) - - #And get the new state - self._state(self._pressure,self._temperature,self._gas_mass,-1) + def evolve(self,dt): + print("Neutral_Hydrogen: 'evolve()' written for testing purposes!") + self.state(P=self._pressure,T=self._temperature,M=self._gas_mass) class Ionized_Hydrogen(gas.Monoatomic): - _particle_mass = (cte.m_p).to(u.g) #Proton mass - - def __init__(self): - print("Warning: code written for testing purposes!") - self.default_settings() - - def default_settings(self): - self._set() - return{ - 'current_mass':0.0*u.solMass, - 'dm_dt': 0.0*(u.solMass/u.yr), - } - - def update_derivatives(self,term): - #DISCLAIMER: THIS IS AN EXAMPLE! - #print("Warning: code written for testing purposes!") - - uV = u.cm*u.cm*u.cm - - self.dm_dt = 0.0*u.solMass/u.yr - dt = 0.0*u.yr - - - for element in term: - - add = 0.0*u.solMass/u.yr - - #Recombination (There is neutral hydrogen) - if (isinstance(element,Neutral_Hydrogen)): - - #ov = 4.1e-10*uV/u.s * (self._temperature/u.K)**(-0.8) - #tau = 1./(ov * (self._number_density).to(1./uV)) - tau = 50.*u.yr - - add = - (self._gas_mass).to(u.solMass) / ((tau).to(u.yr)) - - self.dm_dt += add - - #print(self.dm_dt) - - def evolve(self,dt): - - #Now evolve (simple Euler scheme) - self._gas_mass += self.dm_dt * (dt).to(u.yr) - - #And get the new state - self._state(self._pressure,self._temperature,self._gas_mass,-1) + _particle_mass = (cte.m_p).to(u.g) #Proton mass + + def update_derivatives(self,term): + print("Ionized_Hydrogen: Work in Progress...") + + def evolve(self,dt): + print("Ionized_Hydrogen: 'evolve()' written for testing purposes!") + self.state(P=self._pressure,T=self._temperature,M=self._gas_mass) + +#-------------- +#DIATOMIC SPECIES +#Needed data: particle masses (array), distance between particles and their fundamental vibration mode. +#-------------- +class Molecular_Hydrogen(gas.Diatomic): + + _atom_mass = [(1.008*cte.u).to(u.g) , (1.008*cte.u).to(u.g)] #Two hydrogen atoms + _atom_distance = 74.14*1e-10 * u.cm #Taken from wikipedia + _wavenumber = 4342.0*(1./u.cm) #wavenumber (k) of fundamental vibration. Taken from here: https://www.chem.purdue.edu/gchelp/vibs/h2.html // + + + def update_derivatives(self,term): + print("Molecular_Hydrogen: Work in Progress...") + + def evolve(self,dt): + print("Molecular_Hydrogen: 'evolve()' written for testing purposes!") + self.state(P=self._pressure,T=self._temperature,M=self._gas_mass) diff --git a/phases/HowToUse_gas.py b/phases/HowToUse_gas.py new file mode 100644 index 0000000..cb993fd --- /dev/null +++ b/phases/HowToUse_gas.py @@ -0,0 +1,63 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +''' +How to use the gas species + +Mario Romero June 2019 +''' + +import astropy.units as u +import astropy.constants as cte +import GasSpecies as species + +''' +Each species has its own class defined in 'GasSpecies.py' +''' + +#(1) Let's declare auxiliary variables + +uV = u.cm*u.cm*u.cm #units +P0 = ((1.e4*u.K/uV)*cte.k_B).to(u.erg/uV) #Common pressure to all phases + + +#(2) Let's declare the initial conditions of hydrogen +#molecular phase +M0_m = 10.*u.solMass +T0_m = 10.*u.K + +#atomic phase +M0_a = 100.*u.solMass +T0_a = 100.*u.K + +#ionized phase +M0_i = 1000.*u.solMass +T0_i = 1000.*u.K + +#(3) Let's declare hydrogen objects +H2 = species.Molecular_Hydrogen() +H2.state(P=P0,T=T0_m,M=M0_m) #Your initial conditions goes here + +HI = species.Neutral_Hydrogen() +HI.state(P=P0,T=T0_a,M=M0_a) + +HII = species.Ionized_Hydrogen() +HII.state(P=P0,T=T0_i,M=M0_i) + +#(4) Do stuff (I'm throwing arbitrary numbers) +dt = 1.0*u.yr +t = 0.0*u.yr +tf = 0.5*u.yr + +while(t < tf): + + #Get the derivatives + H2.update_derivatives([H2,HI,HII]) + HI.update_derivatives([H2,HI,HII]) + HII.update_derivatives([H2,HI,HII]) + + #Use your favourite integrator + H2.evolve(dt) + HI.evolve(dt) + HII.evolve(dt) + + t += dt From 561cf23ea7a58c8c9da5d2c4e5e3d78eb718e6d9 Mon Sep 17 00:00:00 2001 From: Mario Romero Date: Fri, 14 Jun 2019 16:41:57 +0200 Subject: [PATCH 03/20] [WIP] Quitando bugs --- phases/GasGenerics.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/phases/GasGenerics.py b/phases/GasGenerics.py index 54d7899..8fab4ee 100644 --- a/phases/GasGenerics.py +++ b/phases/GasGenerics.py @@ -175,10 +175,10 @@ def _Z_rot(self): Z = 1./theta dZ_db = - cte.k_B * self._temperature / theta else: #All rotational modes are disabled - Z = 1.0 - dZ_db = 0.0 + Z = 1.0*no_units + dZ_db = 0.0*u.erg - return [Z*no_units , dZ_db*(u.erg)] + return [Z , dZ_db] ''' Note that this is a simplification. The exact method has Z = sum (2l+1)*exp(-theta*l(l+1)) ; for l=0 to infinity. @@ -195,13 +195,13 @@ def _Z_vib(self): Z = None dZ_db = None #Partial derivative of Z with respect to beta = 1./kT if xi > 1: #All vibrational modes are disabled - Z = 1.0 - dZ_db = 0.0 + Z = 1.0*no_units + dZ_db = 0.0*u.erg else: Z = 1./xi dZ_db = - cte.k_B * self._temperature / xi - return [Z*no_units , dZ_db*(u.erg)] + return [Z , dZ_db] ''' Again, this is a simplification. The exact method has Z = sum exp(-n*xi) = 1./(1-exp(-xi)) . From 8fc1feb64fea4996c0dcf8237ba59c685a652513 Mon Sep 17 00:00:00 2001 From: Mario Romero Date: Wed, 26 Jun 2019 12:44:02 +0200 Subject: [PATCH 04/20] Cambios recientes --- phases/GasGenerics.py | 22 ++++++++++++++-------- phases/GasSpecies.py | 2 -- phases/HowToUse_gas.py | 1 + phases/__init__.py | 4 ++++ settings.example.yml | 17 ++++++++++++++++- settings.py | 1 + 6 files changed, 36 insertions(+), 11 deletions(-) diff --git a/phases/GasGenerics.py b/phases/GasGenerics.py index 8fab4ee..0cf027e 100644 --- a/phases/GasGenerics.py +++ b/phases/GasGenerics.py @@ -19,11 +19,17 @@ class Gas(basic.MultiphaseMedium): #--------------------- #DEFAULT SETTINGS #--------------------- + def __init__(self, params): + self.params = {**self.default_settings(), **params} + self._state() #(!) Mirar si hay que meter parametros o no + def default_settings(self): - self._set() return{ - 'current_mass':0.0*u.solMass, + 'mass':0.0*u.solMass, 'dm_dt': 0.0*(u.solMass/u.yr), + + 'temperature': 0.0*u.K, + 'pressure':0.0*(u.erg/(u.cm*u.cm*u.cm)) } #--------------------- @@ -34,13 +40,13 @@ def _set(self): uV = (u.cm*u.cm*u.cm) #volume units no_units = u.m/u.m #dimensionless units #Internal variables - self._pressure = 0.0*(u.erg/u.s) - self._temperature = 0.0*u.K + self._pressure = self.params["pressure"].to(u.erg/u.s) + self._temperature = self.params["temperature"].to(u.K) self._number_particles = 0.0*no_units self._thermal_energy = 0.0*u.erg self._fugacity = 0.0*no_units self._volume = 0.0*uV - self._gas_mass = 0.0*u.solMass + self._gas_mass = self.params["mass"].to(u.solMass) self._chemical_potential = 0.0*u.erg self._number_density = 0.0 / uV self._mass_density = 0.0*u.g/ uV @@ -61,9 +67,9 @@ def _compute_variables_indepedent_of_state(self): #--------------------- #CONSTRUCTORS #--------------------- - def __init__(self): - print("Warning: '__init__()' code written for testing purposes!") - self.default_settings() + #def __init__(self): + # print("Warning: '__init__()' code written for testing purposes!") + # self.default_settings() #--------------------- #OUTPUTS diff --git a/phases/GasSpecies.py b/phases/GasSpecies.py index 3efc6d5..e4a4ea3 100644 --- a/phases/GasSpecies.py +++ b/phases/GasSpecies.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- ''' Ideal gas classes diff --git a/phases/HowToUse_gas.py b/phases/HowToUse_gas.py index cb993fd..564cd66 100644 --- a/phases/HowToUse_gas.py +++ b/phases/HowToUse_gas.py @@ -34,6 +34,7 @@ T0_i = 1000.*u.K #(3) Let's declare hydrogen objects +#(!) Cambiar esto. H2 = species.Molecular_Hydrogen() H2.state(P=P0,T=T0_m,M=M0_m) #Your initial conditions goes here diff --git a/phases/__init__.py b/phases/__init__.py index 9cfa595..71f00df 100644 --- a/phases/__init__.py +++ b/phases/__init__.py @@ -1,10 +1,14 @@ from . import basic +from . import GasSpecies from .dust import Dust def select_phase(name): available_phases = { "simple": basic.Phase, "multiphase": basic.MultiphaseMedium, + "ionized_hydrogen": GasSpecies.Ionized_Hydrogen, + "neutral_hydrogen": GasSpecies.Neutral_Hydrogen, + "molecular_hydrogen": GasSpecies.Molecular_Hydrogen, "dust": Dust, } diff --git a/settings.example.yml b/settings.example.yml index d96dcaf..7e80414 100644 --- a/settings.example.yml +++ b/settings.example.yml @@ -10,9 +10,24 @@ dtd_sn: rlp phases: gas: - type: simple + type: simple #NOMBRE DE UNA CLASE params: mass: 2.34 + type: ionized_hydrogen + params: + mass: 2.34 + temperature: 1.0e6 + pressure: 1.0e4 + type: neutral_hydrogen + params: + mass: 2.34 + temperature: 1.0e4 + pressure: 1.0e4 + type: molecular_hydrogen + params: + mass: 2.34 + temperature: 1.0e2 + pressure: 1.0e4 dust: type: dust params: diff --git a/settings.py b/settings.py index 54735f5..fed87e3 100644 --- a/settings.py +++ b/settings.py @@ -22,6 +22,7 @@ valid_values = { 'imf': ['salpeter', 'starburst', 'chabrier', 'ferrini', 'kroupa', 'miller_scalo', 'maschberger'], 'dtd_sn': ['rlp', 'mdvp'], # rlp = Ruiz-Lapuente, mdvp = Mannucci, Della Valle, Panagia (2006) + 'phases': ['ionized_hydrogen', 'neutral_hydrogen', 'molecular_hydrogen'], } def validate(params): From 94c625270d2f8dfd980f0ff98f713dac60014dc1 Mon Sep 17 00:00:00 2001 From: Mario Romero Date: Fri, 12 Jul 2019 11:08:22 +0200 Subject: [PATCH 05/20] needed for rebase --- phases/GasGenerics.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/phases/GasGenerics.py b/phases/GasGenerics.py index 0cf027e..301f614 100644 --- a/phases/GasGenerics.py +++ b/phases/GasGenerics.py @@ -20,8 +20,8 @@ class Gas(basic.MultiphaseMedium): #DEFAULT SETTINGS #--------------------- def __init__(self, params): - self.params = {**self.default_settings(), **params} - self._state() #(!) Mirar si hay que meter parametros o no + self.params = {**self.default_settings(), **params} #Falta 'model' + self._state(P=self.params["pressure"],T=self.params["temperature"],M=self.params["mass"]) def default_settings(self): return{ @@ -217,7 +217,7 @@ def _Z_vib(self): #--------------------- #STATE METHOD #--------------------- - def state(self,P=-1,T=-1,M=-1,N=-1):#inputs: pressure, temperature, number of particles, thermal energy, mass + def state(self,P=-1,T=-1,M=-1,N=-1):#inputs: pressure, temperature, number of particles, mass uV = (u.cm*u.cm*u.cm) #volume units no_units = u.m/u.m #dimensionless units From 6c001863f88871cc98fa7ca161d55f45de97953c Mon Sep 17 00:00:00 2001 From: Yago Date: Thu, 23 May 2019 11:44:55 +0200 Subject: [PATCH 06/20] New SF branch (constant efficiency) --- .../star_formation/constant_efficiency.py | 24 ++++++++++++------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/processes/star_formation/constant_efficiency.py b/processes/star_formation/constant_efficiency.py index bdd2331..c235415 100644 --- a/processes/star_formation/constant_efficiency.py +++ b/processes/star_formation/constant_efficiency.py @@ -6,15 +6,23 @@ @author: yago """ -tau = 3. # gas depletion time, in Gyr +# ----------------------------------------------------------------------------- +class Process(object): -def init(params): - global tau - tau = float(params['tau_SF']) + def compute_derivatives(): + raise("ERROR: compute_derivatives method not implemented") -def compute_derivatives(): - SFR = gas.mass()/tau - gas.update_derivatives(-SFR) - stars.update_derivatives(SFR) +# ----------------------------------------------------------------------------- +class Constant_efficiency(Process): + + def __init__(self, gas, stars, tau_SF): + self.gas = gas # gas phase that is converted into stars + self.tau = tau_SF # gas depletion time, in Gyr + self.stars = stars # destination stellar phase + + def compute_derivatives(self): + SFR = self.gas.mass()/self.tau + self.gas.update_derivatives(-SFR) + self.stars.update_derivatives(SFR) From 973d230979d17edd197cd26d19a44ffe3d2a156d Mon Sep 17 00:00:00 2001 From: Yago Date: Wed, 26 Jun 2019 14:00:12 +0200 Subject: [PATCH 07/20] Updated structure of phases and processes --- .gitignore | 3 + SPiCE.py | 68 ++++++++++-------- model.example.yml | 26 +++++++ phases/__init__.py | 10 +-- phases/basic.py | 32 ++++++--- phases/example/__init__.py | 1 + phases/{ => example}/dust.py | 2 +- processes/__init__.py | 6 +- processes/basic.py | 30 ++++++++ processes/registry.py | 10 --- processes/star_formation.py | 23 ++++++ processes/star_formation/__init__.py | 1 - .../__pycache__/__init__.cpython-36.pyc | Bin 193 -> 0 bytes .../constant_efficiency.cpython-36.pyc | Bin 566 -> 0 bytes settings.example.yml | 41 ----------- 15 files changed, 155 insertions(+), 98 deletions(-) create mode 100644 model.example.yml create mode 100644 phases/example/__init__.py rename phases/{ => example}/dust.py (92%) create mode 100644 processes/basic.py delete mode 100644 processes/registry.py create mode 100644 processes/star_formation.py delete mode 100644 processes/star_formation/__init__.py delete mode 100644 processes/star_formation/__pycache__/__init__.cpython-36.pyc delete mode 100644 processes/star_formation/__pycache__/constant_efficiency.cpython-36.pyc delete mode 100644 settings.example.yml diff --git a/.gitignore b/.gitignore index 196888f..54b37fd 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,6 @@ +# SPiCE-specific stuff +model.yml + # Byte-compiled / optimized / DLL files __pycache__/ *.py[cod] diff --git a/SPiCE.py b/SPiCE.py index a6d644f..84dd4d5 100755 --- a/SPiCE.py +++ b/SPiCE.py @@ -1,6 +1,8 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- """ +Self-consistent Photo-ionisation and Chemical Evolution (SPiCE) + Created on Tue Mar 26 11:56:10 2019 @author: yago @@ -8,53 +10,70 @@ from __future__ import print_function, division -import argparse, yaml +import argparse +import yaml -import SPiCE import settings import phases -from phases import select_phase import processes __version__ = "0.0.1-alpha" + +def get_class(base_module, class_name): + selected_class = base_module + for name in class_name.split('.'): + selected_class = selected_class.__dict__[name] + return selected_class + + # ----------------------------------------------------------------------------- class Model(phases.basic.MultiphaseMedium): - def __init__(self, context): + def read_config_file(self, name): + print(name) + with open(name, "r") as settings_file: + user_settings = yaml.safe_load(settings_file) + self.context = settings.validate(user_settings) + + def __init__(self, config_file): + if config_file is None: + config_file = 'model.yml' + self.read_config_file(config_file) + self.phases = {} - for phase_name in context['phases'].keys(): - phase = context['phases'][phase_name] - self.phases[phase_name] = select_phase(phase['type'])(phase['params']) + for phase_name in self.context['phases']: + phase = self.context['phases'][phase_name] + self.phases[phase_name] = phases.registry[phase['type']](self, phase['params']) self.processes = {} - for process_type in context['processes'].keys(): - process = context['processes'][process_type] - process_name = process['name'] - self.processes[process_type] = processes.registry.select_process(process_type, process_name) - self.processes[process_type].init(process['params']) + for process_name in self.context['processes'].keys(): + process = self.context['processes'][process_name] + self.processes[process_name] = processes.registry[process['type']](self, process['params']) + + self.integrator = self.context['integrator'] + +# self.run() + def update_derivatives(self, term): print("This should not happen!") raise(-1) -def main(): + + +if __name__ == "__main__": parser = argparse.ArgumentParser( prog="SPiCE", description="Command line tools for the SPiCE chemical evolution code.", ) - parser.add_argument("-v", "--version", action="version", version=SPiCE.__version__) + parser.add_argument("-v", "--version", action="version", version=__version__) parser.add_argument("--config", metavar="FILENAME", help="configuration file to use containing model initial params") args = parser.parse_args() - user_settings = {} - if args.config != None : user_settings = read_config_file(args.config) - - context = settings.validate(user_settings) - - model = Model(context) + model = Model(args.config) print("\nMasses:") for phase in model.phases.keys(): @@ -63,12 +82,3 @@ def main(): print("\nProcesses:") for process in model.processes.keys(): print(' ', process, model.processes[process].tau) - -def read_config_file(name): - with open(name, "r") as settings_file: - user_settings = yaml.safe_load(settings_file) - return user_settings - - -if __name__ == "__main__": - main() diff --git a/model.example.yml b/model.example.yml new file mode 100644 index 0000000..52e29f5 --- /dev/null +++ b/model.example.yml @@ -0,0 +1,26 @@ +imf: kroupa +binary_star_rates: 0.40 +dtd_sn: rlp + +integrator: + relative_accuracy: 1e-6 + initial_time_Gyr: 0. + final_time_Gyr: 13.7 + +phases: + gas: + type: basic + params: + initial_mass_Msun: 2.34 + stars: + type: basic + params: + initial_mass_Msun: 0. + +processes: + star_formation: + type: timescale_constant + params: + input_phase: gas + output_phase: stars + timescale_Gyr: 3.478 diff --git a/phases/__init__.py b/phases/__init__.py index 71f00df..9ac0973 100644 --- a/phases/__init__.py +++ b/phases/__init__.py @@ -1,15 +1,15 @@ from . import basic from . import GasSpecies -from .dust import Dust +from . import example -def select_phase(name): - available_phases = { - "simple": basic.Phase, +registry = { + "basic": basic.Phase, "multiphase": basic.MultiphaseMedium, "ionized_hydrogen": GasSpecies.Ionized_Hydrogen, "neutral_hydrogen": GasSpecies.Neutral_Hydrogen, "molecular_hydrogen": GasSpecies.Molecular_Hydrogen, - "dust": Dust, + "dust": example.dust.Dust, } return available_phases[name] + diff --git a/phases/basic.py b/phases/basic.py index 9743509..cf13ec3 100644 --- a/phases/basic.py +++ b/phases/basic.py @@ -8,31 +8,43 @@ from __future__ import print_function, division +from astropy import units as u + # ----------------------------------------------------------------------------- class Phase: - def __init__(self, params): + def __init__(self, model, params): + self.model = model self.params = {**self.default_settings(), **params} + self.current_mass = self.params['initial_mass_Msun']*u.Msun + self.dm_dt = self.params['dm_dt_Msun_yr'] * u.Msun/u.yr + + def default_settings(self): + return { + 'initial_mass_Msun': 0., + 'dm_dt_Msun_yr': 0. + } def mass(self): - return self.params['current_mass'] + return self.current_mass def update_derivatives(self, term): - self.params['dm_dt'] += term + self.dm_dt += term - def default_settings(self): - return { - 'current_mass': 0., - 'dm_dt': 0. - } + def get_timestep(self): + self.dm_dt += self.model.integrator['relative_accuracy']*self.current_mass/self.dm_dt + + def update_mass(self, timestep): + self.current_mass += self.dm_dt * timestep # ----------------------------------------------------------------------------- class MultiphaseMedium(Phase): - def __init__(self, phase_dict={}): - self.phases = phase_dict + def __init__(self, model, params): +# self.phases = phase_dict + raise("TO DO: Multiphase medium to be implemented") def mass(self): m = 0. diff --git a/phases/example/__init__.py b/phases/example/__init__.py new file mode 100644 index 0000000..99d379e --- /dev/null +++ b/phases/example/__init__.py @@ -0,0 +1 @@ +from . import dust diff --git a/phases/dust.py b/phases/example/dust.py similarity index 92% rename from phases/dust.py rename to phases/example/dust.py index e85f6c1..520557c 100644 --- a/phases/dust.py +++ b/phases/example/dust.py @@ -1,4 +1,4 @@ -from . import basic +from .. import basic class Dust(basic.Phase): diff --git a/processes/__init__.py b/processes/__init__.py index 7246323..ab526f2 100644 --- a/processes/__init__.py +++ b/processes/__init__.py @@ -1 +1,5 @@ -from . import registry +from . import basic + +registry = { + "timescale_constant": basic.Constant_timescale, + } diff --git a/processes/basic.py b/processes/basic.py new file mode 100644 index 0000000..c267680 --- /dev/null +++ b/processes/basic.py @@ -0,0 +1,30 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Jun 26 08:39:06 2019 + +@author: yago +""" + +from astropy import units as u + + +# ----------------------------------------------------------------------------- +class Process(object): + + def compute_derivatives(): + raise("ERROR: compute_derivatives method not implemented") + + +# ----------------------------------------------------------------------------- +class Constant_timescale(Process): + + def __init__(self, model, params): + self.input = model.phases[params['input_phase']] + self.output = model.phases[params['output_phase']] + self.tau = params['timescale_Gyr']*u.Gyr + + def compute_derivatives(self): + flux = self.input.mass()/self.tau + self.input.update_derivatives(-flux) + self.output.update_derivatives(flux) diff --git a/processes/registry.py b/processes/registry.py deleted file mode 100644 index 9a20242..0000000 --- a/processes/registry.py +++ /dev/null @@ -1,10 +0,0 @@ -from . import star_formation - -available_processes = { - "star_formation": { - "constant": star_formation.constant_efficiency, - } -} - -def select_process(process_type, selection): - return available_processes[process_type][selection] \ No newline at end of file diff --git a/processes/star_formation.py b/processes/star_formation.py new file mode 100644 index 0000000..2ba27fd --- /dev/null +++ b/processes/star_formation.py @@ -0,0 +1,23 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Mar 27 10:30:28 2019 + +@author: yago +""" + +from . import basic + + +# ----------------------------------------------------------------------------- +class Constant_efficiency(basic.Process): + + def __init__(self, gas, stars, tau_SF): + self.gas = gas # gas phase that is converted into stars + self.tau = tau_SF # gas depletion time, in Gyr + self.stars = stars # destination stellar phase + + def compute_derivatives(self): + SFR = self.gas.mass()/self.tau + self.gas.update_derivatives(-SFR) + self.stars.update_derivatives(SFR) diff --git a/processes/star_formation/__init__.py b/processes/star_formation/__init__.py deleted file mode 100644 index e055d96..0000000 --- a/processes/star_formation/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from . import constant_efficiency diff --git a/processes/star_formation/__pycache__/__init__.cpython-36.pyc b/processes/star_formation/__pycache__/__init__.cpython-36.pyc deleted file mode 100644 index 7c31b0c7593b41cc56c69bdb36f6cb5150bd67e7..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 193 zcmXr!<>fMATpy#vz`*brh~a<<$Z`PUVgVqL!jQt4!;s4u#mLBz!W7J)$^4QLD6GkN zOE@_{uec;JuOvP-EiE%SGc_-{(od7=7DExpfRzkIEIcSdy8aryn1mnU`4-AFo$Xd5gm) RH$SB`C)EyQeK8O-007aKG8+H@ diff --git a/processes/star_formation/__pycache__/constant_efficiency.cpython-36.pyc b/processes/star_formation/__pycache__/constant_efficiency.cpython-36.pyc deleted file mode 100644 index f80917948ad821eabdbf2aab2e1b6880b32e9028..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 566 zcmY*V%}&EG40h6Xn+iMdbAv;hcR%%2&de8hcV?o572wO zF@m`Hau^PjU`FsrmeDP`E2g6EQMBLb9dvr#f$<6CTz+lCC|0O>Pb z`vucld01MnS2K?1Ou10`97)ZEt^>E#XpsTH2c_*SB?ANM1dOGrOy6=2PEm1wTx7ff4`! diff --git a/settings.example.yml b/settings.example.yml deleted file mode 100644 index 7e80414..0000000 --- a/settings.example.yml +++ /dev/null @@ -1,41 +0,0 @@ -initial_gas_mass: 1.0 -initial_stellar_mass: 0.0 -tau_SF: 1 -initial_time: 0.0 -final_time: 13.7 -integrator_relative_accuracy: 1e-6 -imf: kroupa -binary_star_rates: 0.40 -dtd_sn: rlp - -phases: - gas: - type: simple #NOMBRE DE UNA CLASE - params: - mass: 2.34 - type: ionized_hydrogen - params: - mass: 2.34 - temperature: 1.0e6 - pressure: 1.0e4 - type: neutral_hydrogen - params: - mass: 2.34 - temperature: 1.0e4 - pressure: 1.0e4 - type: molecular_hydrogen - params: - mass: 2.34 - temperature: 1.0e2 - pressure: 1.0e4 - dust: - type: dust - params: - sticky_coeficient: 0.6254 - molecular_cloud_density: 5671 - -processes: - star_formation: - name: constant - params: - tau_SF: 3.478 \ No newline at end of file From 0a599d436efe831a20d3541b3ad99ad4aa144e3e Mon Sep 17 00:00:00 2001 From: Yago Date: Wed, 26 Jun 2019 14:14:54 +0200 Subject: [PATCH 08/20] Store mass history as 'mass' and add 'current_mass' method --- SPiCE.py | 7 ++++++- phases/basic.py | 17 ++++++++++------- processes/basic.py | 2 +- 3 files changed, 17 insertions(+), 9 deletions(-) diff --git a/SPiCE.py b/SPiCE.py index 84dd4d5..373ffac 100755 --- a/SPiCE.py +++ b/SPiCE.py @@ -12,11 +12,13 @@ import argparse import yaml +from astropy import units as u import settings import phases import processes + __version__ = "0.0.1-alpha" @@ -53,7 +55,10 @@ def __init__(self, config_file): self.integrator = self.context['integrator'] -# self.run() + self.run() + + def run(self): + self.time = [self.integrator['initial_time_Gyr']*u.Gyr] def update_derivatives(self, term): diff --git a/phases/basic.py b/phases/basic.py index cf13ec3..e997e09 100644 --- a/phases/basic.py +++ b/phases/basic.py @@ -8,16 +8,16 @@ from __future__ import print_function, division +import numpy as np from astropy import units as u - # ----------------------------------------------------------------------------- class Phase: def __init__(self, model, params): self.model = model self.params = {**self.default_settings(), **params} - self.current_mass = self.params['initial_mass_Msun']*u.Msun + self.mass_history = [self.params['initial_mass_Msun']*u.Msun] self.dm_dt = self.params['dm_dt_Msun_yr'] * u.Msun/u.yr def default_settings(self): @@ -26,8 +26,11 @@ def default_settings(self): 'dm_dt_Msun_yr': 0. } - def mass(self): - return self.current_mass + def current_mass(self): + return self.mass_history[-1] + + def mass(self, t): + return np.interp(t, self.model.time, self.mass_history[-1]) def update_derivatives(self, term): self.dm_dt += term @@ -49,14 +52,14 @@ def __init__(self, model, params): def mass(self): m = 0. for phase in self.phases.values(): - m += phase.mass() + m += phase.current_mass() return m def m(self, phase='total'): if(phase == 'total'): - return self.mass() + return self.current_mass() else: - return self.phases[phase].mass() + return self.phases[phase].current_mass() def update_derivatives(self, term): print("TO DO: estimate mass fractions") diff --git a/processes/basic.py b/processes/basic.py index c267680..5111073 100644 --- a/processes/basic.py +++ b/processes/basic.py @@ -25,6 +25,6 @@ def __init__(self, model, params): self.tau = params['timescale_Gyr']*u.Gyr def compute_derivatives(self): - flux = self.input.mass()/self.tau + flux = self.input.current_mass()/self.tau self.input.update_derivatives(-flux) self.output.update_derivatives(flux) From e6b56686470b86fdc893c04a4ca171d06844ed12 Mon Sep 17 00:00:00 2001 From: Yago Date: Wed, 26 Jun 2019 15:39:48 +0200 Subject: [PATCH 09/20] First succesful SPiCE run git add .! :tada: --- SPiCE.py | 28 +++++++++++++++++++++++----- phases/basic.py | 29 ++++++++++++++++------------- processes/basic.py | 6 ++---- settings.py | 14 +++++++------- 4 files changed, 48 insertions(+), 29 deletions(-) diff --git a/SPiCE.py b/SPiCE.py index 373ffac..848d1b0 100755 --- a/SPiCE.py +++ b/SPiCE.py @@ -12,7 +12,9 @@ import argparse import yaml -from astropy import units as u +import numpy as np + +from matplotlib import pyplot as plt import settings import phases @@ -54,12 +56,25 @@ def __init__(self, config_file): self.processes[process_name] = processes.registry[process['type']](self, process['params']) self.integrator = self.context['integrator'] - self.run() def run(self): - self.time = [self.integrator['initial_time_Gyr']*u.Gyr] - + print(self.integrator) + self.time_Gyr = [self.integrator['initial_time_Gyr']] + while self.time_Gyr[-1] < self.integrator['final_time_Gyr']: + current_time_Gyr = self.time_Gyr[-1] + for phase in self.phases.values(): + phase.reset_timestep() + for process in self.processes.values(): + process.compute_derivatives() + timesteps_Gyr = [self.integrator['final_time_Gyr'] - current_time_Gyr] + for phase in self.phases.values(): + timesteps_Gyr.append(phase.get_timestep_Gyr()) + timestep_Gyr = np.nanmax([np.nanmin(timesteps_Gyr), + self.integrator['minimum_timestep_Gyr']]) + for phase in self.phases.values(): + phase.update_mass(timestep_Gyr) + self.time_Gyr.append(current_time_Gyr + timestep_Gyr) def update_derivatives(self, term): print("This should not happen!") @@ -86,4 +101,7 @@ def update_derivatives(self, term): print("\nProcesses:") for process in model.processes.keys(): - print(' ', process, model.processes[process].tau) + print(' ', process, model.processes[process].tau_Gyr) + +# plt.plot(model.time_Gyr, model['gas'].mass_history_Msun, 'b-') +# plt.plot(model.time_Gyr, model['stars'].mass_history_Msun, 'r-') diff --git a/phases/basic.py b/phases/basic.py index e997e09..e35ac3c 100644 --- a/phases/basic.py +++ b/phases/basic.py @@ -9,7 +9,7 @@ from __future__ import print_function, division import numpy as np -from astropy import units as u + # ----------------------------------------------------------------------------- class Phase: @@ -17,29 +17,32 @@ class Phase: def __init__(self, model, params): self.model = model self.params = {**self.default_settings(), **params} - self.mass_history = [self.params['initial_mass_Msun']*u.Msun] - self.dm_dt = self.params['dm_dt_Msun_yr'] * u.Msun/u.yr + self.mass_history_Msun = [float(self.params['initial_mass_Msun'])] def default_settings(self): return { 'initial_mass_Msun': 0., - 'dm_dt_Msun_yr': 0. } - def current_mass(self): - return self.mass_history[-1] + def current_mass_Msun(self): + return self.mass_history_Msun[-1] def mass(self, t): return np.interp(t, self.model.time, self.mass_history[-1]) + def reset_timestep(self): + self.dm_dt_Msun_Gyr = 0. + def update_derivatives(self, term): - self.dm_dt += term + self.dm_dt_Msun_Gyr += term - def get_timestep(self): - self.dm_dt += self.model.integrator['relative_accuracy']*self.current_mass/self.dm_dt + def get_timestep_Gyr(self): + return np.abs((self.model.integrator['relative_accuracy'] + * self.current_mass_Msun() / self.dm_dt_Msun_Gyr)) - def update_mass(self, timestep): - self.current_mass += self.dm_dt * timestep + def update_mass(self, timestep_Gyr): + self.mass_history_Msun.append(self.current_mass_Msun() + + self.dm_dt_Msun_Gyr*timestep_Gyr) # ----------------------------------------------------------------------------- @@ -52,14 +55,14 @@ def __init__(self, model, params): def mass(self): m = 0. for phase in self.phases.values(): - m += phase.current_mass() + m += phase.current_mass_Msun() return m def m(self, phase='total'): if(phase == 'total'): return self.current_mass() else: - return self.phases[phase].current_mass() + return self.phases[phase].current_mass_Msun() def update_derivatives(self, term): print("TO DO: estimate mass fractions") diff --git a/processes/basic.py b/processes/basic.py index 5111073..cd87e46 100644 --- a/processes/basic.py +++ b/processes/basic.py @@ -6,8 +6,6 @@ @author: yago """ -from astropy import units as u - # ----------------------------------------------------------------------------- class Process(object): @@ -22,9 +20,9 @@ class Constant_timescale(Process): def __init__(self, model, params): self.input = model.phases[params['input_phase']] self.output = model.phases[params['output_phase']] - self.tau = params['timescale_Gyr']*u.Gyr + self.tau_Gyr = float(params['timescale_Gyr']) def compute_derivatives(self): - flux = self.input.current_mass()/self.tau + flux = self.input.current_mass_Msun()/self.tau_Gyr self.input.update_derivatives(-flux) self.output.update_derivatives(flux) diff --git a/settings.py b/settings.py index fed87e3..caf2ff0 100644 --- a/settings.py +++ b/settings.py @@ -6,15 +6,15 @@ """ default = { - 'initial_gas_mass': 1.0, # in solar masses - 'initial_stellar_mass': 0.0, # in solar masses - 'tau_SF': 1, # in Gyrs - 'initial_time': 0.0, # in Gyr - 'final_time': 13.7, # inGyr - 'integrator_relative_accuracy': 1e-6, # Maximum error accepted for the integration 'imf': "kroupa", # Initial mass function 'binary_star_rates': 0.40, # Fraction of binary systems - 'dtd_sn': "rlp", # Delay Time Distribution for supernova rates + 'dtd_sn': "rlp", + 'integrator': { + 'initial_time_Gyr': 0.0, + 'final_time_Gyr': 13.7, + 'relative_accuracy': 1e-6, + 'minimum_timestep_Gyr': 1e-6 + }, 'phases': {}, 'processes': {} } From ff5694f9ae5b8a8b2ea25643677573e69a558067 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juanjo=20Baz=C3=A1n?= Date: Thu, 11 Jul 2019 11:36:54 +0200 Subject: [PATCH 10/20] rename example file --- model.example.yml => model.yml.example | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename model.example.yml => model.yml.example (100%) diff --git a/model.example.yml b/model.yml.example similarity index 100% rename from model.example.yml rename to model.yml.example From cb7de0a51b71a8406bb46b0c75bffdcb6f4854de Mon Sep 17 00:00:00 2001 From: Mario Romero Date: Wed, 26 Jun 2019 12:44:02 +0200 Subject: [PATCH 11/20] Cambios recientes --- phases/GasGenerics.py | 6 +++--- phases/__init__.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/phases/GasGenerics.py b/phases/GasGenerics.py index 301f614..1f92e45 100644 --- a/phases/GasGenerics.py +++ b/phases/GasGenerics.py @@ -20,8 +20,8 @@ class Gas(basic.MultiphaseMedium): #DEFAULT SETTINGS #--------------------- def __init__(self, params): - self.params = {**self.default_settings(), **params} #Falta 'model' - self._state(P=self.params["pressure"],T=self.params["temperature"],M=self.params["mass"]) + self.params = {**self.default_settings(), **params} + self._state() #(!) Mirar si hay que meter parametros o no def default_settings(self): return{ @@ -257,4 +257,4 @@ def state(self,P=-1,T=-1,M=-1,N=-1):#inputs: pressure, temperature, number of pa #Compute the other variables - self._compute_variables_indepedent_of_state() \ No newline at end of file + self._compute_variables_indepedent_of_state() diff --git a/phases/__init__.py b/phases/__init__.py index 9ac0973..d829b00 100644 --- a/phases/__init__.py +++ b/phases/__init__.py @@ -1,6 +1,7 @@ from . import basic from . import GasSpecies from . import example +from . import GasSpecies registry = { "basic": basic.Phase, @@ -12,4 +13,3 @@ } return available_phases[name] - From aa30eb19d1efdb412ec8fbad4eebce88eca44a2e Mon Sep 17 00:00:00 2001 From: Mario Romero Date: Fri, 12 Jul 2019 15:31:08 +0200 Subject: [PATCH 12/20] some debugging --- SPiCE.py | 2 +- phases/GasGenerics.py | 6 +++--- phases/GasSpecies.py | 3 ++- phases/__init__.py | 4 +--- 4 files changed, 7 insertions(+), 8 deletions(-) diff --git a/SPiCE.py b/SPiCE.py index 848d1b0..6e4beeb 100755 --- a/SPiCE.py +++ b/SPiCE.py @@ -42,7 +42,7 @@ def read_config_file(self, name): def __init__(self, config_file): if config_file is None: - config_file = 'model.yml' + config_file = 'model.yml.example' self.read_config_file(config_file) self.phases = {} diff --git a/phases/GasGenerics.py b/phases/GasGenerics.py index 1f92e45..a0c705f 100644 --- a/phases/GasGenerics.py +++ b/phases/GasGenerics.py @@ -7,7 +7,7 @@ import numpy as np import astropy.units as u import astropy.constants as cte -import basic +from . import basic ''' PARENT CLASSES @@ -20,8 +20,8 @@ class Gas(basic.MultiphaseMedium): #DEFAULT SETTINGS #--------------------- def __init__(self, params): - self.params = {**self.default_settings(), **params} - self._state() #(!) Mirar si hay que meter parametros o no + self.params = {**self.default_settings(), **params} #Falta 'models', mirate el backup + self._state(P=self.params["pressure"],T=self.params["temperature"],M=self.params["mass"]) def default_settings(self): return{ diff --git a/phases/GasSpecies.py b/phases/GasSpecies.py index e4a4ea3..9937e0d 100644 --- a/phases/GasSpecies.py +++ b/phases/GasSpecies.py @@ -7,7 +7,8 @@ import numpy as np import astropy.units as u import astropy.constants as cte -import GasGenerics as gas +#import GasGenerics as gas +from . import GasGenerics as gas ''' diff --git a/phases/__init__.py b/phases/__init__.py index d829b00..ed039c5 100644 --- a/phases/__init__.py +++ b/phases/__init__.py @@ -1,7 +1,6 @@ from . import basic from . import GasSpecies from . import example -from . import GasSpecies registry = { "basic": basic.Phase, @@ -10,6 +9,5 @@ "neutral_hydrogen": GasSpecies.Neutral_Hydrogen, "molecular_hydrogen": GasSpecies.Molecular_Hydrogen, "dust": example.dust.Dust, - } - return available_phases[name] + } \ No newline at end of file From f7cb0cdf7fd2cdc8014270fa3330e61eedbbf65e Mon Sep 17 00:00:00 2001 From: Mario Romero Date: Mon, 15 Jul 2019 14:10:50 +0200 Subject: [PATCH 13/20] Debugging. Code succesfully runs --- model.yml.example | 4 +++- settings.py | 6 +++--- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/model.yml.example b/model.yml.example index 52e29f5..4596dcc 100644 --- a/model.yml.example +++ b/model.yml.example @@ -3,9 +3,11 @@ binary_star_rates: 0.40 dtd_sn: rlp integrator: - relative_accuracy: 1e-6 +#For some reason, you have to give ALL variables of 'integrator' defined in 'settings.py' to make the code work + relative_accuracy: 1.0e-6 initial_time_Gyr: 0. final_time_Gyr: 13.7 + minimum_timestep_Gyr: 1.0e-3 phases: gas: diff --git a/settings.py b/settings.py index caf2ff0..f81050c 100644 --- a/settings.py +++ b/settings.py @@ -12,8 +12,8 @@ 'integrator': { 'initial_time_Gyr': 0.0, 'final_time_Gyr': 13.7, - 'relative_accuracy': 1e-6, - 'minimum_timestep_Gyr': 1e-6 + 'relative_accuracy': 1.0e-6, + 'minimum_timestep_Gyr': 1.0e-6 }, 'phases': {}, 'processes': {} @@ -22,7 +22,7 @@ valid_values = { 'imf': ['salpeter', 'starburst', 'chabrier', 'ferrini', 'kroupa', 'miller_scalo', 'maschberger'], 'dtd_sn': ['rlp', 'mdvp'], # rlp = Ruiz-Lapuente, mdvp = Mannucci, Della Valle, Panagia (2006) - 'phases': ['ionized_hydrogen', 'neutral_hydrogen', 'molecular_hydrogen'], + #'phases': ['ionized_hydrogen', 'neutral_hydrogen', 'molecular_hydrogen'], } def validate(params): From 48b0322a15995bf2e986598e603295ab5a91433e Mon Sep 17 00:00:00 2001 From: Mario Romero Date: Mon, 22 Jul 2019 12:16:34 +0200 Subject: [PATCH 14/20] Code now compatible with master --- SPiCE.py | 2 +- gas_example.yml | 30 +++++ phases/GasGenerics.py | 254 +++++++++++++++++++---------------------- phases/GasSpecies.py | 47 +++----- phases/HowToUse_gas.py | 64 ----------- phases/basic.py | 2 +- settings.py | 2 +- 7 files changed, 161 insertions(+), 240 deletions(-) create mode 100644 gas_example.yml delete mode 100644 phases/HowToUse_gas.py diff --git a/SPiCE.py b/SPiCE.py index 6e4beeb..a5e9e47 100755 --- a/SPiCE.py +++ b/SPiCE.py @@ -42,7 +42,7 @@ def read_config_file(self, name): def __init__(self, config_file): if config_file is None: - config_file = 'model.yml.example' + config_file = 'gas_example.yml' self.read_config_file(config_file) self.phases = {} diff --git a/gas_example.yml b/gas_example.yml new file mode 100644 index 0000000..1770da7 --- /dev/null +++ b/gas_example.yml @@ -0,0 +1,30 @@ +imf: kroupa +binary_star_rates: 0.40 +dtd_sn: rlp + +integrator: +#For some reason, you have to give ALL variables of 'integrator' defined in 'settings.py' to make the code work + relative_accuracy: 1.0e-6 + initial_time_Gyr: 0. + final_time_Gyr: 13.7 + minimum_timestep_Gyr: 1.0e-3 + +phases: + gas: + type: molecular_hydrogen + params: + initial_mass_Msun: 2.34 + initial_temperature_K: 100. + initial_pressure_cgs: 1.38e-12 #erg/cm^3 + stars: + type: basic + params: + initial_mass_Msun: 0. + +processes: + star_formation: + type: timescale_constant + params: + input_phase: gas + output_phase: stars + timescale_Gyr: 3.478 diff --git a/phases/GasGenerics.py b/phases/GasGenerics.py index a0c705f..a0a8cab 100644 --- a/phases/GasGenerics.py +++ b/phases/GasGenerics.py @@ -13,101 +13,110 @@ PARENT CLASSES classes defined here should not be defined as objects in the main code. ''' -class Gas(basic.MultiphaseMedium): +class Gas(basic.Phase): #class Gas(object): #--------------------- #DEFAULT SETTINGS #--------------------- - def __init__(self, params): - self.params = {**self.default_settings(), **params} #Falta 'models', mirate el backup - self._state(P=self.params["pressure"],T=self.params["temperature"],M=self.params["mass"]) + def __init__(self, model, params): + #What everyone does + self.model = model + self.params = {**self.default_settings(), **params} + self.mass_history_Msun = [float(self.params['initial_mass_Msun'])] + self.temperature_history_K = [float(self.params['initial_temperature_K'])] + self.pressure_history_cgs = [float(self.params['initial_pressure_cgs'])] + + self._constants() + self._set() + self._state() + + def _constants(self): + #Useful conversions and constants (this will save us LOTS of parentesis, trust me!) + self._solMass_to_g = ((1*u.solMass).to(u.g)).value + self._kB = ((cte.k_B).to(u.erg/u.K)).value #Boltzmann constant in erg/K + self._hP = ((cte.h).to(u.erg*u.s)).value #Planck constant in erg*s + self._hbar = self._hP/(2.*np.pi) + self._c = ((cte.c).to(u.cm/u.s)).value #Speed of light + + def _set(self): + #All internal variables will be in cgs + self._temperature = self.temperature_history_K[-1] + self._pressure = self.pressure_history_cgs[-1] + self._gas_mass = self.mass_history_Msun[-1]*self._solMass_to_g + def current_temperature_K(self): + return self.temperature_history_K[-1] + def current_pressure_cgs(self): + return self.pressure_history_cgs[-1] + + def update_mass(self, timestep_Gyr): + #Actually, this should be called 'update_all', but is called this way to keep compatibility + self.mass_history_Msun.append(self.current_mass_Msun() + self.dm_dt_Msun_Gyr*timestep_Gyr) + + #If you don't want to store the pressure history because, for example, it does not change over time + # you only have to comment the relevant line of these two. + self.temperature_history_K.append(self.current_temperature_K()) + self.pressure_history_cgs.append(self.current_pressure_cgs()) + + self._set() + self._state() + def default_settings(self): return{ - 'mass':0.0*u.solMass, - 'dm_dt': 0.0*(u.solMass/u.yr), - - 'temperature': 0.0*u.K, - 'pressure':0.0*(u.erg/(u.cm*u.cm*u.cm)) + 'initial_mass_Msun':0.0, + 'initial_temperature_K': 0.0, + 'initial_pressure_cgs':0.0 } #--------------------- #COMMON ATRIBUTES/METHODS TO ALL DERIVED CLASSES #--------------------- - def _set(self): - - uV = (u.cm*u.cm*u.cm) #volume units - no_units = u.m/u.m #dimensionless units - #Internal variables - self._pressure = self.params["pressure"].to(u.erg/u.s) - self._temperature = self.params["temperature"].to(u.K) - self._number_particles = 0.0*no_units - self._thermal_energy = 0.0*u.erg - self._fugacity = 0.0*no_units - self._volume = 0.0*uV - self._gas_mass = self.params["mass"].to(u.solMass) - self._chemical_potential = 0.0*u.erg - self._number_density = 0.0 / uV - self._mass_density = 0.0*u.g/ uV - self._thermal_energy_density = 0.0*u.erg/uV - self._gamma = 0.0*no_units - def _compute_variables_indepedent_of_state(self): - uV = (u.cm*u.cm*u.cm) #volume units - no_units = u.m/u.m #dimensionless units - - self._gas_mass = (self._number_particles * self._particle_mass).to(u.solMass) - self._number_density = (self._number_particles / self._volume).to(1./uV) - self._mass_density = (self._number_density * self._particle_mass).to(u.g/uV) - self._thermal_energy_density = (self._thermal_energy / self._volume).to(u.erg/uV) - self._chemical_potential = (cte.k_B*self._temperature*np.log(self._fugacity)).to(u.erg) - self._gamma = (1.+(self._pressure/self._thermal_energy_density)).to(no_units) - #--------------------- - #CONSTRUCTORS - #--------------------- - #def __init__(self): - # print("Warning: '__init__()' code written for testing purposes!") - # self.default_settings() + self._number_density = self._number_particles / self._volume + self._mass_density = self._number_density * self._particle_mass_g + self._thermal_energy_density = self._thermal_energy / self._volume + self._chemical_potential = self._kB*self._temperature*np.log(self._fugacity) + self._gamma = 1.+(self._pressure/self._thermal_energy_density) #--------------------- #OUTPUTS #--------------------- - def mass(self,units=u.solMass): #Polymorph from Multiphase - return (self._gas_mass).to(units) + #def mass(...): comes from basic phase. So I do not redefine here def pressure(self,units=u.erg/(u.cm*u.cm*u.cm)): - return (self._pressure).to(units) + return ((self._pressure)*units).value def temperature(self,units=u.K): - return (self._temperature).to(units) + return ((self._temperature)*units).value def number_particles(self,units=u.m/u.m): - return (self._number_particles).to(units) + return ((self._number_particles)*units).value def thermal_energy(self,units=u.erg): - return (self._thermal_energy).to(units) + return ((self._thermal_energy)*units).value def fugacity(self,units=u.m/u.m): - return (self._fugacity).to(units) + return ((self._fugacity)*units).value def volume(self,units=(u.cm*u.cm*u.cm)): - return (self._volume).to(units) + return ((self._volume)*units).value def chemical_potential(self,units=u.erg): return (self._chemical_potential).to(units) def number_density(self,units=1./(u.cm*u.cm*u.cm)): - return (self._number_density).to(units) + return ((self._number_density)*units).value def mass_density(self,units=u.g/(u.cm*u.cm*u.cm)): - return (self._mass_density).to(units) + return ((self._mass_density)*units).value def thermal_energy_density(self,units=u.erg/(u.cm*u.cm*u.cm)): - return (self._thermal_energy_density(self)).to(units) + return ((self._thermal_energy_density(self))*(units)).value def adiabatic_constant(self,units=u.m/u.m): - return (self._gamma).to(units) + return ((self._gamma)*units).value #--------------------- #POLYMORPHIC METHODS #--------------------- #I.e.: These are the ONLY functions you have to redefine in each subclass - def state(self,P=-1,T=-1,M=-1,N=-1):#inputs: pressure, temperature, number of particles, thermal energy, mass + def _state(self,P=-1,T=-1,M=-1,N=-1):#inputs: pressure, temperature, number of particles, thermal energy, mass + #print("Warning: Using a parent class. Should not happen!") raise NameError("Called by parent class. Should not happen!") - def update_derivatives(self,term): - raise NameError("Called by parent class. Should not happen!") + #def update_derivatives(self,term): #Comes from basic + # raise NameError("Called by parent class. Should not happen!") ''' CHILDS OF GAS @@ -118,33 +127,15 @@ class Monoatomic(Gas): #--------------------- #STATE METHOD #--------------------- - def state(self,P=-1,T=-1,M=-1,N=-1):#inputs: pressure, temperature, number of particles, thermal energy, mass - uV = (u.cm*u.cm*u.cm) #volume units - no_units = u.m/u.m #dimensionless units - - #Convert mass, if given - if M != -1: - N = (M.to(u.g))/(self._particle_mass).to(u.g) - - #You MUST give pressure, temperature and number of particles (or total mass) - if (P == -1) or (T == -1) or (N == -1): - raise NameError("Gas._state() bad initialized! You must give pressure, temperature and total mass/number or particles!") - else: - #Initialize thermodynamic variables - self._pressure = P.to(u.erg/uV) - self._number_particles = N.to(no_units) - - self._temperature = T.to(u.K) - - #We compute energy first - self._thermal_energy = (1.5*self._number_particles * (cte.k_B) * self._temperature).to(u.erg) - - #Volume... - self._volume = ((2.*self._thermal_energy) / (3.*self._pressure)).to(uV) + def _state(self):#inputs: pressure, temperature, number of particles, thermal energy, mass + #Here we compute the other thermodynamic variables + + self._number_particles = (self._gas_mass)/(self._particle_mass_g) + self._thermal_energy = 1.5*self._number_particles * self._kB * self._temperature + self._volume = (2.*self._thermal_energy) / (3.*self._pressure) - #And fugacity - self._fugacity = (self._pressure * (2.*np.pi*self._particle_mass / (cte.h*cte.h) )**(-1.5) * (cte.k_B*self._temperature)**(-2.5)).to(no_units) - + #Fugacity = exp(chemical_potential / kT) . It is another way to describe the chemical potential + self._fugacity = self._pressure * (2.*np.pi*self._particle_mass_g / (self._hP*self._hP) )**(-1.5) * (self._kB*self._temperature)**(-2.5) #Compute the other variables self._compute_variables_indepedent_of_state() @@ -154,35 +145,33 @@ class Diatomic(Gas): #AUXILIARY METHODS #--------------------- def _molecule_mass(self): - return (self._atom_mass[0]+self._atom_mass[1]).to(u.g) + return self._atom_mass_g[0]+self._atom_mass_g[1] #ESTOY AQUÍ def _moment_of_inertia(self): #Get reduced mass - - m_r = ( self._atom_mass[0]*self._atom_mass[1] / (self._atom_mass[0]+self._atom_mass[1]) ).to(u.g) - return m_r*self._atom_distance*self._atom_distance + m_r = self._atom_mass_g[0]*self._atom_mass_g[1] / (self._atom_mass_g[0]+self._atom_mass_g[1]) + return m_r*self._atom_distance_cm*self._atom_distance_cm def _angular_frequency(self): - return (self._wavenumber * cte.c ).to(1./u.s) + return self._wavenumber_cm_minus1 * self._c #--------------------- #PARTITION FUNCTION METHODS #--------------------- def _Z_rot(self): - - no_units = u.cm/u.cm + #We compute Z for only rotations + #Compute a parameter to check if quantum effects are important - theta = ( 0.5 * cte.hbar*cte.hbar / (cte.k_B * self._temperature * self._moment_of_inertia()) ).to(no_units) - + theta = 0.5 * self._hbar*self._hbar / (self._kB * self._temperature * self._moment_of_inertia()) #Now do stuff Z = None - dZ_db = None #Partial derivative of Z with respect to beta = 1./kT + dZ_db = None #Partial derivative of Z with respect to beta = 1./kT (units: erg) if theta <= 1: #All rotational modes are enabled (=classical) Z = 1./theta - dZ_db = - cte.k_B * self._temperature / theta + dZ_db = - self._kB * self._temperature / theta else: #All rotational modes are disabled - Z = 1.0*no_units - dZ_db = 0.0*u.erg + Z = 1.0 + dZ_db = 0.0 return [Z , dZ_db] ''' @@ -192,69 +181,56 @@ def _Z_rot(self): ''' def _Z_vib(self): + #We compute Z for only molecule vibrations - no_units = u.cm/u.cm + #no_units = u.cm/u.cm #Compute a parameter to check if quantum effects are important - xi = ( cte.hbar * self._angular_frequency() / (cte.k_B * self._temperature) ).to(no_units) - + xi = self._hbar * self._angular_frequency() / (self._kB * self._temperature) #Now do stuff Z = None - dZ_db = None #Partial derivative of Z with respect to beta = 1./kT + dZ_db = None #Partial derivative of Z with respect to beta = 1./kT (units: erg) if xi > 1: #All vibrational modes are disabled - Z = 1.0*no_units - dZ_db = 0.0*u.erg + Z = 1.0 + dZ_db = 0.0 else: Z = 1./xi - dZ_db = - cte.k_B * self._temperature / xi + dZ_db = - self._kB * self._temperature / xi return [Z , dZ_db] ''' Again, this is a simplification. The exact method has Z = sum exp(-n*xi) = 1./(1-exp(-xi)) . - Note that the zero of energies is set as all vibrations at the fundamental state of the harmonic oscillator. + Also note that the zero of energies is set as all vibrations at the fundamental state of the harmonic oscillator. ''' #--------------------- #STATE METHOD #--------------------- - def state(self,P=-1,T=-1,M=-1,N=-1):#inputs: pressure, temperature, number of particles, mass + def _state(self):#inputs: pressure, temperature, number of particles, mass - uV = (u.cm*u.cm*u.cm) #volume units - no_units = u.m/u.m #dimensionless units + #uV = (u.cm*u.cm*u.cm) #volume units + #no_units = u.m/u.m #dimensionless units #Because you give each atom mass individually, you need to construct the particle mass first. - self._particle_mass = self._molecule_mass() + self._particle_mass_g = self._molecule_mass() + #Get the total number of particles + self._number_particles = (self._gas_mass)/(self._particle_mass_g) + #Let's define an useful variable + NkT = self._number_particles * (self._kB) * self._temperature - #Convert mass, if given - if M != -1: - N = (M.to(u.g))/(self._particle_mass).to(u.g) + #Partition functions (they are arrays, [0] is the Z itself, and [1] the derivative with respect to beta) + Zrot = self._Z_rot() + Zvib = self._Z_vib() + + #Now we compute the remaining variables + self._volume = NkT / self._pressure + #Fugacity = exp(chemical_potential / kT) . It is another way to describe the chemical potential + self._fugacity = self._pressure * (2.*np.pi*self._particle_mass_g / (self._hP*self._hP) )**(-1.5) * (self._kB*self._temperature)**(-2.5) / (Zrot[0]*Zvib[0]) + #Thermal energy + #const_factor = 1.5 - (Zrot[1]/( Zrot[0] * self._kB * self._temperature )) - (Zvib[1]/( Zvib[0] * self._kB * self._temperature )) + const_factor = 1.5 #Translational part + const_factor -= (Zrot[1]/( Zrot[0] * self._kB * self._temperature )) #Rotational part + const_factor -= (Zvib[1]/( Zvib[0] * self._kB * self._temperature )) #Vibrational part + self._thermal_energy = NkT * const_factor - #You MUST give pressure, temperature and number of particles (or total mass) - if (P == -1) or (T == -1) or (N == -1): - raise NameError("Gas._state() bad initialized! You must give pressure, temperature and total mass/number or particles!") - else: - #Initialize thermodynamic variables - self._pressure = P.to(u.erg/uV) - self._number_particles = N.to(no_units) - self._temperature = T.to(u.K) - - #Useful variable - NkT = (self._number_particles * (cte.k_B).to(u.erg/u.K) * self._temperature).to(u.erg) - - #Partition functions (they are arrays, [0] is the Z itself, and [1] the derivative with respect to beta) - Zrot = self._Z_rot() - Zvib = self._Z_vib() - - #We compute volume first (it's the ideal gas relation) - self._volume = ( NkT / self._pressure ).to(uV) - - #Fugacity... - self._fugacity = (self._pressure * (2.*np.pi*self._particle_mass / (cte.h*cte.h) )**(-1.5) * (cte.k_B*self._temperature)**(-2.5) / (Zrot[0]*Zvib[0]) ).to(no_units) - - #Thermal energy... - const_factor = 1.5 - (Zrot[1]/( Zrot[0] * cte.k_B * self._temperature )).to(no_units) - (Zvib[1]/( Zvib[0] * cte.k_B * self._temperature )).to(no_units) - self._thermal_energy = NkT * const_factor - - - #Compute the other variables self._compute_variables_indepedent_of_state() diff --git a/phases/GasSpecies.py b/phases/GasSpecies.py index 9937e0d..efa9aa5 100644 --- a/phases/GasSpecies.py +++ b/phases/GasSpecies.py @@ -4,12 +4,17 @@ Mario Romero May 2019 ''' -import numpy as np +#import numpy as np import astropy.units as u import astropy.constants as cte -#import GasGenerics as gas -from . import GasGenerics as gas +import GasGenerics as gas +''' +SOME CONSTANTS +''' + +atomic_mass_g = ((cte.u).to(u.g)).value +proton_mass_g = ((cte.m_p).to(u.g)).value ''' SPECIFIC GASES @@ -21,42 +26,16 @@ #Required data = particle mass #-------------- class Neutral_Hydrogen(gas.Monoatomic): - - _particle_mass = (1.008*cte.u).to(u.g) #1.008*atomic_mass - - def update_derivatives(self,term): - print("Neutral_Hydrogen: Work in Progress...") - - def evolve(self,dt): - print("Neutral_Hydrogen: 'evolve()' written for testing purposes!") - self.state(P=self._pressure,T=self._temperature,M=self._gas_mass) - + _particle_mass_g = 1.008*atomic_mass_g #1.008*atomic_mass class Ionized_Hydrogen(gas.Monoatomic): - - _particle_mass = (cte.m_p).to(u.g) #Proton mass - - def update_derivatives(self,term): - print("Ionized_Hydrogen: Work in Progress...") - - def evolve(self,dt): - print("Ionized_Hydrogen: 'evolve()' written for testing purposes!") - self.state(P=self._pressure,T=self._temperature,M=self._gas_mass) + _particle_mass_g = proton_mass_g #Proton mass #-------------- #DIATOMIC SPECIES #Needed data: particle masses (array), distance between particles and their fundamental vibration mode. #-------------- class Molecular_Hydrogen(gas.Diatomic): - - _atom_mass = [(1.008*cte.u).to(u.g) , (1.008*cte.u).to(u.g)] #Two hydrogen atoms - _atom_distance = 74.14*1e-10 * u.cm #Taken from wikipedia - _wavenumber = 4342.0*(1./u.cm) #wavenumber (k) of fundamental vibration. Taken from here: https://www.chem.purdue.edu/gchelp/vibs/h2.html // - - - def update_derivatives(self,term): - print("Molecular_Hydrogen: Work in Progress...") - - def evolve(self,dt): - print("Molecular_Hydrogen: 'evolve()' written for testing purposes!") - self.state(P=self._pressure,T=self._temperature,M=self._gas_mass) + _atom_mass_g = [1.008*atomic_mass_g , 1.008*atomic_mass_g ] #Two hydrogen atoms + _atom_distance_cm = 74.14*1e-10 #Taken from wikipedia + _wavenumber_cm_minus1 = 4342.0 #wavenumber (k) of fundamental vibration. Taken from here: https://www.chem.purdue.edu/gchelp/vibs/h2.html // diff --git a/phases/HowToUse_gas.py b/phases/HowToUse_gas.py deleted file mode 100644 index 564cd66..0000000 --- a/phases/HowToUse_gas.py +++ /dev/null @@ -1,64 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -''' -How to use the gas species - -Mario Romero June 2019 -''' - -import astropy.units as u -import astropy.constants as cte -import GasSpecies as species - -''' -Each species has its own class defined in 'GasSpecies.py' -''' - -#(1) Let's declare auxiliary variables - -uV = u.cm*u.cm*u.cm #units -P0 = ((1.e4*u.K/uV)*cte.k_B).to(u.erg/uV) #Common pressure to all phases - - -#(2) Let's declare the initial conditions of hydrogen -#molecular phase -M0_m = 10.*u.solMass -T0_m = 10.*u.K - -#atomic phase -M0_a = 100.*u.solMass -T0_a = 100.*u.K - -#ionized phase -M0_i = 1000.*u.solMass -T0_i = 1000.*u.K - -#(3) Let's declare hydrogen objects -#(!) Cambiar esto. -H2 = species.Molecular_Hydrogen() -H2.state(P=P0,T=T0_m,M=M0_m) #Your initial conditions goes here - -HI = species.Neutral_Hydrogen() -HI.state(P=P0,T=T0_a,M=M0_a) - -HII = species.Ionized_Hydrogen() -HII.state(P=P0,T=T0_i,M=M0_i) - -#(4) Do stuff (I'm throwing arbitrary numbers) -dt = 1.0*u.yr -t = 0.0*u.yr -tf = 0.5*u.yr - -while(t < tf): - - #Get the derivatives - H2.update_derivatives([H2,HI,HII]) - HI.update_derivatives([H2,HI,HII]) - HII.update_derivatives([H2,HI,HII]) - - #Use your favourite integrator - H2.evolve(dt) - HI.evolve(dt) - HII.evolve(dt) - - t += dt diff --git a/phases/basic.py b/phases/basic.py index e35ac3c..262036c 100644 --- a/phases/basic.py +++ b/phases/basic.py @@ -28,7 +28,7 @@ def current_mass_Msun(self): return self.mass_history_Msun[-1] def mass(self, t): - return np.interp(t, self.model.time, self.mass_history[-1]) + return np.interp(t, self.model.time, self.mass_history_Msun[-1]) def reset_timestep(self): self.dm_dt_Msun_Gyr = 0. diff --git a/settings.py b/settings.py index f81050c..98f33fe 100644 --- a/settings.py +++ b/settings.py @@ -22,7 +22,7 @@ valid_values = { 'imf': ['salpeter', 'starburst', 'chabrier', 'ferrini', 'kroupa', 'miller_scalo', 'maschberger'], 'dtd_sn': ['rlp', 'mdvp'], # rlp = Ruiz-Lapuente, mdvp = Mannucci, Della Valle, Panagia (2006) - #'phases': ['ionized_hydrogen', 'neutral_hydrogen', 'molecular_hydrogen'], + #'phases': ['gas','ionized_hydrogen', 'neutral_hydrogen', 'molecular_hydrogen'], } def validate(params): From 0ded34a7e96681e38e6c67ae0ccbfc1d607c0f62 Mon Sep 17 00:00:00 2001 From: Mario Romero Date: Mon, 22 Jul 2019 16:29:31 +0200 Subject: [PATCH 15/20] First HI->H->H2->Stars successful run! :tada: --- SPiCE.py | 2 +- gas_example.yml | 34 ++++++++++++++++++--- phases/GasGenerics.py | 12 +++++--- phases/GasSpecies.py | 4 +-- phases/__init__.py | 1 + phases/basic.py | 5 ++- processes/MoleculeFormation.py | 56 ++++++++++++++++++++++++++++++++++ processes/Recombination.py | 51 +++++++++++++++++++++++++++++++ processes/__init__.py | 4 +++ 9 files changed, 156 insertions(+), 13 deletions(-) create mode 100644 processes/MoleculeFormation.py create mode 100644 processes/Recombination.py diff --git a/SPiCE.py b/SPiCE.py index a5e9e47..aacba1e 100755 --- a/SPiCE.py +++ b/SPiCE.py @@ -21,7 +21,7 @@ import processes -__version__ = "0.0.1-alpha" +__version__ = "0.0.2-alpha" def get_class(base_module, class_name): diff --git a/gas_example.yml b/gas_example.yml index 1770da7..789716d 100644 --- a/gas_example.yml +++ b/gas_example.yml @@ -7,24 +7,48 @@ integrator: relative_accuracy: 1.0e-6 initial_time_Gyr: 0. final_time_Gyr: 13.7 - minimum_timestep_Gyr: 1.0e-3 + minimum_timestep_Gyr: 1.0e-4 phases: - gas: + gas_HI: + type: ionized_hydrogen + params: + initial_mass_Msun: 100. + initial_temperature_K: 1.0e4 + initial_pressure_cgs: 1.38e-12 #erg/cm^3 + gas_H: + type: neutral_hydrogen + params: + initial_mass_Msun: 0.0 + initial_temperature_K: 100. + initial_pressure_cgs: 1.38e-12 #erg/cm^3 + gas_H2: type: molecular_hydrogen params: - initial_mass_Msun: 2.34 + initial_mass_Msun: 0.0 initial_temperature_K: 100. initial_pressure_cgs: 1.38e-12 #erg/cm^3 stars: type: basic params: - initial_mass_Msun: 0. + initial_mass_Msun: 0.0 processes: + recombination: + type: H_recombination + params: + input_phase: gas_HI + output_phase: gas_H + molecule_formation: + type: H2_formation + params: + input_phase: gas_H + output_phase: gas_H2 + agent_phase: gas_HI + metallicity: 0.02 #Placeholder! one should put something like 'agent_phase2: dust' or something for better models star_formation: type: timescale_constant params: - input_phase: gas + input_phase: gas_H2 output_phase: stars timescale_Gyr: 3.478 diff --git a/phases/GasGenerics.py b/phases/GasGenerics.py index a0a8cab..5eb5fd9 100644 --- a/phases/GasGenerics.py +++ b/phases/GasGenerics.py @@ -14,7 +14,6 @@ classes defined here should not be defined as objects in the main code. ''' class Gas(basic.Phase): -#class Gas(object): #--------------------- #DEFAULT SETTINGS @@ -74,11 +73,16 @@ def default_settings(self): #--------------------- def _compute_variables_indepedent_of_state(self): - self._number_density = self._number_particles / self._volume + try: + self._number_density = self._number_particles / self._volume + self._thermal_energy_density = self._thermal_energy / self._volume + self._gamma = 1.+(self._pressure/self._thermal_energy_density) + except ZeroDivisionError: + self._number_density = 0.0 + self._thermal_energy_density = 0.0 + self._gamma = np.Infinity self._mass_density = self._number_density * self._particle_mass_g - self._thermal_energy_density = self._thermal_energy / self._volume self._chemical_potential = self._kB*self._temperature*np.log(self._fugacity) - self._gamma = 1.+(self._pressure/self._thermal_energy_density) #--------------------- #OUTPUTS diff --git a/phases/GasSpecies.py b/phases/GasSpecies.py index efa9aa5..c88a74f 100644 --- a/phases/GasSpecies.py +++ b/phases/GasSpecies.py @@ -7,7 +7,7 @@ #import numpy as np import astropy.units as u import astropy.constants as cte -import GasGenerics as gas +from . import GasGenerics as gas ''' SOME CONSTANTS @@ -26,7 +26,7 @@ #Required data = particle mass #-------------- class Neutral_Hydrogen(gas.Monoatomic): - _particle_mass_g = 1.008*atomic_mass_g #1.008*atomic_mass + _particle_mass_g = 1.008*atomic_mass_g class Ionized_Hydrogen(gas.Monoatomic): _particle_mass_g = proton_mass_g #Proton mass diff --git a/phases/__init__.py b/phases/__init__.py index ed039c5..1491a6b 100644 --- a/phases/__init__.py +++ b/phases/__init__.py @@ -1,4 +1,5 @@ from . import basic +#from . import GasGenerics from . import GasSpecies from . import example diff --git a/phases/basic.py b/phases/basic.py index 262036c..2689fbd 100644 --- a/phases/basic.py +++ b/phases/basic.py @@ -37,8 +37,11 @@ def update_derivatives(self, term): self.dm_dt_Msun_Gyr += term def get_timestep_Gyr(self): - return np.abs((self.model.integrator['relative_accuracy'] + try: + return np.abs((self.model.integrator['relative_accuracy'] * self.current_mass_Msun() / self.dm_dt_Msun_Gyr)) + except ZeroDivisionError: + return np.Infinity def update_mass(self, timestep_Gyr): self.mass_history_Msun.append(self.current_mass_Msun() diff --git a/processes/MoleculeFormation.py b/processes/MoleculeFormation.py new file mode 100644 index 0000000..fe050ed --- /dev/null +++ b/processes/MoleculeFormation.py @@ -0,0 +1,56 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Processes related with formation of molecules + +Mario Romero July 2019 +""" + +import numpy as np +import astropy.units as u +from . import basic + +''' +Hydrogen molecule creation +H -> H2 process +''' +class MolecularHydrogen_Formation(basic.Process): + + #--------------------- + #INIT + #--------------------- + def __init__(self, model, params): + self.input = model.phases[params['input_phase']] #Neutral Hydrogen + self.output = model.phases[params['output_phase']] #Molecular Hydrogen + self.agent = model.phases[params['agent_phase']] #Ionized Hydrogen + self._Z = float(params['metallicity']) + + self._constants() + self.tau_Gyr = self._formation_timescale() + + def compute_derivatives(self): + flux = self.input.current_mass_Msun()/self.tau_Gyr + self.input.update_derivatives(-flux) + self.output.update_derivatives(flux) + #--------------------- + #GETTING PHYSICAL PARAMETERS + #--------------------- + def _dust_density(self): + total_H_density = 0.5*self.agent.number_density() + self.input.number_density() + 2.*self.output.number_density() + return self._Z * total_H_density + def _formation_timescale(self): + #Taking Ascasibar+(in prep) formula again, ignoring the effective metallicity thing + mean_ov_cm3_s = 6.e-17 *(self.input.temperature()/100.)**(0.5) #Neutral H or H2 temperature???? + + tau_s = None #Not the correct tau unit-wise + try: + tau_s = 0.5/(mean_ov_cm3_s * self._dust_density()) + except ZeroDivisionError: + tau_s = np.Infinity + + return tau_s*self._s_to_Gyr #Correct units. + #--------------------- + #DEFINING USEFUL VARIABLES + #--------------------- + def _constants(self): + self._s_to_Gyr = ((1.*u.s).to(u.Gyr)).value \ No newline at end of file diff --git a/processes/Recombination.py b/processes/Recombination.py new file mode 100644 index 0000000..8a5b905 --- /dev/null +++ b/processes/Recombination.py @@ -0,0 +1,51 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Processes related with recombination + +Mario Romero July 2019 +""" + +import numpy as np +import astropy.units as u +from . import basic + +''' +Hydrogen Recombination +HI -> H process +''' +class Hydrogen_Recombination(basic.Process): + + #--------------------- + #INIT + #--------------------- + def __init__(self, model, params): + self.input = model.phases[params['input_phase']] #Ionized Hydrogen + self.output = model.phases[params['output_phase']] #Neutral Hydrogen + + self._constants() + self.tau_Gyr = self._recombination_timescale() + + def compute_derivatives(self): + flux = self.input.current_mass_Msun()/self.tau_Gyr + self.input.update_derivatives(-flux) + self.output.update_derivatives(flux) + #--------------------- + #GETTING THE TAU + #--------------------- + def _recombination_timescale(self): + #Taking Ascasibar+(in prep) formula + mean_ov_cm3_s = 4.1e-10 * (self.input.temperature())**(-0.8) #Case B recombination cross-section (Verner&Ferland 96) + + tau_s = None #Not the correct tau unit-wise + try: + tau_s = 2./(mean_ov_cm3_s * self.input.number_density()) + except ZeroDivisionError: + tau_s = np.Infinity + + return tau_s*self._s_to_Gyr #Correct units. + #--------------------- + #DEFINING USEFUL VARIABLES + #--------------------- + def _constants(self): + self._s_to_Gyr = ((1.*u.s).to(u.Gyr)).value \ No newline at end of file diff --git a/processes/__init__.py b/processes/__init__.py index ab526f2..d6256c6 100644 --- a/processes/__init__.py +++ b/processes/__init__.py @@ -1,5 +1,9 @@ from . import basic +from . import Recombination +from . import MoleculeFormation registry = { "timescale_constant": basic.Constant_timescale, + "H_recombination": Recombination.Hydrogen_Recombination, + "H2_formation": MoleculeFormation.MolecularHydrogen_Formation, } From ccc473cb2941f7b746744462f880ea760ff94cc7 Mon Sep 17 00:00:00 2001 From: Mario Romero Date: Tue, 30 Jul 2019 12:28:38 +0200 Subject: [PATCH 16/20] Hotfixes --- SPiCE.py | 1 + gas_example.yml | 37 ++++++++++++++++-------- processes/Ionization.py | 51 ++++++++++++++++++++++++++++++++++ processes/MoleculeFormation.py | 15 ++++++++-- processes/Recombination.py | 4 ++- processes/__init__.py | 3 ++ 6 files changed, 95 insertions(+), 16 deletions(-) create mode 100644 processes/Ionization.py diff --git a/SPiCE.py b/SPiCE.py index aacba1e..2d99ce6 100755 --- a/SPiCE.py +++ b/SPiCE.py @@ -75,6 +75,7 @@ def run(self): for phase in self.phases.values(): phase.update_mass(timestep_Gyr) self.time_Gyr.append(current_time_Gyr + timestep_Gyr) + #print(self.time_Gyr) def update_derivatives(self, term): print("This should not happen!") diff --git a/gas_example.yml b/gas_example.yml index 789716d..eee453f 100644 --- a/gas_example.yml +++ b/gas_example.yml @@ -3,7 +3,6 @@ binary_star_rates: 0.40 dtd_sn: rlp integrator: -#For some reason, you have to give ALL variables of 'integrator' defined in 'settings.py' to make the code work relative_accuracy: 1.0e-6 initial_time_Gyr: 0. final_time_Gyr: 13.7 @@ -13,20 +12,20 @@ phases: gas_HI: type: ionized_hydrogen params: - initial_mass_Msun: 100. + initial_mass_Msun: 99.8 initial_temperature_K: 1.0e4 initial_pressure_cgs: 1.38e-12 #erg/cm^3 gas_H: type: neutral_hydrogen params: - initial_mass_Msun: 0.0 - initial_temperature_K: 100. + initial_mass_Msun: 0.1 + initial_temperature_K: 1000. initial_pressure_cgs: 1.38e-12 #erg/cm^3 gas_H2: type: molecular_hydrogen params: - initial_mass_Msun: 0.0 - initial_temperature_K: 100. + initial_mass_Msun: 0.1 + initial_temperature_K: 10. initial_pressure_cgs: 1.38e-12 #erg/cm^3 stars: type: basic @@ -34,6 +33,12 @@ phases: initial_mass_Msun: 0.0 processes: + star_formation: + type: timescale_constant + params: + input_phase: gas_H2 + output_phase: stars + timescale_Gyr: 3.478 recombination: type: H_recombination params: @@ -46,9 +51,17 @@ processes: output_phase: gas_H2 agent_phase: gas_HI metallicity: 0.02 #Placeholder! one should put something like 'agent_phase2: dust' or something for better models - star_formation: - type: timescale_constant - params: - input_phase: gas_H2 - output_phase: stars - timescale_Gyr: 3.478 +# photoionization: +# type: H_photoionization +# params: +# input_phase: gas_H +# output_phase: gas_HI +# agent_process: star_formation +# efficiency: 955.29 +# photodissociation: +# type: H_photodissociation +# params: +# input_phase: gas_H2 +# output_phase: gas_H +# agent_process: star_formation +# efficiency: 380.93 \ No newline at end of file diff --git a/processes/Ionization.py b/processes/Ionization.py new file mode 100644 index 0000000..c9dea67 --- /dev/null +++ b/processes/Ionization.py @@ -0,0 +1,51 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Processes related with ionization proccesses + +Mario Romero July 2019 +""" + +import numpy as np +import astropy.units as u +from . import basic + +''' +Hydrogen Photoionization +H -> HI process +''' +class Hydrogen_Photoionization(basic.Process): + + #--------------------- + #INIT + #--------------------- + def __init__(self, model, params): + self.input = model.phases[params['output_phase']] #Neutral Hydrogen + self.output = model.phases[params['input_phase']] #Ionized Hydrogen + + self._constants() + self.tau_Gyr = self._recombination_timescale() + + def compute_derivatives(self): + flux = self.input.current_mass_Msun()/self.tau_Gyr + self.input.update_derivatives(-flux) + self.output.update_derivatives(flux) + #--------------------- + #GETTING THE TAU + #--------------------- + def _recombination_timescale(self): + #Taking Ascasibar+(in prep) formula + mean_ov_cm3_s = 4.1e-10 * (self.input.temperature())**(-0.8) #Case B recombination cross-section (Verner&Ferland 96) + + tau_s = None #Not the correct tau unit-wise + try: + tau_s = 2./(mean_ov_cm3_s * self.input.number_density()) + except ZeroDivisionError: + tau_s = np.Infinity + + return tau_s*self._s_to_Gyr #Correct units. + #--------------------- + #DEFINING USEFUL VARIABLES + #--------------------- + def _constants(self): + self._s_to_Gyr = ((1.*u.s).to(u.Gyr)).value \ No newline at end of file diff --git a/processes/MoleculeFormation.py b/processes/MoleculeFormation.py index fe050ed..3706ca5 100644 --- a/processes/MoleculeFormation.py +++ b/processes/MoleculeFormation.py @@ -32,23 +32,32 @@ def compute_derivatives(self): flux = self.input.current_mass_Msun()/self.tau_Gyr self.input.update_derivatives(-flux) self.output.update_derivatives(flux) + self.tau_Gyr = self._formation_timescale() + #--------------------- #GETTING PHYSICAL PARAMETERS #--------------------- def _dust_density(self): - total_H_density = 0.5*self.agent.number_density() + self.input.number_density() + 2.*self.output.number_density() - return self._Z * total_H_density + total_H_density_cm3 = 0.5*self.agent.number_density() + self.input.number_density() + 2.*self.output.number_density() + print(total_H_density_cm3,self.agent.number_density(),self.input.number_density(),self.output.number_density()) + return self._Z * total_H_density_cm3 def _formation_timescale(self): #Taking Ascasibar+(in prep) formula again, ignoring the effective metallicity thing - mean_ov_cm3_s = 6.e-17 *(self.input.temperature()/100.)**(0.5) #Neutral H or H2 temperature???? + mean_ov_cm3_s = 6.e-17 * ((self.input.temperature()/100.)**(0.5)) #Neutral H or H2 temperature???? + #mean_ov_cm3_Gyr = 0.189*(self.input.temperature()/100.)**(0.5) tau_s = None #Not the correct tau unit-wise try: tau_s = 0.5/(mean_ov_cm3_s * self._dust_density()) + #tau_Gyr = 0.5/(mean_ov_cm3_Gyr * self._dust_density()) except ZeroDivisionError: tau_s = np.Infinity + #tau_Gyr = np.Infinity + #print(tau_s*self._s_to_Gyr) + #print(self._s_to_Gyr) return tau_s*self._s_to_Gyr #Correct units. + #return tau_Gyr #--------------------- #DEFINING USEFUL VARIABLES #--------------------- diff --git a/processes/Recombination.py b/processes/Recombination.py index 8a5b905..895dd6a 100644 --- a/processes/Recombination.py +++ b/processes/Recombination.py @@ -30,12 +30,14 @@ def compute_derivatives(self): flux = self.input.current_mass_Msun()/self.tau_Gyr self.input.update_derivatives(-flux) self.output.update_derivatives(flux) + self.tau_Gyr = self._recombination_timescale() + #print(self.tau_Gyr) #--------------------- #GETTING THE TAU #--------------------- def _recombination_timescale(self): #Taking Ascasibar+(in prep) formula - mean_ov_cm3_s = 4.1e-10 * (self.input.temperature())**(-0.8) #Case B recombination cross-section (Verner&Ferland 96) + mean_ov_cm3_s = 4.1e-10 * ((self.input.temperature())**(-0.8)) #Case B recombination cross-section (Verner&Ferland 96) tau_s = None #Not the correct tau unit-wise try: diff --git a/processes/__init__.py b/processes/__init__.py index d6256c6..4aa6d83 100644 --- a/processes/__init__.py +++ b/processes/__init__.py @@ -1,9 +1,12 @@ from . import basic from . import Recombination from . import MoleculeFormation +from . import Ionization registry = { "timescale_constant": basic.Constant_timescale, "H_recombination": Recombination.Hydrogen_Recombination, "H2_formation": MoleculeFormation.MolecularHydrogen_Formation, + #"H_photoionization": Ionization.Hydrogen_Photoionization, + #"H_photodissociation": Ionization.Hydrogen_Photodissociation, } From 7b55d9af5d9f42c51e5238cc2ab59f9727455fbb Mon Sep 17 00:00:00 2001 From: Mario Romero Date: Tue, 30 Jul 2019 17:02:38 +0200 Subject: [PATCH 17/20] volume-averaged number densities in MoleculeFormation.py fixed --- gas_example.yml | 20 ++++++++++---------- phases/GasGenerics.py | 13 +++++++++++-- processes/MoleculeFormation.py | 17 ++++++++--------- processes/Recombination.py | 1 + 4 files changed, 30 insertions(+), 21 deletions(-) diff --git a/gas_example.yml b/gas_example.yml index eee453f..eb0670d 100644 --- a/gas_example.yml +++ b/gas_example.yml @@ -9,22 +9,22 @@ integrator: minimum_timestep_Gyr: 1.0e-4 phases: - gas_HI: + gas_HII: type: ionized_hydrogen params: - initial_mass_Msun: 99.8 + initial_mass_Msun: 99.98 initial_temperature_K: 1.0e4 initial_pressure_cgs: 1.38e-12 #erg/cm^3 - gas_H: + gas_HI: type: neutral_hydrogen params: - initial_mass_Msun: 0.1 - initial_temperature_K: 1000. + initial_mass_Msun: 0.01 + initial_temperature_K: 100. initial_pressure_cgs: 1.38e-12 #erg/cm^3 gas_H2: type: molecular_hydrogen params: - initial_mass_Msun: 0.1 + initial_mass_Msun: 0.01 initial_temperature_K: 10. initial_pressure_cgs: 1.38e-12 #erg/cm^3 stars: @@ -42,14 +42,14 @@ processes: recombination: type: H_recombination params: - input_phase: gas_HI - output_phase: gas_H + input_phase: gas_HII + output_phase: gas_HI molecule_formation: type: H2_formation params: - input_phase: gas_H + input_phase: gas_HI output_phase: gas_H2 - agent_phase: gas_HI + agent_phase: gas_HII metallicity: 0.02 #Placeholder! one should put something like 'agent_phase2: dust' or something for better models # photoionization: # type: H_photoionization diff --git a/phases/GasGenerics.py b/phases/GasGenerics.py index 5eb5fd9..8cf0182 100644 --- a/phases/GasGenerics.py +++ b/phases/GasGenerics.py @@ -26,6 +26,7 @@ def __init__(self, model, params): self.temperature_history_K = [float(self.params['initial_temperature_K'])] self.pressure_history_cgs = [float(self.params['initial_pressure_cgs'])] + #self._count = 0 self._constants() self._set() self._state() @@ -40,9 +41,14 @@ def _constants(self): def _set(self): #All internal variables will be in cgs + self._temperature = self.temperature_history_K[-1] self._pressure = self.pressure_history_cgs[-1] self._gas_mass = self.mass_history_Msun[-1]*self._solMass_to_g + #print(self.mass_history_Msun[-1],self._solMass_to_g,self.mass_history_Msun[-1]*self._solMass_to_g) + #self._count +=1 + #if self._count > 2: + # raise def current_temperature_K(self): return self.temperature_history_K[-1] @@ -119,8 +125,11 @@ def _state(self,P=-1,T=-1,M=-1,N=-1):#inputs: pressure, temperature, number of p #print("Warning: Using a parent class. Should not happen!") raise NameError("Called by parent class. Should not happen!") - #def update_derivatives(self,term): #Comes from basic - # raise NameError("Called by parent class. Should not happen!") + #--------------------- + #DEBUGGING + #--------------------- + def debug(self): + print([self._gas_mass,self._particle_mass_g]) ''' CHILDS OF GAS diff --git a/processes/MoleculeFormation.py b/processes/MoleculeFormation.py index 3706ca5..5c655fd 100644 --- a/processes/MoleculeFormation.py +++ b/processes/MoleculeFormation.py @@ -30,6 +30,7 @@ def __init__(self, model, params): def compute_derivatives(self): flux = self.input.current_mass_Msun()/self.tau_Gyr + #print(self.input.current_mass_Msun(),self.tau_Gyr,flux,'\n') self.input.update_derivatives(-flux) self.output.update_derivatives(flux) self.tau_Gyr = self._formation_timescale() @@ -38,26 +39,24 @@ def compute_derivatives(self): #GETTING PHYSICAL PARAMETERS #--------------------- def _dust_density(self): - total_H_density_cm3 = 0.5*self.agent.number_density() + self.input.number_density() + 2.*self.output.number_density() - print(total_H_density_cm3,self.agent.number_density(),self.input.number_density(),self.output.number_density()) + #total_H_density_cm3 = 0.5*self.agent.number_density() + self.input.number_density() + 2.*self.output.number_density() + total_V_cm3 = self.agent.volume() + self.input.volume() + self.output.volume() + total_H_density_cm3 = (0.5*self.agent.volume()/total_V_cm3)*self.agent.number_density() + (self.input.volume()/total_V_cm3)*self.input.number_density()+ (2.*self.output.volume()/total_V_cm3)*self.output.number_density() + #print(total_H_density_cm3) return self._Z * total_H_density_cm3 def _formation_timescale(self): #Taking Ascasibar+(in prep) formula again, ignoring the effective metallicity thing - mean_ov_cm3_s = 6.e-17 * ((self.input.temperature()/100.)**(0.5)) #Neutral H or H2 temperature???? + mean_ov_cm3_per_s = 6.e-17 * ((self.input.temperature()/100.)**(0.5)) #Neutral H or H2 temperature???? #mean_ov_cm3_Gyr = 0.189*(self.input.temperature()/100.)**(0.5) tau_s = None #Not the correct tau unit-wise try: - tau_s = 0.5/(mean_ov_cm3_s * self._dust_density()) - #tau_Gyr = 0.5/(mean_ov_cm3_Gyr * self._dust_density()) + tau_s = 0.5/(mean_ov_cm3_per_s * self._dust_density()) except ZeroDivisionError: tau_s = np.Infinity - #tau_Gyr = np.Infinity - #print(tau_s*self._s_to_Gyr) - #print(self._s_to_Gyr) return tau_s*self._s_to_Gyr #Correct units. - #return tau_Gyr + #--------------------- #DEFINING USEFUL VARIABLES #--------------------- diff --git a/processes/Recombination.py b/processes/Recombination.py index 895dd6a..fce0347 100644 --- a/processes/Recombination.py +++ b/processes/Recombination.py @@ -28,6 +28,7 @@ def __init__(self, model, params): def compute_derivatives(self): flux = self.input.current_mass_Msun()/self.tau_Gyr + #print(self.input.current_mass_Msun(),self.tau_Gyr,flux,'\n') self.input.update_derivatives(-flux) self.output.update_derivatives(flux) self.tau_Gyr = self._recombination_timescale() From 8914466cbe547d41499384eb4a9e6ecd1431211b Mon Sep 17 00:00:00 2001 From: Mario Romero Date: Thu, 1 Aug 2019 15:22:21 +0200 Subject: [PATCH 18/20] Photoionization processes added! --- SPiCE.py | 7 ++++-- gas_example.yml | 32 ++++++++++++++-------------- phases/GasGenerics.py | 1 + phases/StarGenerics.py | 26 +++++++++++++++++++++++ phases/__init__.py | 2 ++ phases/basic.py | 1 + processes/Ionization.py | 39 +++++++++++++--------------------- processes/MoleculeFormation.py | 1 - processes/__init__.py | 3 +-- processes/basic.py | 2 ++ processes/star_formation.py | 3 +-- 11 files changed, 70 insertions(+), 47 deletions(-) create mode 100644 phases/StarGenerics.py diff --git a/SPiCE.py b/SPiCE.py index fd20ff0..e38cf09 100755 --- a/SPiCE.py +++ b/SPiCE.py @@ -21,7 +21,7 @@ import processes -__version__ = "0.0.1-alpha" +__version__ = "0.0.2-alpha" def get_class(base_module, class_name): @@ -42,7 +42,7 @@ def read_config_file(self, name): def __init__(self, config_file): if config_file is None: - config_file = 'model.yml.example' + config_file = 'gas_example.yml' self.read_config_file(config_file) @@ -73,8 +73,11 @@ def run(self): timesteps_Gyr.append(phase.get_timestep_Gyr()) timestep_Gyr = np.nanmax([np.nanmin(timesteps_Gyr), self.integrator['minimum_timestep_Gyr']]) + #total_mass_sum= 0 for phase in self.phases.values(): + # total_mass_sum += phase.current_mass_Msun() phase.update_mass(timestep_Gyr) + #print(total_mass_sum) self.time_Gyr.append(current_time_Gyr + timestep_Gyr) def update_derivatives(self, term): diff --git a/gas_example.yml b/gas_example.yml index eb0670d..e9b2b4d 100644 --- a/gas_example.yml +++ b/gas_example.yml @@ -6,7 +6,7 @@ integrator: relative_accuracy: 1.0e-6 initial_time_Gyr: 0. final_time_Gyr: 13.7 - minimum_timestep_Gyr: 1.0e-4 + minimum_timestep_Gyr: 1.0e-4 #Should not be higher with the current integrator we have phases: gas_HII: @@ -28,7 +28,7 @@ phases: initial_temperature_K: 10. initial_pressure_cgs: 1.38e-12 #erg/cm^3 stars: - type: basic + type: star params: initial_mass_Msun: 0.0 @@ -51,17 +51,17 @@ processes: output_phase: gas_H2 agent_phase: gas_HII metallicity: 0.02 #Placeholder! one should put something like 'agent_phase2: dust' or something for better models -# photoionization: -# type: H_photoionization -# params: -# input_phase: gas_H -# output_phase: gas_HI -# agent_process: star_formation -# efficiency: 955.29 -# photodissociation: -# type: H_photodissociation -# params: -# input_phase: gas_H2 -# output_phase: gas_H -# agent_process: star_formation -# efficiency: 380.93 \ No newline at end of file + photoionization: + type: basic_photoionization + params: + input_phase: gas_HI + output_phase: gas_HII + agent_phase: stars + efficiency: 955.29 + photodissociation: + type: basic_photoionization + params: + input_phase: gas_H2 + output_phase: gas_HI + agent_phase: stars + efficiency: 380.93 \ No newline at end of file diff --git a/phases/GasGenerics.py b/phases/GasGenerics.py index 8cf0182..6d88b26 100644 --- a/phases/GasGenerics.py +++ b/phases/GasGenerics.py @@ -56,6 +56,7 @@ def current_pressure_cgs(self): return self.pressure_history_cgs[-1] def update_mass(self, timestep_Gyr): + #Euler integrator is becoming a bottleneck, it should be replaced with a more efficient one. #Actually, this should be called 'update_all', but is called this way to keep compatibility self.mass_history_Msun.append(self.current_mass_Msun() + self.dm_dt_Msun_Gyr*timestep_Gyr) diff --git a/phases/StarGenerics.py b/phases/StarGenerics.py new file mode 100644 index 0000000..cf25106 --- /dev/null +++ b/phases/StarGenerics.py @@ -0,0 +1,26 @@ +''' +Star classes + +Mario Romero May 2019 +''' + +import numpy as np +import astropy.units as u +import astropy.constants as cte +from . import basic + +''' +PLACEHOLDER CLASS +Defined here to make photo-ionization and photo-dissociation processes work +''' +class Star(basic.Phase): + + #--------------------- + #DEFAULT SETTINGS + #--------------------- + def __init__(self, model, params): + #What everyone does + self.model = model + self.params = {**self.default_settings(), **params} + self.mass_history_Msun = [float(self.params['initial_mass_Msun'])] + self.SFR_history_Msun_per_Gyr = [] #It updates when a star formation process is called \ No newline at end of file diff --git a/phases/__init__.py b/phases/__init__.py index 768a69d..ad2f126 100644 --- a/phases/__init__.py +++ b/phases/__init__.py @@ -1,5 +1,6 @@ from . import basic # from . import GasGenerics +from . import StarGenerics from . import GasSpecies from . import example @@ -7,6 +8,7 @@ registry = { "basic": basic.Phase, "multiphase": basic.MultiphaseMedium, + "star": StarGenerics.Star, #PLACEHOLDER! "ionized_hydrogen": GasSpecies.Ionized_Hydrogen, "neutral_hydrogen": GasSpecies.Neutral_Hydrogen, "molecular_hydrogen": GasSpecies.Molecular_Hydrogen, diff --git a/phases/basic.py b/phases/basic.py index 2689fbd..c65e34d 100644 --- a/phases/basic.py +++ b/phases/basic.py @@ -44,6 +44,7 @@ def get_timestep_Gyr(self): return np.Infinity def update_mass(self, timestep_Gyr): + #Euler integrator is becoming a bottleneck, it should be replaced with a more efficient one. self.mass_history_Msun.append(self.current_mass_Msun() + self.dm_dt_Msun_Gyr*timestep_Gyr) diff --git a/processes/Ionization.py b/processes/Ionization.py index c9dea67..613f851 100644 --- a/processes/Ionization.py +++ b/processes/Ionization.py @@ -12,40 +12,31 @@ ''' Hydrogen Photoionization -H -> HI process +HI -> HII process ''' -class Hydrogen_Photoionization(basic.Process): +class Basic_Photoionization(basic.Process): #--------------------- #INIT #--------------------- def __init__(self, model, params): - self.input = model.phases[params['output_phase']] #Neutral Hydrogen - self.output = model.phases[params['input_phase']] #Ionized Hydrogen - - self._constants() - self.tau_Gyr = self._recombination_timescale() + self.input = model.phases[params['input_phase']] + self.output = model.phases[params['output_phase']] + self.agent = model.phases[params['agent_phase']] #Stars + self._eff = float(params['efficiency']) #=Mass of photoionized gas per Mass of stars def compute_derivatives(self): - flux = self.input.current_mass_Msun()/self.tau_Gyr + #flux = self.input.current_mass_Msun()/self.tau_Gyr + try: + flux = self._eff * self.agent.SFR_history_Msun_per_Gyr[-1] + except IndexError: + flux = 0.0 + #print(self.agent.SFR_history_Msun_per_Gyr[0]) self.input.update_derivatives(-flux) self.output.update_derivatives(flux) - #--------------------- - #GETTING THE TAU - #--------------------- - def _recombination_timescale(self): - #Taking Ascasibar+(in prep) formula - mean_ov_cm3_s = 4.1e-10 * (self.input.temperature())**(-0.8) #Case B recombination cross-section (Verner&Ferland 96) - tau_s = None #Not the correct tau unit-wise try: - tau_s = 2./(mean_ov_cm3_s * self.input.number_density()) + self.tau_Gyr = self.input.current_mass_Msun() / flux except ZeroDivisionError: - tau_s = np.Infinity - - return tau_s*self._s_to_Gyr #Correct units. - #--------------------- - #DEFINING USEFUL VARIABLES - #--------------------- - def _constants(self): - self._s_to_Gyr = ((1.*u.s).to(u.Gyr)).value \ No newline at end of file + self.tau_Gyr = np.Infinity + #print(self.model.input) \ No newline at end of file diff --git a/processes/MoleculeFormation.py b/processes/MoleculeFormation.py index 5c655fd..5db9272 100644 --- a/processes/MoleculeFormation.py +++ b/processes/MoleculeFormation.py @@ -30,7 +30,6 @@ def __init__(self, model, params): def compute_derivatives(self): flux = self.input.current_mass_Msun()/self.tau_Gyr - #print(self.input.current_mass_Msun(),self.tau_Gyr,flux,'\n') self.input.update_derivatives(-flux) self.output.update_derivatives(flux) self.tau_Gyr = self._formation_timescale() diff --git a/processes/__init__.py b/processes/__init__.py index 4aa6d83..a34cf86 100644 --- a/processes/__init__.py +++ b/processes/__init__.py @@ -7,6 +7,5 @@ "timescale_constant": basic.Constant_timescale, "H_recombination": Recombination.Hydrogen_Recombination, "H2_formation": MoleculeFormation.MolecularHydrogen_Formation, - #"H_photoionization": Ionization.Hydrogen_Photoionization, - #"H_photodissociation": Ionization.Hydrogen_Photodissociation, + "basic_photoionization": Ionization.Basic_Photoionization, } diff --git a/processes/basic.py b/processes/basic.py index cd87e46..c62119b 100644 --- a/processes/basic.py +++ b/processes/basic.py @@ -26,3 +26,5 @@ def compute_derivatives(self): flux = self.input.current_mass_Msun()/self.tau_Gyr self.input.update_derivatives(-flux) self.output.update_derivatives(flux) + #PLACEHOLDER! This will fail if two processes give an SFR as main result + self.output.SFR_history_Msun_per_Gyr.append(flux) diff --git a/processes/star_formation.py b/processes/star_formation.py index 2ba27fd..b9708af 100644 --- a/processes/star_formation.py +++ b/processes/star_formation.py @@ -1,5 +1,3 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ Created on Wed Mar 27 10:30:28 2019 @@ -21,3 +19,4 @@ def compute_derivatives(self): SFR = self.gas.mass()/self.tau self.gas.update_derivatives(-SFR) self.stars.update_derivatives(SFR) + From a30ce8a3ffea735ec37cbaf706b65347c4c45ecf Mon Sep 17 00:00:00 2001 From: Mario Romero Date: Thu, 1 Aug 2019 16:07:25 +0200 Subject: [PATCH 19/20] QoL changes: Code will crash if finds x/0, overflow, underflow and NaN --- SPiCE.py | 2 ++ gas_example.yml | 6 +++--- phases/GasGenerics.py | 4 +++- phases/basic.py | 3 ++- processes/Ionization.py | 2 +- 5 files changed, 11 insertions(+), 6 deletions(-) diff --git a/SPiCE.py b/SPiCE.py index e38cf09..0195231 100755 --- a/SPiCE.py +++ b/SPiCE.py @@ -41,6 +41,8 @@ def read_config_file(self, name): self.context = settings.validate(user_settings) def __init__(self, config_file): + #The code crash with 'FloatingPointError' for the cases: "x/0", overflow, underflow and nan ("0/0" or "inf/inf") + np.seterr(divide='raise',over='raise',under='raise',invalid='raise') if config_file is None: config_file = 'gas_example.yml' diff --git a/gas_example.yml b/gas_example.yml index e9b2b4d..7e6d7f9 100644 --- a/gas_example.yml +++ b/gas_example.yml @@ -12,19 +12,19 @@ phases: gas_HII: type: ionized_hydrogen params: - initial_mass_Msun: 99.98 + initial_mass_Msun: 100. initial_temperature_K: 1.0e4 initial_pressure_cgs: 1.38e-12 #erg/cm^3 gas_HI: type: neutral_hydrogen params: - initial_mass_Msun: 0.01 + initial_mass_Msun: 0.0 initial_temperature_K: 100. initial_pressure_cgs: 1.38e-12 #erg/cm^3 gas_H2: type: molecular_hydrogen params: - initial_mass_Msun: 0.01 + initial_mass_Msun: 0.0 initial_temperature_K: 10. initial_pressure_cgs: 1.38e-12 #erg/cm^3 stars: diff --git a/phases/GasGenerics.py b/phases/GasGenerics.py index 6d88b26..ad1ed4e 100644 --- a/phases/GasGenerics.py +++ b/phases/GasGenerics.py @@ -58,6 +58,7 @@ def current_pressure_cgs(self): def update_mass(self, timestep_Gyr): #Euler integrator is becoming a bottleneck, it should be replaced with a more efficient one. #Actually, this should be called 'update_all', but is called this way to keep compatibility + #print(self.current_mass_Msun()) self.mass_history_Msun.append(self.current_mass_Msun() + self.dm_dt_Msun_Gyr*timestep_Gyr) #If you don't want to store the pressure history because, for example, it does not change over time @@ -84,12 +85,13 @@ def _compute_variables_indepedent_of_state(self): self._number_density = self._number_particles / self._volume self._thermal_energy_density = self._thermal_energy / self._volume self._gamma = 1.+(self._pressure/self._thermal_energy_density) - except ZeroDivisionError: + except (FloatingPointError,ZeroDivisionError): self._number_density = 0.0 self._thermal_energy_density = 0.0 self._gamma = np.Infinity self._mass_density = self._number_density * self._particle_mass_g self._chemical_potential = self._kB*self._temperature*np.log(self._fugacity) + #print(self._thermal_energy_density , self._thermal_energy , self._volume) #--------------------- #OUTPUTS diff --git a/phases/basic.py b/phases/basic.py index c65e34d..225f96c 100644 --- a/phases/basic.py +++ b/phases/basic.py @@ -40,7 +40,8 @@ def get_timestep_Gyr(self): try: return np.abs((self.model.integrator['relative_accuracy'] * self.current_mass_Msun() / self.dm_dt_Msun_Gyr)) - except ZeroDivisionError: + except (FloatingPointError,ZeroDivisionError): + #How should the previous expression be solved if gives 0/0? return np.Infinity def update_mass(self, timestep_Gyr): diff --git a/processes/Ionization.py b/processes/Ionization.py index 613f851..fa76177 100644 --- a/processes/Ionization.py +++ b/processes/Ionization.py @@ -37,6 +37,6 @@ def compute_derivatives(self): try: self.tau_Gyr = self.input.current_mass_Msun() / flux - except ZeroDivisionError: + except (FloatingPointError,ZeroDivisionError): self.tau_Gyr = np.Infinity #print(self.model.input) \ No newline at end of file From cee383f76be2498f39993eb1137dab930d9da78f Mon Sep 17 00:00:00 2001 From: Mario Romero Date: Fri, 2 Aug 2019 14:13:23 +0200 Subject: [PATCH 20/20] Gas pressure now vary with time using Elmegreen89 (Leroy+08) Formula --- gas_example.yml | 2 +- phases/GasGenerics.py | 138 +++++++++++++++++++++++++++++++++++++----- 2 files changed, 123 insertions(+), 17 deletions(-) diff --git a/gas_example.yml b/gas_example.yml index 7e6d7f9..736c839 100644 --- a/gas_example.yml +++ b/gas_example.yml @@ -5,7 +5,7 @@ dtd_sn: rlp integrator: relative_accuracy: 1.0e-6 initial_time_Gyr: 0. - final_time_Gyr: 13.7 + final_time_Gyr: 13.7 minimum_timestep_Gyr: 1.0e-4 #Should not be higher with the current integrator we have phases: diff --git a/phases/GasGenerics.py b/phases/GasGenerics.py index ad1ed4e..32bd441 100644 --- a/phases/GasGenerics.py +++ b/phases/GasGenerics.py @@ -26,11 +26,19 @@ def __init__(self, model, params): self.temperature_history_K = [float(self.params['initial_temperature_K'])] self.pressure_history_cgs = [float(self.params['initial_pressure_cgs'])] - #self._count = 0 self._constants() self._set() self._state() + def reset_timestep(self): + self.dm_dt_Msun_Gyr = 0. + self.dP_dt_cgs_Gyr = 0. + + def current_temperature_K(self): + return self.temperature_history_K[-1] + def current_pressure_cgs(self): + return self.pressure_history_cgs[-1] + def _constants(self): #Useful conversions and constants (this will save us LOTS of parentesis, trust me!) self._solMass_to_g = ((1*u.solMass).to(u.g)).value @@ -38,6 +46,7 @@ def _constants(self): self._hP = ((cte.h).to(u.erg*u.s)).value #Planck constant in erg*s self._hbar = self._hP/(2.*np.pi) self._c = ((cte.c).to(u.cm/u.s)).value #Speed of light + self._G = ((cte.G).to(u.cm*u.cm*u.cm/(u.g*u.s*u.s))).value #Gravitational constant def _set(self): #All internal variables will be in cgs @@ -45,26 +54,31 @@ def _set(self): self._temperature = self.temperature_history_K[-1] self._pressure = self.pressure_history_cgs[-1] self._gas_mass = self.mass_history_Msun[-1]*self._solMass_to_g - #print(self.mass_history_Msun[-1],self._solMass_to_g,self.mass_history_Msun[-1]*self._solMass_to_g) - #self._count +=1 - #if self._count > 2: - # raise - - def current_temperature_K(self): - return self.temperature_history_K[-1] - def current_pressure_cgs(self): - return self.pressure_history_cgs[-1] def update_mass(self, timestep_Gyr): #Euler integrator is becoming a bottleneck, it should be replaced with a more efficient one. #Actually, this should be called 'update_all', but is called this way to keep compatibility - #print(self.current_mass_Msun()) + #------------------------------ + + #First, I deal with Elmegreen89 Formula to get the pressure (or the box area if first iteration) + new_pressure = self.current_pressure_cgs() + try: + self._box_area_cm2 #It will go to the except block for the first iteration + #self._compute_pressure_derivative() + new_pressure = self._compute_pressure_directly() + except AttributeError: + self._set_box_area() + #Update mass self.mass_history_Msun.append(self.current_mass_Msun() + self.dm_dt_Msun_Gyr*timestep_Gyr) + #Update pressure + #self.pressure_history_cgs.append(self.current_pressure_cgs() + self.dP_dt_cgs_Gyr*timestep_Gyr) + self.pressure_history_cgs.append(new_pressure) + #print(self.current_pressure_cgs()/self._kB,self.current_mass_Msun(),type(self)) + #If you don't want to store the pressure history because, for example, it does not change over time # you only have to comment the relevant line of these two. self.temperature_history_K.append(self.current_temperature_K()) - self.pressure_history_cgs.append(self.current_pressure_cgs()) self._set() self._state() @@ -85,14 +99,100 @@ def _compute_variables_indepedent_of_state(self): self._number_density = self._number_particles / self._volume self._thermal_energy_density = self._thermal_energy / self._volume self._gamma = 1.+(self._pressure/self._thermal_energy_density) + self._chemical_potential = self._kB*self._temperature*np.log(self._fugacity) except (FloatingPointError,ZeroDivisionError): self._number_density = 0.0 self._thermal_energy_density = 0.0 self._gamma = np.Infinity + self._chemical_potential = -np.Infinity self._mass_density = self._number_density * self._particle_mass_g - self._chemical_potential = self._kB*self._temperature*np.log(self._fugacity) #print(self._thermal_energy_density , self._thermal_energy , self._volume) - + + def _set_box_area(self): + #Here we use the Elmegreen89 (approximate) formula for pressure: + #P_gas = (pi/2) * G * (S_gas + S_Total) + #Where S is the superficial density = Mass/Area. We don't have the area, so we compute it here as a 'constant' + #Here we get the Area + + #Get total mass, and total GAS mass + total_mass = 0.0 + total_gas_mass = 0.0 + total_gas_pressure = 0.0 + for phase in self.model.phases.values(): + total_mass += phase.mass_history_Msun[0] #We need only the first value, and it may be updated, so using 'current_mass' is a bad idea + if( issubclass(type(phase),Gas) ): + total_gas_mass += phase.mass_history_Msun[0] + total_gas_pressure += phase.pressure_history_cgs[0] + #Convert to cgs + total_mass *= self._solMass_to_g + total_gas_mass *= self._solMass_to_g + #Get the result + self._box_area_cm2 = np.sqrt( np.pi*self._G*total_gas_mass*total_mass / (2.*total_gas_pressure) ) + + def _compute_pressure_directly(self): + #Here we use the Elmegreen89 (approximate) formula for pressure: + #P_i = (pi/2) * G * (S_i + S_Total) + #Where S is the superficial density = Mass/Area. We don't have the area, so we compute it here as a 'constant' + #Here we compute P_i using the formula. + + #Here we need the total mass again. + #Total mass of a previous phase may be updated before in the main routine. We have to cover this! + total_mass = 0.0 + thisPhase_historylength = len(self.mass_history_Msun) + for phase in self.model.phases.values(): + currentPhase_historylength = len(phase.mass_history_Msun) + if( thisPhase_historylength == currentPhase_historylength ): + #currentPhase has not been updated yet, use current value + total_mass += phase.mass_history_Msun[-1] + elif( thisPhase_historylength == currentPhase_historylength-1 ): + #currentPhase has been updated in this timestep before this call + total_mass += phase.mass_history_Msun[-2] + else: + print("This should not happen!") + raise(-1) + #Make the conversion Msun->g + total_mass = total_mass * self._solMass_to_g + current_mass = self.current_mass_Msun() * self._solMass_to_g + #Get the new pressure + return current_mass*total_mass*( np.pi*self._G / (2.*self._box_area_cm2*self._box_area_cm2) ) + + + ''' + def _compute_pressure_derivative(self): + #Here we use the Elmegreen89 (approximate) formula for pressure: + #P_i = (pi/2) * G * (S_i + S_Total) + #Where S is the superficial density = Mass/Area. We don't have the area, so we compute it here as a 'constant' + #Here we compute the derivative of P_i + + #print( (np.sqrt(self._box_area_cm2/np.pi)*u.cm).to(u.pc) ) #Need to check if we have the correct numbers! + + #Here we need the sum of all derivatives of dm_dt, and again the total mass. + #Total mass of a previous phase may be updated before in the main routine. We have to cover this! + total_dm_dt_Msun_Gyr = 0.0 + total_mass = 0.0 + thisPhase_historylength = len(self.mass_history_Msun) + for phase in self.model.phases.values(): + currentPhase_historylength = len(phase.mass_history_Msun) + if( thisPhase_historylength == currentPhase_historylength ): + #currentPhase has not been updated yet, use current value + total_mass += phase.mass_history_Msun[-1] + elif( thisPhase_historylength == currentPhase_historylength-1 ): + #currentPhase has been updated in this timestep before this call + total_mass += phase.mass_history_Msun[-2] + else: + print("This should not happen!") + raise(-1) + #Sum derivatives + total_dm_dt_Msun_Gyr += phase.dm_dt_Msun_Gyr + #Make the conversion Msun->g + total_dm_dt_g_Gyr = total_dm_dt_Msun_Gyr * self._solMass_to_g + curr_dm_dt_g_Gyr = self.dm_dt_Msun_Gyr * self._solMass_to_g + total_mass = total_mass * self._solMass_to_g + current_mass = self.current_mass_Msun() * self._solMass_to_g + #Get the pressure derivative + self.dP_dt_cgs_Gyr = curr_dm_dt_g_Gyr*total_mass + current_mass*total_dm_dt_g_Gyr + self.dP_dt_cgs_Gyr *= np.pi*self._G / (2.*self._box_area_cm2*self._box_area_cm2) + ''' #--------------------- #OUTPUTS #--------------------- @@ -148,7 +248,10 @@ def _state(self):#inputs: pressure, temperature, number of particles, thermal en self._number_particles = (self._gas_mass)/(self._particle_mass_g) self._thermal_energy = 1.5*self._number_particles * self._kB * self._temperature - self._volume = (2.*self._thermal_energy) / (3.*self._pressure) + try: + self._volume = (2.*self._thermal_energy) / (3.*self._pressure) + except (FloatingPointError,ZeroDivisionError): + self._volume = 0.0 #Fugacity = exp(chemical_potential / kT) . It is another way to describe the chemical potential self._fugacity = self._pressure * (2.*np.pi*self._particle_mass_g / (self._hP*self._hP) )**(-1.5) * (self._kB*self._temperature)**(-2.5) @@ -239,7 +342,10 @@ def _state(self):#inputs: pressure, temperature, number of particles, mass Zvib = self._Z_vib() #Now we compute the remaining variables - self._volume = NkT / self._pressure + try: + self._volume = NkT / self._pressure + except (FloatingPointError,ZeroDivisionError): + self._volume = 0.0 #Fugacity = exp(chemical_potential / kT) . It is another way to describe the chemical potential self._fugacity = self._pressure * (2.*np.pi*self._particle_mass_g / (self._hP*self._hP) )**(-1.5) * (self._kB*self._temperature)**(-2.5) / (Zrot[0]*Zvib[0]) #Thermal energy