From e4527a46078e23ff68b88f8760b35b55d6133255 Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Wed, 22 May 2019 17:24:55 +0200
Subject: [PATCH] Add missed files

---
 docs/RST/examples.rst                 |  10 ++
 docs/RST/examples/cooling.rst         |  27 ++++
 examples/cooling_box/cooling.yml      |  28 ++++
 examples/cooling_box/plot_cooling.py  | 136 +++++++++++++++++
 examples/cooling_box/run.sh           |   3 +
 examples/cooling_rate/cooling.yml     |  41 +++++
 examples/cooling_rate/plot_cooling.py | 206 ++++++++++++++++++++++++++
 examples/cooling_rate/run.sh          |   3 +
 8 files changed, 454 insertions(+)
 create mode 100644 docs/RST/examples.rst
 create mode 100644 docs/RST/examples/cooling.rst
 create mode 100644 examples/cooling_box/cooling.yml
 create mode 100644 examples/cooling_box/plot_cooling.py
 create mode 100644 examples/cooling_box/run.sh
 create mode 100644 examples/cooling_rate/cooling.yml
 create mode 100644 examples/cooling_rate/plot_cooling.py
 create mode 100644 examples/cooling_rate/run.sh

diff --git a/docs/RST/examples.rst b/docs/RST/examples.rst
new file mode 100644
index 0000000..0912efd
--- /dev/null
+++ b/docs/RST/examples.rst
@@ -0,0 +1,10 @@
+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
diff --git a/docs/RST/examples/cooling.rst b/docs/RST/examples/cooling.rst
new file mode 100644
index 0000000..1cc6f4d
--- /dev/null
+++ b/docs/RST/examples/cooling.rst
@@ -0,0 +1,27 @@
+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>`_
diff --git a/examples/cooling_box/cooling.yml b/examples/cooling_box/cooling.yml
new file mode 100644
index 0000000..415f22c
--- /dev/null
+++ b/examples/cooling_box/cooling.yml
@@ -0,0 +1,28 @@
+# 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
diff --git a/examples/cooling_box/plot_cooling.py b/examples/cooling_box/plot_cooling.py
new file mode 100644
index 0000000..13c2e4e
--- /dev/null
+++ b/examples/cooling_box/plot_cooling.py
@@ -0,0 +1,136 @@
+#!/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)
diff --git a/examples/cooling_box/run.sh b/examples/cooling_box/run.sh
new file mode 100644
index 0000000..ca4ed09
--- /dev/null
+++ b/examples/cooling_box/run.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+./plot_cooling.py
diff --git a/examples/cooling_rate/cooling.yml b/examples/cooling_rate/cooling.yml
new file mode 100644
index 0000000..166fa04
--- /dev/null
+++ b/examples/cooling_rate/cooling.yml
@@ -0,0 +1,41 @@
+# 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
diff --git a/examples/cooling_rate/plot_cooling.py b/examples/cooling_rate/plot_cooling.py
new file mode 100644
index 0000000..55d21ee
--- /dev/null
+++ b/examples/cooling_rate/plot_cooling.py
@@ -0,0 +1,206 @@
+#!/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)
diff --git a/examples/cooling_rate/run.sh b/examples/cooling_rate/run.sh
new file mode 100644
index 0000000..ca4ed09
--- /dev/null
+++ b/examples/cooling_rate/run.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+./plot_cooling.py
-- 
GitLab