diff --git a/examples/cooling_box/cooling.yml b/examples/cooling_box/cooling.yml deleted file mode 100644 index 66fa4b67f5245eb33a07fae082c870ea3b47535f..0000000000000000000000000000000000000000 --- 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 6d8485ec05c2f292d82c0772c08db97681bec21a..18bbb26903a25e87667887dc5beeacc8f3feb976 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 8fa0e720004ba4101d496108db29ef32f3722ef2..0000000000000000000000000000000000000000 --- 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 8ea0e9ac66756031d151f467edf1bf48cf5c5a2c..884f0592da44acc0444e14e41ea63dcba67bf4eb 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 8c6a5748ae78c94a83aee7c1d808247ac8fd102e..f72623890a318801cc59e92993cbfb17f422ec56 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 94e3a46917f74659c8b8502f2922d316f14d5c41..0000000000000000000000000000000000000000 --- 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 0000000000000000000000000000000000000000..576823abc6a284f8c79281e412f0d0de69dce657 --- /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 0000000000000000000000000000000000000000..64c489befb53c8f6eb912ef194cdc865bdffc24b --- /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