From cdf5e10f915eb5997d16afae0dbb033d9b9659be Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Mon, 27 May 2019 11:40:08 +0200
Subject: [PATCH] Add ParameterFile
---
examples/cooling_box/cooling.yml | 28 ------
examples/cooling_box/plot_cooling.py | 50 +++++-----
examples/cooling_rate/cooling.yml | 41 --------
examples/cooling_rate/plot_cooling.py | 42 ++++-----
pyswiftsim/__init__.py | 131 ++++++++++++++++++++++++++
test/test_cooling.py | 86 -----------------
tests/generate_parameter_file.py | 85 +++++++++++++++++
tests/test_cooling.py | 91 ++++++++++++++++++
{test => tests}/test_cooling.yml | 0
9 files changed, 355 insertions(+), 199 deletions(-)
delete mode 100644 examples/cooling_box/cooling.yml
delete mode 100644 examples/cooling_rate/cooling.yml
delete mode 100644 test/test_cooling.py
create mode 100644 tests/generate_parameter_file.py
create mode 100644 tests/test_cooling.py
rename {test => tests}/test_cooling.yml (100%)
diff --git a/examples/cooling_box/cooling.yml b/examples/cooling_box/cooling.yml
deleted file mode 100644
index 66fa4b6..0000000
--- a/examples/cooling_box/cooling.yml
+++ /dev/null
@@ -1,28 +0,0 @@
-# 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
diff --git a/examples/cooling_box/plot_cooling.py b/examples/cooling_box/plot_cooling.py
index 6d8485e..18bbb26 100644
--- a/examples/cooling_box/plot_cooling.py
+++ b/examples/cooling_box/plot_cooling.py
@@ -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)
diff --git a/examples/cooling_rate/cooling.yml b/examples/cooling_rate/cooling.yml
deleted file mode 100644
index 8fa0e72..0000000
--- a/examples/cooling_rate/cooling.yml
+++ /dev/null
@@ -1,41 +0,0 @@
-# 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
diff --git a/examples/cooling_rate/plot_cooling.py b/examples/cooling_rate/plot_cooling.py
index 8ea0e9a..884f059 100644
--- a/examples/cooling_rate/plot_cooling.py
+++ b/examples/cooling_rate/plot_cooling.py
@@ -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)
diff --git a/pyswiftsim/__init__.py b/pyswiftsim/__init__.py
index 8c6a574..f726238 100644
--- a/pyswiftsim/__init__.py
+++ b/pyswiftsim/__init__.py
@@ -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)
diff --git a/test/test_cooling.py b/test/test_cooling.py
deleted file mode 100644
index 94e3a46..0000000
--- a/test/test_cooling.py
+++ /dev/null
@@ -1,86 +0,0 @@
-#!/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.libcooling import coolingRate
-from pyswiftsim import downloadCoolingTable
-import numpy as np
-from astropy import units, constants
-import yaml
-
-# define parameters
-filename = "test_cooling.yml"
-# adiabatic index
-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))
diff --git a/tests/generate_parameter_file.py b/tests/generate_parameter_file.py
new file mode 100644
index 0000000..576823a
--- /dev/null
+++ b/tests/generate_parameter_file.py
@@ -0,0 +1,85 @@
+#!/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))
diff --git a/tests/test_cooling.py b/tests/test_cooling.py
new file mode 100644
index 0000000..64c489b
--- /dev/null
+++ b/tests/test_cooling.py
@@ -0,0 +1,91 @@
+#!/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.libcooling import coolingRate
+import pyswiftsim
+import numpy as np
+from astropy import units, constants
+import yaml
+
+# define parameters
+filename = "test_cooling.yml"
+# adiabatic index
+gamma = 5. / 3.
+
+# download cooling table
+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)
diff --git a/test/test_cooling.yml b/tests/test_cooling.yml
similarity index 100%
rename from test/test_cooling.yml
rename to tests/test_cooling.yml
--
GitLab