diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 6251ddd24e369f246e64c15b92c399bcf5dca2a0..8098ff3743c52eb01cb9537a59ab734e07929c58 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,11 +1,40 @@
+# you can find different images here https://hub.docker.com/
+image: ubuntu:latest
+
+
before_script:
- # compile swift
- - git clone https://gitlab.cosma.dur.ac.uk/swift/swiftsim.git
+ - apt-get update -qq && apt-get install -y -qq build-essential libhdf5-serial-dev gfortran csh git libtool-bin autoconf
+ - export LHOME=`pwd`
+# Compile grackle
+ - mkdir grackle
+ - git clone https://github.com/grackle-project/grackle.git
+ - cd grackle
+ - ./configure
+ - cd src/clib
+ - sed -i 's/LOCAL_HDF5_INSTALL = \/usr\/local\/hdf5\/1.8.2s/LOCAL_HDF5_INSTALL = \/usr\/lib\/x86_64-linux-gnu\/hdf5\/serial\/lib/' Make.mach.linux-gnu
+ - sed -i 's/MACH_INSTALL_PREFIX = $(HOME)\/local/MACH_INSTALL_PREFIX = $(LHOME)\/grackle/' Make.mach.linux-gnu
+ - make machine-linux-gnu
+ - make && make install
+ - export GRACKLE_ROOT=$LHOME/grackle
+ - cd $LHOME
+
+# compile swift
+ - git clone https://gitlab.cosma.dur.ac.uk/swift/swiftsim.git --single-branch --branch master
- cd swiftsim
- - ./autogen
- - ./configure --with-subgrid=GEAR
+ - ./autogen.sh
+ - ./configure --with-subgrid=GEAR --with-grackle=$GRACKLE_ROOT
- make -j 4
+ - cd ..
+
+# compile pyswiftsim
+ - python setup.py install --user
+
+# compile the docs
+ - cd docs
+ - make html
-test_build:
+test:
script:
- - python setup.py build
\ No newline at end of file
+ - cd test
+ - ./getCoolingTable.sh
+ - ./test_cooling.py
\ No newline at end of file
diff --git a/README b/README
deleted file mode 120000
index 602d87dfd9405502dcd988b9dab5bb724dc4471d..0000000000000000000000000000000000000000
--- a/README
+++ /dev/null
@@ -1 +0,0 @@
-Readme.md
\ No newline at end of file
diff --git a/Readme.md b/Readme.md
index 35a4f470363f3463f9600b28f85ce3bf5a6d001c..4f4c47bc88a31b1aadd99adc92d102415c503904 100644
--- a/Readme.md
+++ b/Readme.md
@@ -1,72 +1,16 @@
+[](https://gitlab.com/clastro/pyswiftsim/commits/master)
+[](https://pyswiftsim.readthedocs.io/en/latest/?badge=latest)
+
Welcome to PySWIFTsim
=====================
This code is a python wrapper around the cosmological hydrodynamical code SWIFT hosted at https://gitlab.cosma.dur.ac.uk/swift/swiftsim.
-
-A few examples are given in test and the API is trying to follow the SWIFT API, therefore a developer of SWIFT should be able to quickly learn how to use it.
+You can find a few examples in test and examples and a documentation in docs.
Install
=======
To install PySWIFT, you can run the following command: `python3 setup.py install --with-swift SWIFT_PATH --user`.
The setup script will read all your parameters in the config file of swift, therefore you need to reinstall PySWIFT with each configuration of SWIFT.
-You may need to remove the build directory to recompile PySWIFT.
-
-
-Adding a new structure
-======================
-
-In pyswiftsim/structure.py, a list of all SWIFT structure are described in python.
-When adding a new structure, a new class should be added to this file using the following example:
-
-```
-class SwiftParams(SwiftStruct):
- """ All object should inherit from SwiftStruct """
- _format = "{sec}s{data}sii{line_size}c".format(
- sec=struct.calcsize(Section._format)*PARSER_MAX_NO_OF_SECTIONS,
- data=struct.calcsize(Parameter._format)*PARSER_MAX_NO_OF_PARAMS,
- line_size=PARSER_MAX_LINE_SIZE
- )
- """ _format describes the type of the object contained in a C structure.
- see struct python package for a list of all possible types.
- For an undefined/complex object, the letter 's' is used and for string 'c'.
- """
-
- _name = [
- "section",
- "data_params",
- "sectionCount",
- "paramCount",
- "filename"
- ]
- """ _name defines the name of each variable """
-
- def __init__(self, data, parent=None):
- """ defines the constructor
- data is a list of bytes and parent is the structure in which this one is included
- """
- super().__init__(self.struct_format, data, parent)
-
- @property
- def struct_substruct(self):
- """ This method should be defined if a substructure is present in this one.
- It returns a dictionary where the key is the name of the object (defined in _name)
- and the value is the array size (1 if this structure is not in an array).
- """
- sec = {
- "class": Section,
- "size": PARSER_MAX_NO_OF_SECTIONS
- }
-
- param = {
- "class": Parameter,
- "size": PARSER_MAX_NO_OF_PARAMS
- }
- return {
- "section": sec,
- "data_params": param
- }
-```
-
-If the new structure is an option (e.g. Part, Cooling, ...) you can use an if to define the _format and _name variables.
\ No newline at end of file
+You may need to remove the build directory to recompile PySWIFT.
\ No newline at end of file
diff --git a/docs/Makefile b/docs/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..298ea9e213e8c4c11f0431077510d4e325733c65
--- /dev/null
+++ b/docs/Makefile
@@ -0,0 +1,19 @@
+# Minimal makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line.
+SPHINXOPTS =
+SPHINXBUILD = sphinx-build
+SOURCEDIR = .
+BUILDDIR = _build
+
+# Put it first so that "make" without argument is like "make help".
+help:
+ @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
+
+.PHONY: help Makefile
+
+# Catch-all target: route all unknown targets to Sphinx using the new
+# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS).
+%: Makefile
+ @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
\ No newline at end of file
diff --git a/docs/RST/examples.rst b/docs/RST/examples.rst
new file mode 100644
index 0000000000000000000000000000000000000000..0912efdbcafde647c6059e69b21fbedb34a6d76f
--- /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 0000000000000000000000000000000000000000..1cc6f4db5c44493cdf7792754084760d5b89954b
--- /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/docs/RST/libcooling.rst b/docs/RST/libcooling.rst
new file mode 100644
index 0000000000000000000000000000000000000000..b83bf4409a206c2f0e27c83b51bb0a1fc7b2382a
--- /dev/null
+++ b/docs/RST/libcooling.rst
@@ -0,0 +1,11 @@
+Libcooling
+==========
+
+.. automodule:: pyswiftsim.libcooling
+ :members:
+ :undoc-members:
+ :private-members:
+ :special-members:
+ :show-inheritance:
+ :inherited-members:
+
diff --git a/docs/RST/pyswiftsim.rst b/docs/RST/pyswiftsim.rst
new file mode 100644
index 0000000000000000000000000000000000000000..0e77bb58ed4a0765bdd419362edebdda15f60178
--- /dev/null
+++ b/docs/RST/pyswiftsim.rst
@@ -0,0 +1,11 @@
+PySWIFTsim
+==========
+
+.. automodule:: pyswiftsim
+ :members:
+ :undoc-members:
+ :private-members:
+ :special-members:
+ :show-inheritance:
+ :inherited-members:
+
diff --git a/docs/conf.py b/docs/conf.py
new file mode 100644
index 0000000000000000000000000000000000000000..8d40b6e2f01cd0c989474cbfab4425561875bd84
--- /dev/null
+++ b/docs/conf.py
@@ -0,0 +1,203 @@
+# -*- coding: utf-8 -*-
+#
+# Configuration file for the Sphinx documentation builder.
+#
+# This file does only contain a selection of the most common options. For a
+# full list see the documentation:
+# http://www.sphinx-doc.org/en/master/config
+
+# -- Path setup --------------------------------------------------------------
+
+# If extensions (or modules to document with autodoc) are in another directory,
+# add these directories to sys.path here. If the directory is relative to the
+# documentation root, use os.path.abspath to make it absolute, like shown here.
+#
+# import os
+# import sys
+# sys.path.insert(0, os.path.abspath('.'))
+
+
+# -- Project information -----------------------------------------------------
+
+project = 'PySWIFTsim'
+copyright = '2019, Loic Hausammann'
+author = 'Loic Hausammann'
+
+# The short X.Y version
+version = ''
+# The full version, including alpha/beta/rc tags
+release = ''
+
+
+# -- General configuration ---------------------------------------------------
+
+# If your documentation needs a minimal Sphinx version, state it here.
+#
+# needs_sphinx = '1.0'
+
+# Add any Sphinx extension module names here, as strings. They can be
+# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
+# ones.
+extensions = [
+ 'sphinx.ext.autodoc',
+ 'sphinx.ext.doctest',
+ 'sphinx.ext.intersphinx',
+ 'sphinx.ext.todo',
+ 'sphinx.ext.coverage',
+ 'sphinx.ext.mathjax',
+ 'sphinx.ext.viewcode',
+ 'numpydoc'
+]
+
+# Add any paths that contain templates here, relative to this directory.
+templates_path = ['_templates']
+
+# The suffix(es) of source filenames.
+# You can specify multiple suffix as a list of string:
+#
+# source_suffix = ['.rst', '.md']
+source_suffix = '.rst'
+
+# The master toctree document.
+master_doc = 'index'
+
+# The language for content autogenerated by Sphinx. Refer to documentation
+# for a list of supported languages.
+#
+# This is also used if you do content translation via gettext catalogs.
+# Usually you set "language" from the command line for these cases.
+language = None
+
+# List of patterns, relative to source directory, that match files and
+# directories to ignore when looking for source files.
+# This pattern also affects html_static_path and html_extra_path.
+exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store']
+
+# The name of the Pygments (syntax highlighting) style to use.
+pygments_style = None
+
+
+# -- Options for HTML output -------------------------------------------------
+
+# The theme to use for HTML and HTML Help pages. See the documentation for
+# a list of builtin themes.
+#
+html_theme = 'sphinx_rtd_theme'
+
+# Theme options are theme-specific and customize the look and feel of a theme
+# further. For a list of options available for each theme, see the
+# documentation.
+#
+# html_theme_options = {}
+
+# Add any paths that contain custom static files (such as style sheets) here,
+# relative to this directory. They are copied after the builtin static files,
+# so a file named "default.css" will overwrite the builtin "default.css".
+html_static_path = ['_static']
+
+# Custom sidebar templates, must be a dictionary that maps document names
+# to template names.
+#
+# The default sidebars (for documents that don't match any pattern) are
+# defined by theme itself. Builtin themes are using these templates by
+# default: ``['localtoc.html', 'relations.html', 'sourcelink.html',
+# 'searchbox.html']``.
+#
+# html_sidebars = {}
+
+
+# -- Options for HTMLHelp output ---------------------------------------------
+
+# Output file base name for HTML help builder.
+htmlhelp_basename = 'PySWIFTsimdoc'
+
+
+# -- Options for LaTeX output ------------------------------------------------
+
+latex_elements = {
+ # The paper size ('letterpaper' or 'a4paper').
+ #
+ # 'papersize': 'letterpaper',
+
+ # The font size ('10pt', '11pt' or '12pt').
+ #
+ # 'pointsize': '10pt',
+
+ # Additional stuff for the LaTeX preamble.
+ #
+ # 'preamble': '',
+
+ # Latex figure (float) alignment
+ #
+ # 'figure_align': 'htbp',
+}
+
+# Grouping the document tree into LaTeX files. List of tuples
+# (source start file, target name, title,
+# author, documentclass [howto, manual, or own class]).
+latex_documents = [
+ (master_doc, 'PySWIFTsim.tex', 'PySWIFTsim Documentation',
+ 'Loic Hausammann', 'manual'),
+]
+
+
+# -- Options for manual page output ------------------------------------------
+
+# One entry per manual page. List of tuples
+# (source start file, name, description, authors, manual section).
+man_pages = [
+ (master_doc, 'pyswiftsim', 'PySWIFTsim Documentation',
+ [author], 1)
+]
+
+
+# -- Options for Texinfo output ----------------------------------------------
+
+# Grouping the document tree into Texinfo files. List of tuples
+# (source start file, target name, title, author,
+# dir menu entry, description, category)
+texinfo_documents = [
+ (master_doc, 'PySWIFTsim', 'PySWIFTsim Documentation',
+ author, 'PySWIFTsim', 'One line description of project.',
+ 'Miscellaneous'),
+]
+
+
+# -- Options for Epub output -------------------------------------------------
+
+# Bibliographic Dublin Core info.
+epub_title = project
+
+# The unique identifier of the text. This can be a ISBN number
+# or the project homepage.
+#
+# epub_identifier = ''
+
+# A unique identification for the text.
+#
+# epub_uid = ''
+
+# A list of files that should not be packed into the epub file.
+epub_exclude_files = ['search.html']
+
+
+# -- Extension configuration -------------------------------------------------
+
+# -- Options for intersphinx extension ---------------------------------------
+
+# Example configuration for intersphinx: refer to the Python standard library.
+intersphinx_mapping = {
+ 'python': ('https://docs.python.org/', None),
+ 'numpy': ('https://docs.scipy.org/doc/numpy', None),
+}
+
+# -- Options for todo extension ----------------------------------------------
+
+# If true, `todo` and `todoList` produce output, else they produce nothing.
+todo_include_todos = True
+
+# -- Options for numpy extension ----------------------------------------------
+
+numpydoc_use_plots = True
+
+numpydoc_xref_param_type = True
diff --git a/docs/index.rst b/docs/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..d11523e49bcf892c03137a516830300e753817cc
--- /dev/null
+++ b/docs/index.rst
@@ -0,0 +1,26 @@
+.. PySWIFTsim documentation master file, created by
+ sphinx-quickstart on Wed May 22 08:34:27 2019.
+ You can adapt this file completely to your liking, but it should at least
+ contain the root `toctree` directive.
+
+Welcome to PySWIFTsim's documentation!
+======================================
+
+PySWIFTsim is a python wrapper around the gravity and SPH solver SWIFT.
+The aim is to provide a python tool that benefits from HPC code and allows to directly use some feature of SWIFT in order to either test it or analyze a simulation.
+
+.. toctree::
+ :maxdepth: 2
+
+ RST/pyswiftsim.rst
+ RST/libcooling.rst
+ RST/examples.rst
+
+
+
+Indices and tables
+==================
+
+* :ref:`genindex`
+* :ref:`modindex`
+* :ref:`search`
diff --git a/docs/make.bat b/docs/make.bat
new file mode 100644
index 0000000000000000000000000000000000000000..27f573b87af11e2cbbd9f54eb1ee285a58550146
--- /dev/null
+++ b/docs/make.bat
@@ -0,0 +1,35 @@
+@ECHO OFF
+
+pushd %~dp0
+
+REM Command file for Sphinx documentation
+
+if "%SPHINXBUILD%" == "" (
+ set SPHINXBUILD=sphinx-build
+)
+set SOURCEDIR=.
+set BUILDDIR=_build
+
+if "%1" == "" goto help
+
+%SPHINXBUILD% >NUL 2>NUL
+if errorlevel 9009 (
+ echo.
+ echo.The 'sphinx-build' command was not found. Make sure you have Sphinx
+ echo.installed, then set the SPHINXBUILD environment variable to point
+ echo.to the full path of the 'sphinx-build' executable. Alternatively you
+ echo.may add the Sphinx directory to PATH.
+ echo.
+ echo.If you don't have Sphinx installed, grab it from
+ echo.http://sphinx-doc.org/
+ exit /b 1
+)
+
+%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS%
+goto end
+
+:help
+%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS%
+
+:end
+popd
diff --git a/examples/cooling_box/cooling.yml b/examples/cooling_box/cooling.yml
new file mode 100644
index 0000000000000000000000000000000000000000..415f22ce7fadc74f042ad91d30079aa8ce0ffbea
--- /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 0000000000000000000000000000000000000000..13c2e4e514fd14b4f749a8a7d7f39a1e19414088
--- /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 0000000000000000000000000000000000000000..ca4ed090d7161d18ce520af16a81caaaebd04cf3
--- /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 0000000000000000000000000000000000000000..166fa0401a5870b941e503152d9843a8e2417e69
--- /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 0000000000000000000000000000000000000000..55d21ee1790177aa1755a7cc2815cdb26414502a
--- /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 0000000000000000000000000000000000000000..ca4ed090d7161d18ce520af16a81caaaebd04cf3
--- /dev/null
+++ b/examples/cooling_rate/run.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+./plot_cooling.py
diff --git a/pyswiftsim/__init__.py b/pyswiftsim/__init__.py
index 526dbce334dd759ae3af7299c78404f02c44fa8d..8c6a5748ae78c94a83aee7c1d808247ac8fd102e 100644
--- a/pyswiftsim/__init__.py
+++ b/pyswiftsim/__init__.py
@@ -1,19 +1,30 @@
-"""
-This file is part of PYSWIFT.
-Copyright (c) 2018 Loic Hausammann (loic.hausammann@epfl.ch)
+# ########################################################################
+# 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/>.
+# ########################################################################
-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.
+import os
+import urllib.request
-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/>.
-"""
+def downloadCoolingTable():
+ url = "http://virgodb.cosma.dur.ac.uk/"
+ url += "swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5"
+ filename = "CloudyData_UVB=HM2012.h5"
-from pyswiftsim import *
+ if not os.path.isfile(filename):
+ # Download the file from `url` and save it locally under `file_name`:
+ urllib.request.urlretrieve(url, filename)
diff --git a/setup.py b/setup.py
index 919adee2849f4fc5dcc55e8dd88d25a6b2ccf337..1c76bc94e8463de8a2bf4af34c8a9c61a2505365 100644
--- a/setup.py
+++ b/setup.py
@@ -118,7 +118,7 @@ data_files = []
##############
c_src = glob("src/*.c")
-ext_modules = Extension("pyswiftsim.wrapper",
+ext_modules = Extension("pyswiftsim.libcooling",
c_src,
include_dirs=include,
libraries=lib,
diff --git a/src/Makefile.am b/src/Makefile.am
deleted file mode 100644
index a77473fd9cf45fb642c0175a3d22490b9b7f9084..0000000000000000000000000000000000000000
--- a/src/Makefile.am
+++ /dev/null
@@ -1,68 +0,0 @@
-# This file is part of PySWIFTsim.
-# Copyright (c) 2019 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 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 General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-# Add the non-standard paths to the included library headers
-AM_CFLAGS = -I$(top_srcdir)/../src $(HDF5_CPPFLAGS) $(GSL_INCS) $(FFTW_INCS) $(NUMA_INCS) $(GRACKLE_INCS) $(PYTHON_INCS)
-
-# Assign a "safe" version number
-AM_LDFLAGS = $(HDF5_LDFLAGS) $(FFTW_LIBS) -version-info 0:0:0
-
-# The git command, if available.
-GIT_CMD = @GIT_CMD@
-
-# Additional dependencies for shared libraries.
-EXTRA_LIBS = $(HDF5_LIBS) $(FFTW_LIBS) $(NUMA_LIBS) $(PROFILER_LIBS) $(TCMALLOC_LIBS) \
- $(JEMALLOC_LIBS) $(TBBMALLOC_LIBS) $(GRACKLE_LIBS) $(GSL_LIBS) $(PYTHON_LIBS) \
- $(top_srcdir)/src/.libs
-
-# MPI libraries.
-MPI_LIBS = $(PARMETIS_LIBS) $(METIS_LIBS) $(MPI_THREAD_LIBS)
-MPI_FLAGS = -DWITH_MPI $(PARMETIS_INCS) $(METIS_INCS)
-
-# Build the libswiftsim library
-lib_LTLIBRARIES = libpyswiftsim.la
-# Build a MPI-enabled version too?
-if HAVEMPI
-lib_LTLIBRARIES += libpyswiftsim_mpi.la
-endif
-
-# List required headers
-include_HEADERS = chemistry_wrapper.h cooling_wrapper.h \
- cosmology_wrapper.h parser_wrapper.h part_wrapper.h \
- pyswiftsim_tools.h units_wrapper.h
-
-
-# Common source files
-AM_SOURCES = chemistry_wrapper.c cooling_wrapper.c cosmology_wrapper.c \
- parser_wrapper.c part_wrapper.c pyswiftsim_tools.c \
- units_wrapper.c wrapper.c
-
-# Include files for distribution, not installation.
-nobase_noinst_HEADERS =
-
-# Sources and flags for regular library
-libpyswiftsim_la_SOURCES = $(AM_SOURCES)
-libpyswiftsim_la_CFLAGS = $(AM_CFLAGS)
-libpyswiftsim_la_LDFLAGS = $(AM_LDFLAGS) $(EXTRA_LIBS) -lswiftsim
-libpyswiftsim_la_LIBADD = $(GRACKLE_LIBS) $(VELOCIRAPTOR_LIBS)
-
-# Sources and flags for MPI library
-libpyswiftsim_mpi_la_SOURCES = $(AM_SOURCES)
-libpyswiftsim_mpi_la_CFLAGS = $(AM_CFLAGS) $(MPI_FLAGS)
-libpyswiftsim_mpi_la_LDFLAGS = $(AM_LDFLAGS) $(MPI_LIBS) $(EXTRA_LIBS) -lswiftsim_mpi
-libpyswiftsim_mpi_la_SHORTNAME = mpi
-libpyswiftsim_mpi_la_LIBADD = $(GRACKLE_LIBS) $(VELOCIRAPTOR_LIBS)
-
diff --git a/src/cooling_wrapper.c b/src/cooling_wrapper.c
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..47dd1ccb52857d1794719eaac9ba2bdcedded1c2 100644
--- a/src/cooling_wrapper.c
+++ b/src/cooling_wrapper.c
@@ -0,0 +1,203 @@
+/*******************************************************************************
+ * This file is part of PYSWIFTSIM.
+ * 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/>.
+ *
+ ******************************************************************************/
+#include "cooling_wrapper.h"
+#include "pyswiftsim_tools.h"
+
+/* TODO Remove this */
+struct cooling_function_data cooling;
+int initialized;
+
+
+/**
+ * @brief Set the cooling element fractions
+ *
+ * @param xp The #xpart to set
+ * @param frac The numpy array containing the fractions (id, element)
+ * @param idx The id (in frac) of the particle to set
+ */
+void pycooling_set_fractions(struct xpart *xp, PyArrayObject* frac, const int idx) {
+#ifdef COOLING_GRACKLE
+ struct cooling_xpart_data *data = &xp->cooling_data;
+ data->metal_frac = *(float*)PyArray_GETPTR2(frac, idx, 0);
+
+#ifdef COOLING_GRACKLE
+#if COOLING_GRACKLE_MODE > 0
+ data->HI_frac = *(float*)PyArray_GETPTR2(frac, idx, 1);
+ data->HII_frac = *(float*)PyArray_GETPTR2(frac, idx, 2);
+ data->HeI_frac = *(float*)PyArray_GETPTR2(frac, idx, 3);
+ data->HeII_frac = *(float*)PyArray_GETPTR2(frac, idx, 4);
+ data->HeIII_frac = *(float*)PyArray_GETPTR2(frac, idx, 5);
+ data->e_frac = *(float*)PyArray_GETPTR2(frac, idx, 6);
+#endif // COOLING_GRACKLE_MODE
+#if COOLING_GRACKLE_MODE > 1
+ data->HM_frac = *(float*)PyArray_GETPTR2(frac, idx, 7);
+ data->H2I_frac = *(float*)PyArray_GETPTR2(frac, idx, 8);
+ data->H2II_frac = *(float*)PyArray_GETPTR2(frac, idx, 9);
+#endif // COOLING_GRACKLE_MODE
+#if COOLING_GRACKLE_MODE > 2
+ data->DI_frac = *(float*)PyArray_GETPTR2(frac, idx, 10);
+ data->DII_frac = *(float*)PyArray_GETPTR2(frac, idx, 11);
+ data->HDI_frac = *(float*)PyArray_GETPTR2(frac, idx, 12);
+#endif // COOLING_GRACKLE_MODE
+#endif // COOLING_GRACKLE
+
+#else
+ error("Not implemented");
+#endif
+}
+
+/**
+ * @brief Compute the cooling rate
+ *
+ * args is expecting:
+ * - a parameter filename,
+ * - a numpy array of the densities,
+ * - a numpy array of the specific energies,
+ * - an optional int if the cosmology should be used,
+ * - an optional float for the time step,
+ * - an optional array with the element fractions.
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return cooling rate
+ */
+PyObject* pycooling_rate(PyObject* self, PyObject* args) {
+ import_array();
+
+
+ char *filename;
+
+ PyArrayObject *rho;
+ PyArrayObject *energy;
+ PyArrayObject *fractions = NULL;
+
+ float dt = 1e-3;
+ int with_cosmo = 0;
+
+ /* parse argument */
+ if (!PyArg_ParseTuple(
+ args, "sOO|ifO", &filename, &rho, &energy,
+ &with_cosmo, &dt, &fractions))
+ return NULL;
+
+ /* check numpy array */
+ if (pytools_check_array(energy, 1, NPY_FLOAT, "Energy") != SWIFT_SUCCESS)
+ return NULL;
+
+ if (pytools_check_array(rho, 1, NPY_FLOAT, "Density") != SWIFT_SUCCESS)
+ return NULL;
+
+ if (fractions != NULL &&
+ pytools_check_array(fractions, 2, NPY_FLOAT, "Fractions") != SWIFT_SUCCESS)
+ return NULL;
+
+ if (PyArray_DIM(energy, 0) != PyArray_DIM(rho, 0))
+ error("Density and energy should have the same dimension");
+
+ if (fractions != NULL &&
+ PyArray_DIM(fractions, 0) != PyArray_DIM(rho,0))
+ error("Fractions should have the same first dimension than the others");
+
+ size_t N = PyArray_DIM(energy, 0);
+
+
+ /* Initialize parser */
+ struct swift_params params;
+ parser_read_file(filename, ¶ms);
+
+ /* Init unit system */
+ struct unit_system us;
+ units_init_from_params(&us, ¶ms, "InternalUnitSystem");
+
+ /* Init physical constant */
+ struct phys_const pconst;
+ phys_const_init(&us, ¶ms, &pconst);
+
+ /* init cooling */
+ if (initialized == 0) {
+ // struct cooling_function_data cooling;
+ cooling_init_backend(¶ms, &us, &pconst, &cooling);
+ initialized = 1;
+ }
+ else if (initialized == 1) {
+ message("WARNING: not reinitializing the cooling, need to be fixed");
+ initialized = 2;
+ }
+
+ /* Init chemistry */
+ struct chemistry_global_data chemistry;
+ chemistry_init_backend(¶ms, &us, &pconst, &chemistry);
+
+ /* Init cosmo */
+ struct cosmology cosmo;
+
+ if (with_cosmo)
+ cosmology_init(¶ms, &us, &pconst, &cosmo);
+ else
+ cosmology_init_no_cosmo(&cosmo);
+
+ /* Init hydro prop */
+ struct hydro_props hydro_props;
+ hydro_props_init(&hydro_props, &pconst, &us, ¶ms);
+
+ /* Init entropy floor */
+ struct entropy_floor_properties floor_props;
+ entropy_floor_init(&floor_props, &pconst, &us,
+ &hydro_props, ¶ms);
+
+ struct part p;
+ struct xpart xp;
+
+ hydro_first_init_part(&p, &xp);
+
+ /* return object */
+ PyArrayObject *rate = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(energy), PyArray_DIMS(energy), NPY_FLOAT);
+
+ /* loop over all particles */
+ for(size_t i = 0; i < N; i++)
+ {
+ /* set particle data */
+ p.rho = *(float*) PyArray_GETPTR1(rho, i);
+ float u = *(float*) PyArray_GETPTR1(energy, i);
+ hydro_set_physical_internal_energy(&p, &xp, &cosmo, u);
+
+ if (fractions != NULL)
+ pycooling_set_fractions(&xp, fractions, i);
+ else
+ cooling_first_init_part(&pconst, &us, &cosmo, &cooling, &p, &xp);
+
+ chemistry_first_init_part(&pconst, &us, &cosmo, &chemistry, &p, &xp);
+
+ hydro_set_physical_internal_energy_dt(&p, &cosmo, 0);
+
+ /* compute cooling rate */
+ cooling_cool_part(&pconst, &us, &cosmo, &hydro_props,
+ &floor_props, &cooling, &p, &xp, dt, dt);
+ float *tmp = PyArray_GETPTR1(rate, i);
+ *tmp = hydro_get_physical_internal_energy_dt(&p, &cosmo);
+
+ }
+
+ /* Cleanup */
+ cosmology_clean(&cosmo);
+
+ return (PyObject *) rate;
+
+}
+
diff --git a/src/cooling_wrapper.h b/src/cooling_wrapper.h
index 86841cff6fcfed69a26a4fa89e42bd0ab9d1ab0e..bd9017f8abb4c31084db4ace60fa52372f943006 100644
--- a/src/cooling_wrapper.h
+++ b/src/cooling_wrapper.h
@@ -21,6 +21,11 @@
#include "pyswiftsim_includes.h"
+/* A few global variables */
+extern struct cooling_function_data cooling;
+extern int initialized;
+
+
PyObject* pycooling_rate(PyObject* self, PyObject* args);
#endif // __PYSWIFTSIM_COOLING_H__
diff --git a/src/wrapper.c b/src/libcooling.c
similarity index 67%
rename from src/wrapper.c
rename to src/libcooling.c
index af15b99499ce9bd9ca7ea085da4ad418cf9bc856..6248d437d0ad9519b033c81acf8f80760f54280b 100644
--- a/src/wrapper.c
+++ b/src/libcooling.c
@@ -22,30 +22,38 @@
#include <math.h>
#include <numpy/arrayobject.h>
-
/* definition of the method table */
-static PyMethodDef wrapper_methods[] = {
+static PyMethodDef libcooling_methods[] = {
{"coolingRate", pycooling_rate, METH_VARARGS,
"Compute the cooling rate.\n\n"
"Parameters\n"
"----------\n\n"
- "pyconst: swift physical constant\n"
- "pyus: swift unit system\n"
- "cooling: swift cooling structure\n"
- "rho: np.array\n"
- "\t Mass density in pyus units\n"
- "energy: np.array\n"
- "\t Internal energy in pyus units\n"
+ "filename: str\n"
+ "\t Parameter file\n"
+ "rho: array\n"
+ "\t Mass density in internal units\n"
+ "energy: array\n"
+ "\t Internal energy in internal units\n"
"dt: float, optional\n"
- "\t Time step in pyus units\n"
- "fractions: np.array, optional\n"
+ "\t Time step in internal units\n"
+ "with_cosmo: int, optional\n"
+ "\t Use cosmology?\n"
+ "fractions: array, optional\n"
"\t Fraction of each cooling element (including metals)\n"
+ "\n"
"Returns\n"
- "-------\n\n"
- "rate: np.array\n"
+ "-------\n"
+ "rate: array\n"
"\t Cooling rate\n"
+ "\n"
+ "Examples\n"
+ "--------\n"
+ ">>> rho = np.array([1])\n"
+ ">>> u = np.array([1000])\n"
+ ">>> filename = 'params.yml'\n"
+ ">>> rate = coolingRate(filename, rho, u)\n"
},
@@ -54,12 +62,12 @@ static PyMethodDef wrapper_methods[] = {
-static struct PyModuleDef wrapper_cmodule = {
+static struct PyModuleDef libcooling_cmodule = {
PyModuleDef_HEAD_INIT,
- "wrapper",
- "Wrapper around the SPH cosmological simulation code SWIFT",
+ "libcooling",
+ "Wrapper around the cooling of SWIFT",
-1,
- wrapper_methods,
+ libcooling_methods,
NULL, /* m_slots */
NULL, /* m_traverse */
NULL, /* m_clear */
@@ -67,7 +75,7 @@ static struct PyModuleDef wrapper_cmodule = {
};
-PyMODINIT_FUNC PyInit_libpyswiftsim(void)
+PyMODINIT_FUNC PyInit_libcooling(void)
{
PyObject *m;
@@ -76,7 +84,7 @@ PyMODINIT_FUNC PyInit_libpyswiftsim(void)
import_array();
Py_Initialize();
- m = PyModule_Create(&wrapper_cmodule);
+ m = PyModule_Create(&libcooling_cmodule);
return m;
}
diff --git a/src/pyswiftsim_tools.c b/src/pyswiftsim_tools.c
index 4fd232542387b913e7604f99cb95eb18ae3f8e70..e9c8f3d3196cbca65265e7f9d62073d902c31529 100644
--- a/src/pyswiftsim_tools.c
+++ b/src/pyswiftsim_tools.c
@@ -35,21 +35,21 @@ int pytools_check_array(PyArrayObject *obj, int dim, int type, const char* name)
/* check if array */
if (!PyArray_Check(obj))
{
- pyerror("Expecting a numpy array for %s", name);
+ error("Expecting a numpy array for %s", name);
}
/* check if required dim */
if (PyArray_NDIM(obj) != dim)
{
- pyerror("Array should be a %i dimensional object for %s", dim, name);
+ error("Array should be a %i dimensional object for %s", dim, name);
}
/* check data type */
if (PyArray_TYPE(obj) != type)
{
- pyerror("Wrong array type for %s", name);
+ error("Wrong array type for %s", name);
}
- return SUCCESS;
+ return SWIFT_SUCCESS;
}
diff --git a/src/pyswiftsim_tools.h b/src/pyswiftsim_tools.h
index b7dd07b2b5972cb8b77d1afe0ffdb0b1dc7855fc..567314f028f17bf34e395ea77292763d1373447e 100644
--- a/src/pyswiftsim_tools.h
+++ b/src/pyswiftsim_tools.h
@@ -19,7 +19,7 @@
#ifndef __PYSWIFTSIM_TOOLS_H__
#define __PYSWIFTSIM_TOOLS_H__
-#include "pyswiftsim.h"
+#include "pyswiftsim_includes.h"
#define DIM 3
#define STRING_SIZE 200
@@ -33,7 +33,14 @@
}
#define IMPORT_ARRAY() IMPORT_ARRAY1(NUMPY_IMPORT_ARRAY_RETVAL);
-
+
+/* error code in pyswiftsim */
+enum error_code {
+ SWIFT_FAIL,
+ SWIFT_SUCCESS,
+};
+
+
int pytools_check_array(PyArrayObject *obj, int dim, int type, const char* name);
diff --git a/test/test_cooling.py b/test/test_cooling.py
new file mode 100644
index 0000000000000000000000000000000000000000..94e3a46917f74659c8b8502f2922d316f14d5c41
--- /dev/null
+++ b/test/test_cooling.py
@@ -0,0 +1,86 @@
+#!/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/test/test_cooling.yml b/test/test_cooling.yml
new file mode 100644
index 0000000000000000000000000000000000000000..a30e9edf8a5fa99564e288e5029c038af1fb2b32
--- /dev/null
+++ b/test/test_cooling.yml
@@ -0,0 +1,25 @@
+# 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
+
+# Cooling with Grackle 3.0
+GrackleCooling:
+ CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
+ 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 # (optional) User provide volumetric heating rates
+ ProvideSpecificHeatingRates: 0 # (optional) User provide specific heating rates
+ SelfShieldingMethod: 0 # (optional) Grackle (<= 3) or Gear self shielding method
+ OutputMode: 0 # (optional) Write in output corresponding primordial chemistry mode
+ MaxSteps: 10000 # (optional) Max number of step when computing the initial composition
+ ConvergenceLimit: 1e-2 # (optional) Convergence threshold (relative) for initial composition
+
+# 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.