Skip to content
Snippets Groups Projects
Commit cdf5e10f authored by Loic Hausammann's avatar Loic Hausammann
Browse files

Add ParameterFile

parent 8d0629e4
No related branches found
No related tags found
1 merge request!9New implementation
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.989e43 # Grams
UnitLength_in_cgs: 3.085e21 # Centimeters
UnitVelocity_in_cgs: 1e5 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
# Cooling with Grackle 2.0
GrackleCooling:
CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table
WithUVbackground: 1 # Enable or not the UV background
Redshift: 0 # Redshift to use (-1 means time based redshift)
WithMetalCooling: 0 # Enable or not the metal cooling
ProvideVolumetricHeatingRates: 0 # User provide volumetric heating rates
ProvideSpecificHeatingRates: 0 # User provide specific heating rates
SelfShieldingMethod: 0 # Grackle (<= 3) or Gear self shielding method
OutputMode: 0 # Write in output corresponding primordial chemistry mode
MaxSteps: 1000
ConvergenceLimit: 1e-2
GearChemistry:
InitialMetallicity: 0.01295
\ No newline at end of file
......@@ -18,7 +18,7 @@
# ########################################################################
from pyswiftsim.libcooling import coolingRate
from pyswiftsim import downloadCoolingTable
import pyswiftsim
import numpy as np
import matplotlib.pyplot as plt
......@@ -28,15 +28,19 @@ from time import strftime, gmtime
plt.rc('text', usetex=True)
downloadCoolingTable()
pyswiftsim.downloadCoolingTable()
# PARAMETERS
# swift params filename
filename = "cooling.yml"
# use cosmology?
with_cosmo = False
# include metal cooling?
with_metals = 1
# particle data
# in atom/cc
part_density = np.array([0.1])
......@@ -123,31 +127,31 @@ def plot_solution(energy, data, us):
plt.show()
def getParser(filename):
with open(filename, "r") as stream:
stream = yaml.load(stream)
return stream
if __name__ == "__main__":
us = getParser(filename)["InternalUnitSystem"]
overwrite = {
"GrackleCooling": {
"WithMetalCooling": with_metals,
}
}
with pyswiftsim.ParameterFile(overwrite) as filename:
us = pyswiftsim.parseYamlFile(filename)["InternalUnitSystem"]
d = generate_initial_condition(us)
d = generate_initial_condition(us)
t = d["time"]
dt = t[1] - t[0]
t = d["time"]
dt = t[1] - t[0]
# du / dt
print("Computing cooling...")
# du / dt
print("Computing cooling...")
N = len(t)
energy = np.zeros(t.shape)
energy[0] = d["energy"]
rho = np.array(d["density"], dtype=np.float32)
for i in range(N-1):
ene = np.array([energy[i]], dtype=np.float32)
rate = coolingRate(filename, rho, ene, with_cosmo, dt)
energy[i+1] = energy[i] + rate * dt
N = len(t)
energy = np.zeros(t.shape)
energy[0] = d["energy"]
rho = np.array(d["density"], dtype=np.float32)
for i in range(N-1):
ene = np.array([energy[i]], dtype=np.float32)
rate = coolingRate(filename, rho, ene, with_cosmo, dt)
energy[i+1] = energy[i] + rate * dt
plot_solution(energy, d, us)
plot_solution(energy, d, us)
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.989e43 # Grams
UnitLength_in_cgs: 3.085e21 # Centimeters
UnitVelocity_in_cgs: 1e5 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
minimal_temperature: 0 # (Optional) Minimal temperature (in internal units) allowed for the gas particles. Value is ignored if set to 0.
# Cosmological parameters
Cosmology:
h: 0.6777 # Reduced Hubble constant
a_begin: 0.0078125 # Initial scale-factor of the simulation
a_end: 1.0 # Final scale factor of the simulation
Omega_m: 0.307 # Matter density parameter
Omega_lambda: 0.693 # Dark-energy density parameter
Omega_b: 0.0455 # Baryon density parameter
Omega_r: 0. # (Optional) Radiation density parameter
w_0: -1.0 # (Optional) Dark-energy equation-of-state parameter at z=0.
w_a: 0. # (Optional) Dark-energy equation-of-state time evolution parameter.
# Cooling with Grackle 2.0
GrackleCooling:
CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table
WithUVbackground: 1 # Enable or not the UV background
Redshift: 0 # Redshift to use (-1 means time based redshift)
WithMetalCooling: 1 # Enable or not the metal cooling
ProvideVolumetricHeatingRates: 0 # User provide volumetric heating rates
ProvideSpecificHeatingRates: 0 # User provide specific heating rates
SelfShieldingMethod: 0 # Grackle (<= 3) or Gear self shielding method
ConvergenceLimit: 1e-2
MaxSteps: 20000
Omega: 0.8
GearChemistry:
InitialMetallicity: 0.01295
\ No newline at end of file
......@@ -18,19 +18,18 @@
# ########################################################################
from pyswiftsim.libcooling import coolingRate
from pyswiftsim import downloadCoolingTable
import pyswiftsim
from copy import deepcopy
import numpy as np
import matplotlib.pyplot as plt
from astropy import units, constants
import yaml
from time import strftime, gmtime
import sys
plt.rc('text', usetex=True)
downloadCoolingTable()
pyswiftsim.downloadCoolingTable()
# PARAMETERS
......@@ -39,8 +38,8 @@ colormap = "RdBu"
# cosmology
with_cosmo = False
# swift params filename
filename = "cooling.yml"
# include metals?
with_metals = 1
# hydrogen mass fraction
h_mass_frac = 0.76
......@@ -114,12 +113,6 @@ def generate_initial_condition(us):
return d
def getParser(filename):
with open(filename, "r") as stream:
stream = yaml.load(stream)
return stream
def plot_solution(rate, data, us):
print("Plotting solution")
energy = data["energy"]
......@@ -195,17 +188,24 @@ def plot_solution(rate, data, us):
if __name__ == "__main__":
parser = getParser(filename)
us = parser["InternalUnitSystem"]
overwrite = {
"GrackleCooling": {
"WithMetalCooling": with_metals,
}
}
with pyswiftsim.ParameterFile(overwrite) as filename:
parser = pyswiftsim.parseYamlFile(filename)
us = parser["InternalUnitSystem"]
d = generate_initial_condition(us)
d = generate_initial_condition(us)
# du / dt
print("Computing cooling...")
# du / dt
print("Computing cooling...")
rate = coolingRate(filename,
d["density"].astype(np.float32),
d["energy"].astype(np.float32),
with_cosmo, d["time_step"])
rate = coolingRate(filename,
d["density"].astype(np.float32),
d["energy"].astype(np.float32),
with_cosmo, d["time_step"])
plot_solution(rate, d, us)
plot_solution(rate, d, us)
......@@ -18,6 +18,8 @@
import os
import urllib.request
from tempfile import NamedTemporaryFile
import yaml
def downloadCoolingTable():
......@@ -28,3 +30,132 @@ def downloadCoolingTable():
if not os.path.isfile(filename):
# Download the file from `url` and save it locally under `file_name`:
urllib.request.urlretrieve(url, filename)
def parseYamlFile(filename):
with open(filename, "r") as stream:
stream = yaml.load(stream)
return stream
class ParameterFile:
default_parameters = {
"InternalUnitSystem": {
"UnitMass_in_cgs": 1.989e43, # 1e10 Msun
"UnitLength_in_cgs": 3.085e21, # Kpc
"UnitVelocity_in_cgs": 1e5, # km / s
"UnitCurrent_in_cgs": 1, # Amperes
"UnitTemp_in_cgs": 1, # Kelvin
},
"SPH": {
"resolution_eta": 1.2348,
"CFL_condition": 0.1,
"minimal_temperature": 0
},
"Cosmology": {
"h": 0.6777,
"a_begin": 0.0078125,
"a_end": 1.0,
"Omega_m": 0.307,
"Omega_lambda": 0.693,
"Omega_b": 0.0455,
"Omega_r": 0.,
"w_0": -1.0,
"w_a": 0.
},
"GrackleCooling": {
"CloudyTable": "CloudyData_UVB=HM2012.h5",
"WithUVbackground": 1,
"Redshift": 0,
"WithMetalCooling": 1,
"ProvideVolumetricHeatingRates": 0,
"ProvideSpecificHeatingRates": 0,
"SelfShieldingMethod": 0,
"ConvergenceLimit": 1e-2,
"MaxSteps": 20000,
"Omega": 0.8
},
"GearChemistry": {
"InitialMetallicity": 0.01295,
}
}
def _write(self, f, txt):
txt = txt.encode()
f.write(txt)
def _writeParameters(self, f, d, overwrite):
_format = " {}: {}\n"
for sec in d.keys():
p = d[sec]
# write section
txt = "{}:\n".format(sec)
self._write(f, txt)
for param in p.keys():
value = p[param]
# test if overwritten
test = overwrite is not None and sec in overwrite
test = test and param in overwrite[sec]
if test:
value = overwrite[sec][param]
# write the parameter
txt = _format.format(param, value)
self._write(f, txt)
# check if in the overwrite dictionary
if overwrite is None or sec not in overwrite:
continue
# ensure that all the elements in overwrite are written
for param in overwrite[sec].keys():
# check if already done
if param in p:
continue
# write the parameter
txt = _format.format(param, overwrite[sec][param])
self._write(f, txt)
# add a space between sections
print()
if overwrite is None:
return
# now do all the sections that do not exist in default
for sec in overwrite.keys():
# check if already done
if sec in d:
continue
# write section
txt = "{}:\n".format(sec)
self._write(f, txt)
for param in overwrite[sec].keys():
# write parameter
txt = _format.format(param, overwrite[sec][param])
self._write(f, txt)
# add space between sections
print()
def __init__(self, overwrite=None):
self.overwrite = overwrite
self.d = self.default_parameters
def __enter__(self):
f = NamedTemporaryFile(delete=False)
self.filename = f.name
self._writeParameters(f, self.default_parameters, self.overwrite)
f.close()
return self.filename
def __exit__(self, *args):
os.remove(self.filename)
#!/usr/bin/env python3
# ########################################################################
# This file is part of PYSWIFT.
# Copyright (c) 2018 Loic Hausammann (loic.hausammann@epfl.ch)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# ########################################################################
from pyswiftsim import ParameterFile, parseYamlFile
import os
# defines a few parameter that will be overwritten
overwrite = {
"InternalUnitSystem": {
"UnitMass_in_cgs": 1,
},
"test": {
"test": -1
}
}
print("Testing the default file")
# Generate and read the default file
with ParameterFile() as params:
default = parseYamlFile(params)
# check if the temporary file is removed
if os.path.isfile(params):
raise Exception("File {} not destroyed".format(params))
print("Testing to overwrite the default file")
# Now do the same but with the overwrite feature
with ParameterFile(overwrite) as params:
over = parseYamlFile(params)
# check if the temporary file is removed
if os.path.isfile(params):
raise Exception("File {} not destroyed".format(params))
print("Comparing results")
# Now check that the correct values are found
for sec in over.keys():
params = over[sec]
for p in params.keys():
# First case not in default parameters
if sec not in default:
assert overwrite[sec][p] == params[p]
continue
# Second case in default parameters and overwritted
if sec in overwrite and p in overwrite[sec]:
assert overwrite[sec][p] == params[p]
continue
# Third case in default
assert params[p] == default[sec][p]
print("Checking if all the elements exist")
# Finally ensure that all the elements are present
for sec in default.keys():
params = default[sec]
for p in params.keys():
if sec not in over:
raise Exception("Section {} not found".format(sec))
if p not in over[sec]:
raise Exception("Parameter {}:{} not found".format(sec, p))
......@@ -18,7 +18,7 @@
# ########################################################################
from pyswiftsim.libcooling import coolingRate
from pyswiftsim import downloadCoolingTable
import pyswiftsim
import numpy as np
from astropy import units, constants
import yaml
......@@ -29,58 +29,63 @@ filename = "test_cooling.yml"
gamma = 5. / 3.
# download cooling table
downloadCoolingTable()
# Read parameter file
f = open(filename, "r")
data = yaml.load(f)
f.close()
# Get the units
swift_units = data["InternalUnitSystem"]
# length unit
u_l = float(swift_units["UnitLength_in_cgs"])
# velocity unit
u_v = float(swift_units["UnitVelocity_in_cgs"])
# mass unit
u_m = float(swift_units["UnitMass_in_cgs"])
# time unit
u_t = u_l / u_v
# energy unit
u_e = u_m * u_l**2 / u_t**2
# boltzman constant
k_B = constants.k_B.to("erg / K") * units.K / (units.erg * u_e)
k_B = float(k_B)
# proton mass
m_p = constants.m_p.to("g") / (units.g * u_m)
m_p = float(m_p)
def internalEnergy(T):
u = (k_B / m_p) * T
u /= (gamma - 1.)
return u
# generates data
# time step in Myr
dt = 1e-3
dt = float(units.Myr.to('s')) / u_t
# density in atom / cm3
density = np.array([1.], dtype=np.float32)
density *= m_p * u_l**3
# Temperature in K
T = np.array([1e4], dtype=np.float32)
u = internalEnergy(T)
# compute the cooling rate
with_cosmo = False
rate = coolingRate(filename, density, u, with_cosmo, dt)
# go back to cgs
rate /= u_e / u_t
print("Cooling done: {} erg/s".format(rate))
pyswiftsim.downloadCoolingTable()
def computeCooling(filename):
# Read parameter file
f = open(filename, "r")
data = yaml.load(f)
f.close()
# Get the units
swift_units = data["InternalUnitSystem"]
# length unit
u_l = float(swift_units["UnitLength_in_cgs"])
# velocity unit
u_v = float(swift_units["UnitVelocity_in_cgs"])
# mass unit
u_m = float(swift_units["UnitMass_in_cgs"])
# time unit
u_t = u_l / u_v
# energy unit
u_e = u_m * u_l**2 / u_t**2
# boltzman constant
k_B = constants.k_B.to("erg / K") * units.K / (units.erg * u_e)
k_B = float(k_B)
# proton mass
m_p = constants.m_p.to("g") / (units.g * u_m)
m_p = float(m_p)
def internalEnergy(T):
u = (k_B / m_p) * T
u /= (gamma - 1.)
return u
# generates data
# time step in Myr
dt = 1e-3
dt = float(units.Myr.to('s')) / u_t
# density in atom / cm3
density = np.array([1.], dtype=np.float32)
density *= m_p * u_l**3
# Temperature in K
T = np.array([1e4], dtype=np.float32)
u = internalEnergy(T)
# compute the cooling rate
with_cosmo = False
rate = coolingRate(filename, density, u, with_cosmo, dt)
# go back to cgs
rate /= u_e / u_t
print("Cooling done: {} erg/s".format(rate))
with pyswiftsim.ParameterFile() as filename:
computeCooling(filename)
File moved
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment