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

Add missed files

parent a25d70c9
Branches
No related tags found
1 merge request!7Change clone to fetch
Examples
========
Here you will find the image generated by every example in the ``examples`` directory.
If you wish to generate them, you will need to first install PySWIFTsim and then run the corresponding example.
.. toctree::
:maxdepth: 2
examples/cooling.rst
Cooling Examples
================
Cooling Rate
~~~~~~~~~~~~
This example is generated in ``examples/cooling_rate`` with the Grackle cooling.
The plot are the same than in [Smith2016]_.
The code has two different modes 1D or 2D.
In the first mode, the code generates particles at different temperatures but same density and computes the cooling rate associated to theses conditions.
.. image:: https://obswww.unige.ch/~lhausamm/PySWIFTsim/examples/cooling_rate.png
In the second mode, the code generates a grid of particles at different density and temperatures and then computes the cooling.
.. image:: https://obswww.unige.ch/~lhausamm/PySWIFTsim/examples/cooling_rate_2d.png
Cooling Box
~~~~~~~~~~~
.. image:: https://obswww.unige.ch/~lhausamm/PySWIFTsim/examples/cooling_box.png
.. [Smith2016] `Grackle: a Chemistry and Cooling Library for Astrophysics <https://arxiv.org/abs/1610.09591>`_
# 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: 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
OutputMode: 0 # Write in output corresponding primordial chemistry mode
MaxSteps: 1000
ConvergenceLimit: 1e-2
GearChemistry:
InitialMetallicity: 0.01295
\ No newline at end of file
#!/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
import matplotlib.pyplot as plt
import yaml
from astropy import units, constants
from time import strftime, gmtime
plt.rc('text', usetex=True)
downloadCoolingTable()
# PARAMETERS
# swift params filename
filename = "cooling.yml"
with_cosmo = False
# particle data
# in atom/cc
part_density = np.array([1.])
part_temperature = np.array([1e4]) # in K
gamma = 5./3.
# time in Myr
t_end = 1.
N_time = 100000
# SCRIPT
def generate_initial_condition(us):
print("Generating default initial conditions")
# get units
unit_length = float(us["UnitLength_in_cgs"])
unit_vel = float(us["UnitVelocity_in_cgs"])
unit_mass = float(us["UnitMass_in_cgs"])
unit_temp = float(us["UnitTemp_in_cgs"])
d = {}
# Deal with units
m_p = constants.m_p.to("g") / (unit_mass * units.g)
rho = part_density * unit_length**3 * m_p
d["density"] = rho.to("")
unit_time = unit_length / unit_vel
unit_energy = unit_mass * unit_length**2
unit_energy /= unit_time**2
k_B = constants.k_B.to("erg / K") / (unit_energy / unit_temp)
k_B *= units.K / units.erg
energy = k_B * part_temperature / unit_temp
energy /= (gamma - 1.) * m_p
d["energy"] = energy
t = t_end * units.Myr.to("s") / unit_time
d["time"] = np.linspace(0, t, N_time)
return d
def plot_solution(energy, data, us):
print("Plotting solution")
rho = data["density"]
time = data["time"]
unit_length = float(us["UnitLength_in_cgs"])
unit_vel = float(us["UnitVelocity_in_cgs"])
unit_mass = float(us["UnitMass_in_cgs"])
unit_time = unit_length / unit_vel
# change units => cgs
rho *= unit_mass / unit_length**3
energy *= unit_length**2 / unit_time**2
time *= unit_time
# do plot
plt.figure()
date = strftime("%d %b %Y", gmtime())
plt.title("Date: {}".format(date))
plt.plot(time, energy)
plt.xlabel("Time [s]")
plt.ylabel("Energy [erg]")
plt.show()
def getParser(filename):
with open(filename, "r") as stream:
stream = yaml.load(stream)
return stream
if __name__ == "__main__":
us = getParser(filename)["InternalUnitSystem"]
d = generate_initial_condition(us)
t = d["time"]
dt = t[1] - t[0]
# 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
plot_solution(energy, d, us)
#!/bin/bash
./plot_cooling.py
# 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-4
MaxSteps: 20000
Omega: 0.8
GearChemistry:
InitialMetallicity: 0.01295
\ No newline at end of file
#!/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
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
plt.rc('text', usetex=True)
downloadCoolingTable()
# PARAMETERS
# cosmology
with_cosmo = False
# swift params filename
filename = "cooling.yml"
# hydrogen mass fraction
h_mass_frac = 0.76
# density in atom / cm3
N_rho = 100
if N_rho == 1:
density = np.array([1.])
else:
density = np.logspace(-6, 4, N_rho)
# temperature in K
N_T = 1000
temperature = np.logspace(1, 9, N_T)
# time step in s
dt = units.Myr * 1e-8
dt = dt.to("s") / units.s
# adiabatic index
gamma = 5. / 3.
def mu(T):
# define constants
T_transition = 1e4
mu_neutral = 4. / (1. + 3. * h_mass_frac)
mu_ion = 4. / (8. - 5. * (1 - h_mass_frac))
# Find correct mu
ind = T > T_transition
mu = np.zeros(T.shape)
mu[ind] = mu_ion
mu[~ind] = mu_neutral
return mu
# SCRIPT
def generate_initial_condition(us):
print("Generating default initial conditions")
# get units
unit_length = float(us["UnitLength_in_cgs"])
unit_vel = float(us["UnitVelocity_in_cgs"])
unit_mass = float(us["UnitMass_in_cgs"])
unit_temp = float(us["UnitTemp_in_cgs"])
d = {}
# generate grid
rho, T = np.meshgrid(density, temperature)
rho = deepcopy(rho.reshape(-1))
T = T.reshape(-1)
d["temperature"] = T
# Deal with units
m_p = constants.m_p.to("g") / (unit_mass * units.g)
rho = rho * unit_length**3 * m_p
d["density"] = rho.to("")
unit_time = unit_length / unit_vel
unit_energy = unit_mass * unit_length**2
unit_energy /= unit_time**2
k_B = constants.k_B.to("erg / K") / (unit_energy / unit_temp)
k_B *= units.K / units.erg
energy = k_B * T / (unit_temp * mu(T))
energy /= (gamma - 1.) * m_p
d["energy"] = energy.to("")
time_step = dt / unit_time
d["time_step"] = time_step
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"]
rho = data["density"]
T = data["temperature"]
T *= float(us["UnitTemp_in_cgs"])
# change units => cgs
u_m = float(us["UnitMass_in_cgs"])
u_l = float(us["UnitLength_in_cgs"])
u_v = float(us["UnitVelocity_in_cgs"])
u_t = u_l / u_v
rho *= u_m / u_l**3
energy *= u_v**2
rate *= u_v**2 / u_t
# lambda cooling
lam = rate * rho
# do plot
if N_rho == 1:
# plot Lambda vs T
plt.figure()
plt.loglog(T, np.abs(lam),
label="SWIFT")
plt.xlabel("Temperature [K]")
plt.ylabel("$\\Lambda$ [erg s$^{-1}$ cm$^{3}$]")
else:
m_p = constants.m_p.to("g") / units.g
rho /= m_p
shape = [N_rho, N_T]
cooling_time = energy / rate
cooling_length = np.sqrt(gamma * (gamma-1.) * energy) * cooling_time
cooling_length = np.log10(np.abs(cooling_length) / units.kpc.to('cm'))
# reshape
rho = rho.reshape(shape)
T = T.reshape(shape)
energy = energy.reshape(shape)
cooling_length = cooling_length.reshape(shape)
_min = -7
_max = 7
N_levels = 100
levels = np.linspace(_min, _max, N_levels)
plt.figure()
plt.contourf(rho, T, cooling_length, levels)
plt.xlabel("Density [atom/cm3]")
plt.ylabel("Temperature [K]")
ax = plt.gca()
ax.set_xscale("log")
ax.set_yscale("log")
cbar = plt.colorbar()
tc = np.arange(_min, _max, 1.5)
cbar.set_ticks(tc)
cbar.set_ticklabels(tc)
cbar.ax.set_ylabel("Log( Cooling Length / kpc )")
date = strftime("%d %b %Y", gmtime())
plt.title("Date: {}".format(date))
plt.show()
if __name__ == "__main__":
parser = getParser(filename)
us = parser["InternalUnitSystem"]
d = generate_initial_condition(us)
# du / dt
print("Computing cooling...")
rate = coolingRate(filename,
d["density"].astype(np.float32),
d["energy"].astype(np.float32),
with_cosmo, d["time_step"])
plot_solution(rate, d, us)
#!/bin/bash
./plot_cooling.py
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment