diff --git a/.gitignore b/.gitignore
index e4e5f6c8b2deb54bf38312dd9e2f53489b60d6a6..68ac6223c076dcd2a3315a51b721e0515a6f2a8a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,2 @@
-*~
\ No newline at end of file
+*~
+*.hdf5
\ No newline at end of file
diff --git a/AUTHORS b/AUTHORS
new file mode 100644
index 0000000000000000000000000000000000000000..d202675ba9a917e0cd9db1c0fd569fc037c444cd
--- /dev/null
+++ b/AUTHORS
@@ -0,0 +1,2 @@
+Loic Hausammann loic.hausammann@epfl.ch
+Yves Revaz yves.revaz@epfl.ch
diff --git a/COPYING b/COPYING
new file mode 120000
index 0000000000000000000000000000000000000000..88798ab60d0071429c816e50551d74789c45a51a
--- /dev/null
+++ b/COPYING
@@ -0,0 +1 @@
+/usr/share/automake-1.15/COPYING
\ No newline at end of file
diff --git a/ChangeLog b/ChangeLog
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/INSTALL b/INSTALL
new file mode 120000
index 0000000000000000000000000000000000000000..ddcdb76b36bd0e6c66dda2e2a0b0f48bfffc7109
--- /dev/null
+++ b/INSTALL
@@ -0,0 +1 @@
+/usr/share/automake-1.15/INSTALL
\ No newline at end of file
diff --git a/Makefile.am b/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..3195eae8221667b0734e984c6c8ee4d21a187856
--- /dev/null
+++ b/Makefile.am
@@ -0,0 +1,21 @@
+# This file is part of PySWIFTsim.
+# Copyright (c) 2019 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/>.
+
+# Automake stuff
+ACLOCAL_AMFLAGS = -I m4
+
+# Show the way...
+SUBDIRS = src
diff --git a/NEWS b/NEWS
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/README b/README
new file mode 120000
index 0000000000000000000000000000000000000000..602d87dfd9405502dcd988b9dab5bb724dc4471d
--- /dev/null
+++ b/README
@@ -0,0 +1 @@
+Readme.md
\ No newline at end of file
diff --git a/Readme.md b/Readme.md
index 4cb9fdfb8e46e71cebb31aff1b97db3b923ae2ee..35a4f470363f3463f9600b28f85ce3bf5a6d001c 100644
--- a/Readme.md
+++ b/Readme.md
@@ -6,6 +6,13 @@ This code is a python wrapper around the cosmological hydrodynamical code SWIFT
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.
+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
======================
diff --git a/configure.ac b/configure.ac
new file mode 100644
index 0000000000000000000000000000000000000000..f9353c57a215829c97255817266d01f20f2d2c5b
--- /dev/null
+++ b/configure.ac
@@ -0,0 +1,145 @@
+# This file is part of SWIFT.
+# Copyright (C) 2012 pedro.gonnet@durham.ac.uk.
+# 2016 p.w.draper@durham.ac.uk.
+#
+# 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/>.
+
+# Init the project.
+AC_INIT([PySWIFTsim],[0.1.0],[https://gitlab.cosma.dur.ac.uk/swift/PySWIFTsim])
+
+swift_config_flags="$*"
+
+# We want to stop when given unrecognised options. No subdirs so this is safe.
+enable_option_checking=${enable_option_checking:-fatal}
+if test -n "$ac_unrecognized_opts"; then
+ case $enable_option_checking in
+ no)
+ ;;
+ fatal)
+ { $as_echo "$as_me: error: unrecognized options: $ac_unrecognized_opts" >&2
+ { (exit 1); exit 1; }; }
+ ;;
+ *)
+ $as_echo "$as_me: WARNING: unrecognized options: $ac_unrecognized_opts" >&2
+ ;;
+ esac
+fi
+
+AC_COPYRIGHT
+AC_CONFIG_SRCDIR([src/wrapper.c])
+AC_CONFIG_AUX_DIR([.])
+AM_INIT_AUTOMAKE([subdir-objects])
+
+# Add local macro collection.
+AC_CONFIG_MACRO_DIR([m4])
+
+# Stop default CFLAGS from anyone except the environment.
+: ${CFLAGS=""}
+
+# Generate header file.
+AM_CONFIG_HEADER(config.h)
+
+# Find and test the compiler.
+AX_CHECK_ENABLE_DEBUG
+AC_PROG_CC
+AM_PROG_CC_C_O
+
+# Enable POSIX and platform extension preprocessor macros.
+AC_USE_SYSTEM_EXTENSIONS
+
+# Check for compiler version and vendor.
+AX_COMPILER_VENDOR
+AX_COMPILER_VERSION
+
+# Restrict support.
+AC_C_RESTRICT
+
+
+# Add libtool support (now that CC is defined).
+LT_INIT
+
+# Need C99 and inline support.
+AC_PROG_CC_C99
+AC_C_INLINE
+
+# Need fortran for grackle and velociraptor
+AC_PROG_FC
+AC_FC_LIBRARY_LDFLAGS
+
+
+# Define HAVE_POSIX_MEMALIGN if it works.
+AX_FUNC_POSIX_MEMALIGN
+
+# Autoconf stuff.
+AC_PROG_INSTALL
+AC_PROG_MAKE_SET
+AC_HEADER_STDC
+
+SWIFT_ADD_LIBRARIES
+
+SWIFT_ADD_SPECIAL_FLAGS
+
+# this one should be at the beginning
+# it overwrites all the required parameters
+SWIFT_ADD_SUBGRID_MODELS
+
+
+
+SWIFT_ADD_GRAVITY
+
+SWIFT_ADD_HYDRO
+
+SWIFT_ADD_STARS
+
+SWIFT_ADD_KERNEL
+
+SWIFT_ADD_DIMENSION
+
+SWIFT_ADD_EQUATION_OF_STATE
+
+SWIFT_ADD_CHEMISTRY
+
+SWIFT_ADD_COOLING
+
+SWIFT_ADD_RIEMANN
+
+SWIFT_ADD_TRACER
+
+PYSWIFT_ADD_PYTHON
+
+# Check for git, needed for revision stamps.
+AC_PATH_PROG([GIT_CMD], [git])
+AC_SUBST([GIT_CMD])
+
+# Make the documentation. Add conditional to handle disable option.
+DX_INIT_DOXYGEN(libswift,doc/Doxyfile,doc/)
+AM_CONDITIONAL([HAVE_DOXYGEN], [test "$ac_cv_path_ac_pt_DX_DOXYGEN" != ""])
+
+# Check if using EAGLE cooling
+AM_CONDITIONAL([HAVEEAGLECOOLING], [test $with_cooling = "EAGLE"])
+
+# Handle .in files.
+AC_CONFIG_FILES([Makefile src/Makefile])
+
+# Save the compilation options
+AC_DEFINE_UNQUOTED([PYSWIFT_CONFIG_FLAGS],["$swift_config_flags"],[Flags passed to configure])
+
+# Make sure the latest git revision string gets included, when we are
+# working in a checked out repository.
+test -d ${srcdir}/.git && touch ${srcdir}/src/version.c
+
+AC_DEFINE([PYSWIFT_PACKAGE_URL],["www.swiftsim.com"], [Package web pages])
+
+# Generate output.
+AC_OUTPUT
diff --git a/examples/cooling_box/CloudyData_UVB=HM2012.h5 b/examples/cooling_box/CloudyData_UVB=HM2012.h5
new file mode 100644
index 0000000000000000000000000000000000000000..1e35a914e33cb6465b48babba8730ba47a92f6f2
Binary files /dev/null and b/examples/cooling_box/CloudyData_UVB=HM2012.h5 differ
diff --git a/examples/cooling_box/cooling.yml b/examples/cooling_box/cooling.yml
new file mode 100644
index 0000000000000000000000000000000000000000..cbcb4c555033d6204cd1fbffcd1a70e613a3d47c
--- /dev/null
+++ b/examples/cooling_box/cooling.yml
@@ -0,0 +1,60 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 2.e33 # 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 governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 5e-2 # The end time of the simulation (in internal units).
+ dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: sedov # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 1e-2 # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-3 # Time between statistics output
+ tasks_per_cell: 1000
+
+# 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.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./sedov.hdf5 # The file to read
+ h_scaling: 3.33
+
+
+# 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
+
+Gravity:
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ theta: 0.7 # Opening angle (Multipole acceptance criterion)
+ epsilon: 0.1 # Softening length (in internal units).
+ a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
+ r_cut_max: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
+ r_cut_min: 0.1 # (Optional) Cut-off in number of top-level cells below which no truncation of FMM forces are performed (this is the default value).
+
+GearChemistry:
+ InitialMetallicity: 0.01295
\ No newline at end of file
diff --git a/examples/cooling_box/generate_grackle_data.py b/examples/cooling_box/generate_grackle_data.py
new file mode 100644
index 0000000000000000000000000000000000000000..e9ed1bb77a4e05a2534773fdbfb8d40d51f8e152
--- /dev/null
+++ b/examples/cooling_box/generate_grackle_data.py
@@ -0,0 +1,262 @@
+#!/usr/bin/env python
+########################################################################
+#
+# Cooling rate data generation. This code is a modified version
+# of a code developed by the Grackle Team.
+#
+#
+# Copyright (c) 2013-2016, Grackle Development Team.
+#
+# Distributed under the terms of the Enzo Public Licence.
+#
+# The full license is in the file LICENSE, distributed with this
+# software.
+########################################################################
+
+"""
+Generate (or update) a hdf5 file containing
+the grackle data for each primordial chemistry
+"""
+
+import numpy as np
+from h5py import File
+
+from pygrackle import \
+ chemistry_data, \
+ setup_fluid_container
+
+from pygrackle.utilities.physical_constants import \
+ mass_hydrogen_cgs
+
+filename = "grackle.hdf5"
+
+# metal_mass_fraction = 1e-10
+metal_mass_fraction = 0.01295
+
+tolerance = 1e-1
+max_iter = 1e8
+
+t_end = 1. # in cooling time
+N_time = 1000
+
+part_density = 1.0 # in atom/cc
+part_temperature = 1e6 # in K
+
+
+def generate_data(primordial_chemistry):
+ """
+ Compute the cooling a single particles and save its evolution
+ to an hdf5 file
+ """
+ current_redshift = 0.
+
+ # Set solver parameters
+ my_chemistry = chemistry_data()
+ my_chemistry.use_grackle = 1
+ my_chemistry.with_radiative_cooling = 1
+ my_chemistry.primordial_chemistry = primordial_chemistry
+ my_chemistry.metal_cooling = 1
+ my_chemistry.UVbackground = 1
+ my_chemistry.self_shielding_method = 0
+ my_chemistry.H2_self_shielding = 0
+ my_chemistry.grackle_data_file = "CloudyData_UVB=HM2012.h5"
+
+ # Set units
+ my_chemistry.comoving_coordinates = 0 # proper units
+ my_chemistry.a_units = 1.0
+ my_chemistry.a_value = 1.0 / (1.0 + current_redshift) / \
+ my_chemistry.a_units
+ my_chemistry.length_units = 3.085e21
+ my_chemistry.density_units = 2.0e33 / my_chemistry.length_units**3
+ my_chemistry.velocity_units = 1e5
+ my_chemistry.time_units = my_chemistry.length_units / \
+ my_chemistry.velocity_units
+
+ density = np.array([part_density * mass_hydrogen_cgs])
+ # Call convenience function for setting up a fluid container.
+ # This container holds the solver parameters, units, and fields.
+ temperature = np.array([part_temperature])
+ fc = setup_fluid_container(my_chemistry,
+ metal_mass_fraction=metal_mass_fraction,
+ temperature=temperature,
+ density=density,
+ converge=True,
+ tolerance=tolerance,
+ max_iterations=max_iter)
+
+ # compute time steps
+ fc.calculate_cooling_time()
+ t_max = t_end * abs(fc["cooling_time"])
+
+ t = np.linspace(0, t_max, N_time)
+
+ f, data = write_input(
+ fc, my_chemistry, primordial_chemistry)
+
+ dt = t[1] - t[0]
+
+ # compute evolution
+ E = np.zeros(t.shape)
+ E[0] = fc["energy"]
+ for i in range(N_time-1):
+ fc.solve_chemistry(dt)
+ E[i+1] = fc["energy"]
+
+ write_output(f, data, fc, t, E)
+
+
+def write_output(f, data, fc, t, E):
+ """
+ Write the output of the simulation
+ """
+ data.create_group("Output")
+ tmp = data["Output"]
+
+ # Energy
+ dset = tmp.create_dataset("Energy", E.shape,
+ dtype=E.dtype)
+ dset[:] = E
+
+ # Time
+ dset = tmp.create_dataset("Time", t.shape,
+ dtype=t.dtype)
+ dset[:] = t
+
+ f.close()
+
+
+def simulationAreEqual(fields, data):
+ """
+ Check if current simulation has already
+ been run
+ """
+ test = True
+ for key in fields.keys():
+ if data.attrs[key] != fields[key]:
+ test = False
+ break
+
+ return test
+
+
+def getSameSimulation(data, density, my_chemistry):
+ """
+ Get indice of the same simulation.
+ returns -1 if not found
+ """
+ fields = {
+ "MetalCooling": my_chemistry.metal_cooling,
+ "UVBackground": my_chemistry.UVbackground,
+ "SelfShieldingMethod": my_chemistry.self_shielding_method,
+ "MetalMassFraction": metal_mass_fraction,
+ "ScaleFactor": my_chemistry.a_value,
+ "Temperature": part_temperature,
+ "Density": density,
+ }
+ ind = 0
+ for d in data:
+ if "Params" in data[d]:
+ tmp = data[d]["Params"]
+ if simulationAreEqual(fields, tmp):
+ return ind
+ ind += 1
+
+ return -1
+
+
+def write_input(fc, my_chemistry, primordial_chemistry):
+ """
+ Write the simulation input in the HDF5 file
+ """
+ f = File(filename, "a")
+ data_name = "PrimordialChemistry%i" % primordial_chemistry
+ if data_name in f:
+ data = f[data_name]
+ else:
+ print("Creating Dataset")
+ data = f.create_group(data_name)
+
+ i = getSameSimulation(data, fc["density"], my_chemistry)
+ if i != -1:
+ print("Updating Simulation")
+ data.pop(data.keys()[i])
+ data = data.create_group("Simulation %i" % len(data))
+
+ # Units
+ data.create_group("Units")
+ tmp = data["Units"]
+ tmp.attrs["Length"] = my_chemistry.length_units
+ tmp.attrs["Density"] = my_chemistry.density_units
+ tmp.attrs["Velocity"] = my_chemistry.velocity_units
+ tmp.attrs["Time"] = my_chemistry.time_units
+
+ # Parameters
+ data.create_group("Params")
+ tmp = data["Params"]
+ tmp.attrs["MetalCooling"] = my_chemistry.metal_cooling
+ tmp.attrs["UVBackground"] = my_chemistry.UVbackground
+ tmp.attrs["SelfShieldingMethod"] = my_chemistry.self_shielding_method
+ tmp.attrs["MetalMassFraction"] = metal_mass_fraction
+ tmp.attrs["ScaleFactor"] = my_chemistry.a_value
+ tmp.attrs["Tolerance"] = tolerance
+ tmp.attrs["MaxIter"] = max_iter
+ tmp.attrs["Density"] = fc["density"]
+ tmp.attrs["Temperature"] = part_temperature
+ tmp.attrs["Energy"] = fc["energy"][0]
+
+ # Inputs
+ data.create_group("Input")
+ tmp = data["Input"]
+
+ write_fractions(tmp, fc, primordial_chemistry)
+
+ return f, data
+
+
+def write_fractions(tmp, fc, primordial_chemistry):
+ """
+ Write the element fractions in the HDF5 file
+ """
+ fields = get_fields(primordial_chemistry)
+ for i in fields:
+ dset = tmp.create_dataset(i, fc[i].shape,
+ dtype=fc[i].dtype)
+ dset[:] = fc[i] / fc["density"]
+
+
+def get_fields(primordial_chemistry):
+ """
+ Return the fields present in grackle
+ """
+ fields = [
+ "metal"
+ ]
+ if primordial_chemistry > 0:
+ fields.extend([
+ "HI",
+ "HII",
+ "HeI",
+ "HeII",
+ "HeIII",
+ "de"])
+
+ if primordial_chemistry > 1:
+ fields.extend([
+ "HM",
+ "H2I",
+ "H2II"
+ ])
+
+ if primordial_chemistry > 2:
+ fields.extend([
+ "DI",
+ "DII",
+ "HDI"
+ ])
+ return fields
+
+
+if __name__ == "__main__":
+ for i in range(4):
+ print("Computing Primordial Chemistry %i" % i)
+ generate_data(i)
diff --git a/examples/cooling_box/getCoolingTable.sh b/examples/cooling_box/getCoolingTable.sh
new file mode 100644
index 0000000000000000000000000000000000000000..809958bfc236e5ab0f7be924c62789fa13b3ac29
--- /dev/null
+++ b/examples/cooling_box/getCoolingTable.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5
diff --git a/examples/cooling_box/plot_cooling.py b/examples/cooling_box/plot_cooling.py
new file mode 100644
index 0000000000000000000000000000000000000000..917ea82d48a324c1a5fbceb076b89f865129a5dd
--- /dev/null
+++ b/examples/cooling_box/plot_cooling.py
@@ -0,0 +1,335 @@
+#!/usr/bin/env python3
+"""
+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/>.
+"""
+
+from pyswiftsim import wrapper
+
+import numpy as np
+import matplotlib.pyplot as plt
+from os.path import isfile
+from h5py import File
+from glob import glob
+
+plt.rc('text', usetex=True)
+
+# PARAMETERS
+
+# grackle primordial chemistry (for comparison)
+primordial_chemistry = 0
+
+# reference data
+grackle_filename = "grackle.hdf5"
+compute_equilibrium = False
+
+# swift params filename
+filename = "cooling.yml"
+
+# particle data
+# in atom/cc if generating
+# if reading => same than in grackle.hdf5
+part_density = 2.4570967948472597e7
+part_temperature = 1e6 # in K
+gamma = 5./3.
+
+# default parameters
+t_end = 1.
+N_time = 100000
+
+# cooling output
+swift_output_filename = "swift.hdf5"
+
+# simulations
+simulations = {
+ "SWIFT": "/home/loikki/data/swiftsim/examples/CoolingBox/coolingBox_*.hdf5"
+}
+# SCRIPT
+
+
+def get_fields(primordial_chemistry):
+ fields = [
+ "metal"
+ ]
+ if primordial_chemistry > 0:
+ fields.extend([
+ "HI",
+ "HII",
+ "HeI",
+ "HeII",
+ "HeIII",
+ "de"])
+
+ if primordial_chemistry > 1:
+ fields.extend([
+ "HM",
+ "H2I",
+ "H2II"
+ ])
+
+ if primordial_chemistry > 2:
+ fields.extend([
+ "DI",
+ "DII",
+ "HDI"
+ ])
+ return fields
+
+
+def generate_default_initial_condition(us, pconst):
+ print("Generating default initial conditions")
+ d = {}
+ # generate grid
+ d["temperature"] = part_temperature
+
+ # Deal with units
+ rho = part_density * us.UnitLength_in_cgs**3 * pconst.const_proton_mass
+ d["density"] = [rho]
+
+ energy = pconst.const_boltzmann_k * part_temperature \
+ / us.UnitTemperature_in_cgs
+ energy /= (gamma - 1.) * pconst.const_proton_mass
+ d["energy"] = energy
+
+ d["time"] = np.linspace(0, t_end, N_time)
+
+ return d
+
+
+def get_io_fields(cooling, chemistry, output):
+ a = 1. / (1. + cooling.redshift)
+ fields = {
+ "MetalCooling": cooling.with_metal_cooling,
+ "UVBackground": cooling.with_uv_background,
+ "SelfShieldingMethod": cooling.self_shielding_method,
+ "MetalMassFraction": chemistry.initial_metallicity,
+ "ScaleFactor": a,
+ "Temperature": part_temperature,
+ "Density": part_density,
+ }
+
+ if output:
+ fields.update({
+ "MaxIter": cooling.max_step,
+ "Tolerance": cooling.convergence_limit,
+ })
+ return fields
+
+
+def getSimulationIndice(data, cooling, chemistry, output):
+
+ fields = get_io_fields(cooling, chemistry, output)
+
+ ind = 0
+ for name in data:
+ test = True
+ d = data[name]
+ if "Params" in d:
+ tmp = d["Params"]
+ for key in fields.keys():
+ check = abs(tmp.attrs[key] - fields[key]) > \
+ 1e-3 * tmp.attrs[key]
+ if check:
+ test = False
+ break
+ if test:
+ break
+ ind += 1
+ if ind == len(data):
+ return -1
+ else:
+ return ind
+
+
+def read_grackle_data(filename, us, pconst, primordial_chemistry,
+ cooling, chemistry):
+ print("Reading initial conditions")
+ f = File(filename, "r")
+ data = f["PrimordialChemistry%i" % primordial_chemistry]
+
+ i = getSimulationIndice(data, cooling, chemistry, output=False)
+ if i < 0:
+ raise IOError("Unable to locate corresponding grackle simulation")
+ else:
+ tmp = list(data.keys())
+ data = data[tmp[i]]
+
+ # read units
+ tmp = data["Units"].attrs
+
+ u_len = tmp["Length"] / us.UnitLength_in_cgs
+ u_time = tmp["Time"] / us.UnitTime_in_cgs
+ u_den = tmp["Density"] * us.UnitLength_in_cgs**3 / us.UnitMass_in_cgs
+
+ # read input
+ tmp = data["Params"]
+ energy = tmp.attrs["Energy"] * u_len**2 / u_time**2
+ d["energy"] = energy
+
+ density = tmp.attrs["Density"] * u_den
+ d["density"] = density
+
+ # # read fractions
+ # tmp = data["Input"]
+ # for i in get_fields(primordial_chemistry):
+ # d[i] = tmp[i][:]
+
+ # read output
+ tmp = data["Output"]
+
+ energy = tmp["Energy"][:] * u_len**2 / u_time**2
+ d["out_energy"] = energy
+
+ t = tmp["Time"][:] * u_time
+ d["time"] = t
+
+ f.close()
+ return d
+
+
+def initialize_swift(filename):
+ print("Initialization of SWIFT")
+ d = {}
+
+ # parse swift params
+ params = wrapper.parserReadFile(filename)
+ d["params"] = params
+
+ # init units
+ us, pconst = wrapper.unitSystemInit(params)
+ d["us"] = us
+ d["pconst"] = pconst
+
+ cosmo = wrapper.cosmologyInit(params, us, pconst)
+ d["cosmo"] = cosmo
+
+ # init cooling
+ cooling = wrapper.coolingInit(params, us, pconst)
+ d["cooling"] = cooling
+
+ # init chemistry
+ chemistry = wrapper.chemistryInit(params, us, pconst)
+ d["chemistry"] = chemistry
+
+ return d
+
+
+def plot_solution(energy, data, us):
+ print("Plotting solution")
+ rho = data["density"]
+ time = data["time"]
+
+ ref = False
+ if "out_energy" in data:
+ ref = True
+ grackle_energy = data["out_energy"]
+
+ # change units => cgs
+ rho *= us.UnitMass_in_cgs / us.UnitLength_in_cgs**3
+
+ energy *= us.UnitLength_in_cgs**2 / us.UnitTime_in_cgs**2
+
+ if ref:
+ grackle_energy *= us.UnitLength_in_cgs**2 / us.UnitTime_in_cgs**2
+
+ time *= us.UnitTime_in_cgs
+
+ # do plot
+
+ plt.figure()
+ plt.plot(time, energy, label="PySWIFTsim, %s" % wrapper.configGetCooling())
+ if ref:
+ plt.plot(time, grackle_energy,
+ label="Grackle, Prim. Chem. %i" % primordial_chemistry)
+
+ for i in simulations:
+ files = glob(simulations[i])
+ files.sort()
+ E = np.zeros([len(files)])
+ t = np.zeros([len(files)])
+ j = 0
+ for f in files:
+ f = File(f)
+ tmp = f["Units"]
+ UL = tmp.attrs["Unit length in cgs (U_L)"]
+ UT = tmp.attrs["Unit time in cgs (U_t)"]
+ UE = UL**2 / UT**2
+ tmp = f["PartType0"]
+ E[j] = np.mean(tmp["InternalEnergy"]) * UE
+ t[j] = f["Header"].attrs["Time"] * UT
+ print(t[j])
+ j += 1
+ cooling_type = f["SubgridScheme"].attrs["Cooling Model"]
+ plt.plot(t, E, label=i + ": " + cooling_type.decode("utf8"))
+ plt.xlabel("Time [s]")
+ plt.ylabel("Energy [erg]")
+ plt.legend()
+ plt.show()
+
+
+if __name__ == "__main__":
+
+ d = initialize_swift(filename)
+ pconst = d["pconst"]
+ us = d["us"]
+ params = d["params"]
+ cosmo = d["cosmo"]
+ cooling = d["cooling"]
+ chemistry = d["chemistry"]
+
+ if isfile(grackle_filename) and False:
+ d = read_grackle_data(grackle_filename, us, pconst,
+ primordial_chemistry, cooling,
+ chemistry)
+ else:
+ d = generate_default_initial_condition(us, pconst)
+
+ 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)
+ if not compute_equilibrium:
+ energy[i+1] = wrapper.doCooling(pconst, us,
+ cosmo,
+ cooling,
+ chemistry,
+ rho,
+ ene,
+ dt)
+ else:
+ fields = get_fields(primordial_chemistry)
+ N = len(fields)
+ frac = np.zeros([1, N])
+ for j in range(N):
+ frac[:, j] = d[fields[j]]
+ energy[i+1] = wrapper.doCooling(pconst, us,
+ cosmo,
+ cooling,
+ chemistry,
+ rho,
+ ene,
+ dt,
+ frac.astype(np.float32))
+
+ 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..d2f01dae43b010839cddca2d11d20d474dd76f0f
--- /dev/null
+++ b/examples/cooling_box/run.sh
@@ -0,0 +1,17 @@
+#!/bin/bash
+
+# Get the Grackle cooling table
+if [ ! -e CloudyData_UVB=HM2012.h5 ]
+then
+ echo "Fetching the Cloudy tables required by Grackle..."
+ ./getCoolingTable.sh
+fi
+
+# generate grackle data
+if [ ! -e grackle.hdf5 ]
+then
+ echo "Generating Grackle data..."
+ ./generate_grackle_data.py
+fi
+
+./plot_cooling.py
diff --git a/examples/cooling_rate/cooling.yml b/examples/cooling_rate/cooling.yml
index 54c4cf3c2ecdfe33e16cb28ebe16243206ecffac..10608d34b8e2fcdb2ad2db441e5b0ae655195541 100644
--- a/examples/cooling_rate/cooling.yml
+++ b/examples/cooling_rate/cooling.yml
@@ -34,6 +34,17 @@ InitialConditions:
file_name: ./sedov.hdf5 # The file to read
h_scaling: 3.33
+# 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:
@@ -44,9 +55,9 @@ GrackleCooling:
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
- ConvergenceLimit: 1e-3
- MaxSteps: 100000
+ ConvergenceLimit: 1e-4
+ MaxSteps: 20000
+ Omega: 0.8
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
@@ -55,3 +66,6 @@ Gravity:
a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
r_cut_max: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
r_cut_min: 0.1 # (Optional) Cut-off in number of top-level cells below which no truncation of FMM forces are performed (this is the default value).
+
+GearChemistry:
+ InitialMetallicity: 0.01295
\ No newline at end of file
diff --git a/examples/cooling_rate/generate_grackle_data.py b/examples/cooling_rate/generate_grackle_data.py
index 3ebb364bf9d227113b781cbf068ea2c86290219f..290c6611dd5a454242de02f94d8a6efa97291cca 100644
--- a/examples/cooling_rate/generate_grackle_data.py
+++ b/examples/cooling_rate/generate_grackle_data.py
@@ -34,8 +34,18 @@ filename = "grackle.hdf5"
debug = 1
+# metal_mass_fraction = 1e-10
+metal_mass_fraction = 0.01295
+
+tolerance = 1e-5
+max_iter = 1e8
+
def generate_data(primordial_chemistry):
+ """
+ Compute the cooling rate of a bunch of particles and save
+ them to an hdf5 file
+ """
current_redshift = 0.
# Set solver parameters
@@ -66,12 +76,12 @@ def generate_data(primordial_chemistry):
# This container holds the solver parameters, units, and fields.
temperature = np.logspace(1, 9, 1000)
fc = setup_fluid_container(my_chemistry,
- metal_mass_fraction=0.01295,
+ metal_mass_fraction=metal_mass_fraction,
temperature=temperature,
density=density,
converge=True,
- tolerance=1e-5,
- max_iterations=1e8)
+ tolerance=tolerance,
+ max_iterations=max_iter)
f, data = write_input(
fc, my_chemistry, temperature, dt, primordial_chemistry)
@@ -84,6 +94,9 @@ def generate_data(primordial_chemistry):
def write_output(f, data, fc, rate):
+ """
+ Write the output of the simulation
+ """
data.create_group("Output")
tmp = data["Output"]
# Rate
@@ -98,17 +111,60 @@ def write_output(f, data, fc, rate):
f.close()
+def simulationAreEqual(fields, data):
+ """
+ Check if current simulation has already
+ been run
+ """
+ test = True
+ for key in fields.keys():
+ if data.attrs[key] != fields[key]:
+ test = False
+ break
+
+ return test
+
+
+def getSameSimulation(data, my_chemistry):
+ """
+ Get indice of the same simulation.
+ returns -1 if not found
+ """
+ fields = {
+ "MetalCooling": my_chemistry.metal_cooling,
+ "UVBackground": my_chemistry.UVbackground,
+ "SelfShieldingMethod": my_chemistry.self_shielding_method,
+ "MetalMassFraction": metal_mass_fraction,
+ "ScaleFactor": my_chemistry.a_value
+ }
+ ind = 0
+ for d in data:
+ if "Params" in data[d]:
+ tmp = data[d]["Params"]
+ if simulationAreEqual(fields, tmp):
+ return ind
+ ind += 1
+
+ return -1
+
+
def write_input(fc, my_chemistry, temperature, dt, primordial_chemistry):
+ """
+ Write the simulation input in the HDF5 file
+ """
f = File(filename, "a")
data_name = "PrimordialChemistry%i" % primordial_chemistry
if data_name in f:
- print("Updating Dataset")
data = f[data_name]
- f.pop(data_name)
else:
print("Creating Dataset")
+ data = f.create_group(data_name)
- data = f.create_group(data_name)
+ i = getSameSimulation(data, my_chemistry)
+ if i != -1:
+ print("Updating Simulation")
+ data.pop(data.keys()[i])
+ data = data.create_group("Simulation %i" % len(data))
# Units
data.create_group("Units")
@@ -124,7 +180,11 @@ def write_input(fc, my_chemistry, temperature, dt, primordial_chemistry):
tmp.attrs["MetalCooling"] = my_chemistry.metal_cooling
tmp.attrs["UVBackground"] = my_chemistry.UVbackground
tmp.attrs["SelfShieldingMethod"] = my_chemistry.self_shielding_method
+ tmp.attrs["MetalMassFraction"] = metal_mass_fraction
+ tmp.attrs["ScaleFactor"] = my_chemistry.a_value
tmp.attrs["TimeStep"] = dt
+ tmp.attrs["Tolerance"] = tolerance
+ tmp.attrs["MaxIter"] = max_iter
# Inputs
data.create_group("Input")
@@ -148,6 +208,9 @@ def write_input(fc, my_chemistry, temperature, dt, primordial_chemistry):
def write_fractions(tmp, fc, primordial_chemistry):
+ """
+ Write the element fractions in the HDF5 file
+ """
fields = get_fields(primordial_chemistry)
for i in fields:
dset = tmp.create_dataset(i, fc[i].shape,
@@ -156,6 +219,9 @@ def write_fractions(tmp, fc, primordial_chemistry):
def get_fields(primordial_chemistry):
+ """
+ Return the fields present in grackle
+ """
fields = [
"metal"
]
diff --git a/examples/cooling_rate/plot_cooling.py b/examples/cooling_rate/plot_cooling.py
index 1cf9aae0de8a26b35591506e3ab289d17933c53a..62bc5377c5de9357a543ac34faeb92c4de7ef063 100644
--- a/examples/cooling_rate/plot_cooling.py
+++ b/examples/cooling_rate/plot_cooling.py
@@ -1,24 +1,44 @@
#!/usr/bin/env python3
+"""
+This file is part of PYSWIFTSIM.
+Copyright (c) 2018 Loic Hausammann (loic.hausammann@epfl.ch)
-from pyswiftsim import wrapper
+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 libpyswiftsim import coolingRate
from copy import deepcopy
import numpy as np
import matplotlib.pyplot as plt
-from astropy import units
+from astropy import units, constants
from os.path import isfile
from h5py import File
+import yaml
plt.rc('text', usetex=True)
# PARAMETERS
+# cosmology
+with_cosmo = 0
-# grackle primordial chemistry
-primordial_chemistry = 1
+# grackle primordial chemistry (for comparison)
+primordial_chemistry = 0
# reference data
grackle_filename = "grackle.hdf5"
-compute_equilibrium = True
+compute_equilibrium = False
# swift params filename
filename = "cooling.yml"
@@ -36,7 +56,7 @@ else:
default_density = np.logspace(-3, 1, N_rho)
# temperature in K
-N_T = 10
+N_T = 1000
default_temperature = np.logspace(1, 9, N_T)
# time step in s
@@ -46,6 +66,9 @@ default_dt = default_dt.to("s") / units.s
# adiabatic index
gamma = 5. / 3.
+# cooling output
+swift_output_filename = "swift.hdf5"
+
# SCRIPT
@@ -78,7 +101,7 @@ def get_fields(primordial_chemistry):
return fields
-def generate_default_initial_condition(us, pconst):
+def generate_default_initial_condition(us):
print("Generating default initial conditions")
d = {}
# generate grid
@@ -88,38 +111,60 @@ def generate_default_initial_condition(us, pconst):
d["temperature"] = T
# Deal with units
- rho *= us.UnitLength_in_cgs**3 * pconst.const_proton_mass
+ m_p = constants.m_p.to("{} g".format(us["UnitLength_in_cgs"]))
+ rho *= us["UnitLength_in_cgs"]**3 * m_p
d["density"] = rho
- energy = pconst.const_boltzmann_k * T / us.UnitTemperature_in_cgs
- energy /= (gamma - 1.) * pconst.const_proton_mass
+ unit_time = us["UnitLength_in_cgs"] / us["UnitVelocity_in_cgs"]
+ unit_energy = us["UnitMass_in_cgs"] * us["UnitLength_in_cgs"]**2
+ unit_energy /= unit_time**2
+ k_B = constants.k_B.to("{} J / {} K".format(
+ unit_energy, us["UnitTemp_in_cgs"]
+ ))
+ energy = k_B * T / us["UnitTemp_in_cgs"]
+ energy /= (gamma - 1.) * m_p
d["energy"] = energy
- dt = default_dt / us.UnitTime_in_cgs
+ dt = default_dt / unit_time
d["time_step"] = dt
return d
-def read_grackle_data(filename, us, primordial_chemistry):
+def read_grackle_data(filename, us):
print("Reading initial conditions")
f = File(filename, "r")
data = f["PrimordialChemistry%i" % primordial_chemistry]
+<<<<<<< HEAD
+ print(cooling)
+ i = getSimulationIndice(data, cooling, chemistry, output=False)
+ if i < 0:
+ raise IOError("Unable to locate corresponding grackle simulation")
+ else:
+ tmp = list(data.keys())
+ data = data[tmp[i]]
+=======
+ data = data[grackle_simulation]
+>>>>>>> 083548de07e9dd7ccb17d962bb3e253d656c0a83
+
# read units
tmp = data["Units"].attrs
- u_len = tmp["Length"] / us.UnitLength_in_cgs
- u_den = tmp["Density"] * us.UnitLength_in_cgs**3 / us.UnitMass_in_cgs
- u_time = tmp["Time"] / us.UnitTime_in_cgs
+ u_len = tmp["Length"] / float(us["UnitLength_in_cgs"])
+ u_den = tmp["Density"] * float(us["UnitLength_in_cgs"])**3
+ u_den /= float(us["UnitMass_in_cgs"])
+ u_time = float(us["UnitLength_in_cgs"]) / float(us["UnitVelocity_in_cgs"])
+ u_time = tmp["Time"] / u_time
# read input
tmp = data["Input"]
+ d = {}
energy = tmp["Energy"][:] * u_len**2 / u_time**2
d["energy"] = energy
- T = tmp["Temperature"][:] / us.UnitTemperature_in_cgs
+ T = tmp["Temperature"][:] / float(us["UnitTemp_in_cgs"])
d["temperature"] = T
density = tmp["Density"][:] * u_den
@@ -145,23 +190,10 @@ def read_grackle_data(filename, us, primordial_chemistry):
return d
-def initialize_swift(filename):
- print("Initialization of SWIFT")
- d = {}
-
- # parse swift params
- params = wrapper.parserReadFile(filename)
- d["params"] = params
-
- # init units
- us, pconst = wrapper.unitSystemInit(params)
- d["us"] = us
- d["pconst"] = pconst
-
- # init cooling
- cooling = wrapper.coolingInit(params, us, pconst)
- d["cooling"] = cooling
- return d
+def getParser(filename):
+ with open(filename, "r") as stream:
+ stream = yaml.load(stream)
+ return stream
def plot_solution(rate, data, us):
@@ -176,16 +208,19 @@ def plot_solution(rate, data, us):
grackle_rate = data["rate"]
# change units => cgs
- rho *= us.UnitMass_in_cgs / us.UnitLength_in_cgs**3
+ u_m = float(us["UnitMass_in_cgs"])
+ u_l = float(us["UnitLength_in_cgs"])
+ u_t = float(us["UnitLength_in_cgs"]) / float(us["UnitVelocity_in_cgs"])
+ rho *= u_m / u_l**3
- T *= us.UnitTemperature_in_cgs
+ T *= float(us["UnitTemp_in_cgs"])
- energy *= us.UnitLength_in_cgs**2 / us.UnitTime_in_cgs**2
+ energy *= u_l**2 / u_t**2
- rate *= us.UnitLength_in_cgs**2 / us.UnitTime_in_cgs**3
+ rate *= u_l**2 / u_t**3
if ref:
- grackle_rate *= us.UnitLength_in_cgs**2 / us.UnitTime_in_cgs**3
+ grackle_rate *= u_l**2 / u_t**3
# lambda cooling
lam = rate * rho
@@ -198,7 +233,7 @@ def plot_solution(rate, data, us):
# plot Lambda vs T
plt.figure()
plt.loglog(T, np.abs(lam),
- label="SWIFT, %s" % wrapper.configGetCooling())
+ label="SWIFT")
if ref:
# plot reference
plt.loglog(T, np.abs(grackle_lam),
@@ -247,32 +282,38 @@ def plot_solution(rate, data, us):
tc = np.arange(_min, _max, 1.5)
cbar.set_ticks(tc)
cbar.set_ticklabels(tc)
+ cbar.ax.set_ylabel("Log( Cooling Length / kpc )")
plt.show()
if __name__ == "__main__":
- d = initialize_swift(filename)
- pconst = d["pconst"]
- us = d["us"]
- params = d["params"]
- cooling = d["cooling"]
+ parser = getParser(filename)
+ us = parser["InternalUnitSystem"]
- if isfile(grackle_filename):
- d = read_grackle_data(grackle_filename, us, primordial_chemistry)
+ if isfile(grackle_filename) and N_rho == 1:
+ d = read_grackle_data(grackle_filename, us)
else:
- d = generate_default_initial_condition(us, pconst)
+ d = generate_default_initial_condition(us)
# du / dt
print("Computing cooling...")
rate = np.zeros(d["density"].shape)
if compute_equilibrium:
- rate = wrapper.coolingRate(pconst, us, cooling,
- d["density"].astype(np.float32),
- d["energy"].astype(np.float32),
- d["time_step"])
+<<<<<<< HEAD
+ rate = wrapper.coolingRate(
+ filename,
+ d["density"].astype(np.float32),
+ d["energy"].astype(np.float32),
+ d["time_step"])
+=======
+ rate = coolingRate(filename, with_cosmo,
+ d["density"].astype(np.float32),
+ d["energy"].astype(np.float32),
+ d["time_step"])
+>>>>>>> 083548de07e9dd7ccb17d962bb3e253d656c0a83
else:
fields = get_fields(primordial_chemistry)
N = len(fields)
@@ -280,10 +321,10 @@ if __name__ == "__main__":
for i in range(N):
frac[:, i] = d[fields[i]]
- rate = wrapper.coolingRate(pconst, us, cooling,
- d["density"].astype(np.float32),
- d["energy"].astype(np.float32),
- d["time_step"],
- frac.astype(np.float32))
+ rate = coolingRate(filename, with_cosmo,
+ d["density"].astype(np.float32),
+ d["energy"].astype(np.float32),
+ d["time_step"],
+ frac.astype(np.float32))
plot_solution(rate, d, us)
diff --git a/pyswiftsim/__init__.py b/pyswiftsim/__init__.py
index 50a5cb6cd4d2eec05cc57e8c426c001cd68b2db3..526dbce334dd759ae3af7299c78404f02c44fa8d 100644
--- a/pyswiftsim/__init__.py
+++ b/pyswiftsim/__init__.py
@@ -1,2 +1,19 @@
+"""
+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 *
diff --git a/pyswiftsim/dev.py b/pyswiftsim/dev.py
index 6646ecf4235022217387f0775c6fdf17ad3dd730..64561b3bcb36adc66df701352e5705d95afead36 100644
--- a/pyswiftsim/dev.py
+++ b/pyswiftsim/dev.py
@@ -1,4 +1,23 @@
#!/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 sys import argv
symbol_table = {
" int ": "i",
@@ -6,6 +25,7 @@ symbol_table = {
" char ": "c",
}
+
def _generateFormatString(filename):
"""
Open a file containing a single struct and generate a format string
@@ -14,7 +34,7 @@ def _generateFormatString(filename):
struct_format = ""
struct_attribute = []
-
+
with open(filename, "r") as f:
line = ""
# skip until reaching struct
@@ -46,21 +66,20 @@ def _generateFormatString(filename):
j = line.index(";")
struct_attribute.append(line[i:j])
-
return struct_format, struct_attribute
-
+
if __name__ == "__main__":
- filename = "chemistry_data.h"
+ filename = argv[-1]
struct_format, struct_attribute = _generateFormatString(filename)
attr = "[\n"
for i in struct_attribute:
attr += "\t'%s',\n" % i
attr += "]"
-
+
txt = """
_format = "{form}"
_name = {attr}
diff --git a/pyswiftsim/structure.py b/pyswiftsim/structure.py
index 78a35b781e6062e2f133bbcdccb5d726ae67b440..184dff7107432e44cf71836d197c732db18a9ef0 100644
--- a/pyswiftsim/structure.py
+++ b/pyswiftsim/structure.py
@@ -1,3 +1,20 @@
+"""
+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 wrapper
import struct
@@ -286,7 +303,7 @@ class UnitSystem(SwiftStruct):
class ChemistryPartData(SwiftStruct):
_format = "f"
_name = [
- "he_density"
+ "initial_metallicity"
]
def __init__(self, data, parent=None):
@@ -330,14 +347,16 @@ class Part(SwiftStruct):
# #
######################################################################
class Parameter(SwiftStruct):
- _format = "{line_size}c{line_size}c".format(
+ _format = "{line_size}c{line_size}cii".format(
line_size=PARSER_MAX_LINE_SIZE
)
_name = [
- "name",
- "value"
- ]
+ "name",
+ "value",
+ "used",
+ "is_default"
+ ]
def __init__(self, data, parent=None):
super().__init__(self.struct_format, data, parent)
@@ -406,14 +425,16 @@ class SwiftParams(SwiftStruct):
# #
######################################################################
class PhysConst(SwiftStruct):
- _format = "dddddddddddddddd"
+ _format = "ddddddddddddddddddddd"
_name = [
"const_newton_G",
"const_speed_light_c",
"const_planck_h",
"const_planck_hbar",
"const_boltzmann_k",
+ "const_avogadro_number",
"const_thomson_cross_section",
+ "const_stefan_boltzmann",
"const_electron_charge",
"const_electron_volt",
"const_electron_mass",
@@ -424,79 +445,9 @@ class PhysConst(SwiftStruct):
"const_light_year",
"const_solar_mass",
"const_earth_mass",
- ]
-
- def __init__(self, data, parent=None):
- super().__init__(self.struct_format, data, parent)
-
-
-class GrackleCodeUnits(SwiftStruct):
- cooling_type = wrapper.configGetCooling()
- _format = "idddddd"
- if cooling_type == "grackle_float":
- _format = _format.replace("d", "f")
- _name = [
- "comoving_coordinates",
- "density_units",
- "length_units",
- "time_units",
- "velocity_units",
- "a_units",
- "a_value"
- ]
-
- def __init__(self, data, parent=None):
- super().__init__(self.struct_format, data, parent)
-
-
-class GrackleChemistryData(SwiftStruct):
- cooling_type = wrapper.configGetCooling()
- _format = "iiiiiPidiidiiiiiiidddiiddiddiiddddddiiiiii"
- if cooling_type == "grackle_float":
- _format = _format.replace("d", "f")
- _name = [
- 'use_grackle',
- 'with_radiative_cooling',
- 'primordial_chemistry',
- 'metal_cooling',
- 'UVbackground',
- 'grackle_data_file',
- 'cmb_temperature_floor',
- 'Gamma',
- 'h2_on_dust',
- 'photoelectric_heating',
- 'photoelectric_heating_rate',
- 'use_volumetric_heating_rate',
- 'use_specific_heating_rate',
- 'three_body_rate',
- 'cie_cooling',
- 'h2_optical_depth_approximation',
- 'ih2co',
- 'ipiht',
- 'HydrogenFractionByMass',
- 'DeuteriumToHydrogenRatio',
- 'SolarMetalFractionByMass',
- 'NumberOfTemperatureBins',
- 'CaseBRecombination',
- 'TemperatureStart',
- 'TemperatureEnd',
- 'NumberOfDustTemperatureBins',
- 'DustTemperatureStart',
- 'DustTemperatureEnd',
- 'Compton_xray_heating',
- 'LWbackground_sawtooth_suppression',
- 'LWbackground_intensity',
- 'UVbackground_redshift_on',
- 'UVbackground_redshift_off',
- 'UVbackground_redshift_fullon',
- 'UVbackground_redshift_drop',
- 'cloudy_electron_fraction_factor',
- 'use_radiative_transfer',
- 'radiative_transfer_coupled_rate_solver',
- 'radiative_transfer_intermediate_step',
- 'radiative_transfer_hydrogen_only',
- 'self_shielding_method',
- 'H2_self_shielding',
+ "const_T_CMB_0",
+ "const_primordial_He_fraction",
+ "const_reduced_hubble"
]
def __init__(self, data, parent=None):
@@ -515,38 +466,98 @@ class CoolingFunctionData(SwiftStruct):
"cooling_tstep_mult"
]
elif "grackle" in cooling_type:
- _format = "200cidd{code_units}s{chemistry}s".format(
- code_units=struct.calcsize(GrackleCodeUnits._format),
- chemistry=struct.calcsize(GrackleChemistryData._format)
+ _format = "200cid{code_units}c{chemistry}ciiiifif".format(
+ code_units=56,
+ chemistry=248
)
_name = [
- "GrackleCloudyTable",
- "UVbackground",
- "GrackleRedshift",
- "GrackleHSShieldingDensityThreshold",
+ "cloudy_table",
+ "with_uv_background",
+ "redshift",
"code_units",
"chemistry",
+ "with_metal_cooling",
+ "provide_volumetric_heating_rate",
+ "provide_specificy_heating_rate",
+ "self_shielding_method",
+ "convergence_limit",
+ "max_step",
+ "omega"
]
- @property
- def struct_substruct(self):
- chemistry = {
- "class": GrackleChemistryData,
- "size": 1
- }
-
- units = {
- "class": GrackleCodeUnits,
- "size": 1
- }
- return {
- "units": units,
- "chemistry": chemistry
- }
-
else:
raise ValueError(
"Cooling Type %s not implemented" % cooling_type)
def __init__(self, data, parent=None):
super().__init__(self.struct_format, data, parent)
+
+
+class ChemistryGlobalData(SwiftStruct):
+ _format = "f"
+ _name = [
+ "initial_metallicity",
+ ]
+
+ def __init__(self, data, parent=None):
+ super().__init__(self.struct_format, data, parent)
+
+######################################################################
+# #
+# PhysConst #
+# #
+######################################################################
+
+
+class Cosmology(SwiftStruct):
+ _format = "ddddddddddddddddddddddddddddddddddddppppdd"
+ _name = [
+ 'a',
+ 'a_inv',
+ 'a2_inv',
+ 'a3_inv',
+ 'a_factor_internal_energy',
+ 'a_factor_pressure',
+ 'a_factor_sound_speed',
+ 'a_factor_mu',
+ 'a_factor_Balsara_eps',
+ 'a_factor_grav_accel',
+ 'a_factor_hydro_accel',
+ 'z',
+ 'H',
+ 'critical_density',
+ 'critical_density_0',
+ 'time_step_factor',
+ 'a_dot',
+ 'time',
+ 'lookback_time',
+ 'w',
+ 'a_begin',
+ 'a_end',
+ 'time_begin',
+ 'time_end',
+ 'time_base',
+ 'time_base_inv',
+ 'h',
+ 'H0',
+ 'Hubble_time',
+ 'Omega_m',
+ 'Omega_b',
+ 'Omega_lambda',
+ 'Omega_r',
+ 'Omega_k',
+ 'w_0',
+ 'w_a',
+ 'log_a_begin',
+ 'log_a_end',
+ 'drift_fac_interp_table',
+ 'grav_kick_fac_interp_table',
+ 'hydro_kick_fac_interp_table',
+ 'hydro_kick_corr_interp_table',
+ 'time_interp_table',
+ 'time_interp_table_offset',
+ 'universe_age_at_present_day',
+ ]
+
+ def __init__(self, data, parent=None):
+ super().__init__(self.struct_format, data, parent)
diff --git a/setup.py b/setup.py
index 6f8712d727e90d1af93e599ed618d67ad68da860..919adee2849f4fc5dcc55e8dd88d25a6b2ccf337 100644
--- a/setup.py
+++ b/setup.py
@@ -17,13 +17,15 @@ cflags = [
"-Werror",
"-Wall",
"-Wextra",
+ "-Wshadow",
# disables some warnings due to python
"-Wno-unused-parameter",
+ "-Wno-maybe-uninitialized",
"-Wno-strict-prototypes",
"-Wno-unused-function",
"-Wno-incompatible-pointer-types",
"-Wno-missing-field-initializers",
- "-fopenmp"
+ "-fopenmp",
]
lflags = [
@@ -82,11 +84,14 @@ if swift_path:
])
# hdf5
+ if "bin" in hdf5_root:
+ hdf5_root = hdf5_root.split("bin")[0]
include.append(hdf5_root + "/include")
# grackle
- grackle_inc = getValueFromMakefile(swift_path, "GRACKLE_INCS = -I")
- include.append(grackle_inc)
+ if parseCmdLine("--with-grackle", store=True):
+ grackle_inc = getValueFromMakefile(swift_path, "GRACKLE_INCS = -I")
+ include.append(grackle_inc)
# C libraries
lib = [
diff --git a/src/Makefile.am b/src/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..a77473fd9cf45fb642c0175a3d22490b9b7f9084
--- /dev/null
+++ b/src/Makefile.am
@@ -0,0 +1,68 @@
+# 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/chemistry_wrapper.c b/src/chemistry_wrapper.c
new file mode 100644
index 0000000000000000000000000000000000000000..9c6cdc09f6f3919de4d69aa5ab1bab516808ed79
--- /dev/null
+++ b/src/chemistry_wrapper.c
@@ -0,0 +1,63 @@
+/*******************************************************************************
+ * 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 "pyswiftsim_tools.h"
+#include "chemistry_wrapper.h"
+
+/**
+ * @brief Initialize the chemistry
+ *
+ * args is expecting pyswiftsim classes in the following order:
+ * SwiftParams, UnitSystem and PhysConst.
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return ChemistryGlobalData
+ */
+PyObject* pychemistry_init(PyObject* self, PyObject* args) {
+ PyObject *pyparams;
+ PyObject *pyus;
+ PyObject *pypconst;
+
+ /* parse arguments */
+ if (!PyArg_ParseTuple(args, "OOO", &pyparams, &pyus, &pypconst))
+ return NULL;
+
+ struct swift_params *params = (struct swift_params*) pytools_construct(pyparams, class_swift_params);
+ if (params == NULL)
+ return NULL;
+
+ struct unit_system *us = (struct unit_system*) pytools_construct(pyus, class_unit_system);
+ if (us == NULL)
+ return NULL;
+
+ struct phys_const *pconst = (struct phys_const*) pytools_construct(pypconst, class_phys_const);
+ if (pconst == NULL)
+ return NULL;
+
+ struct chemistry_global_data chemistry;
+
+ /* init cooling */
+ chemistry_init_backend(params, us, pconst, &chemistry);
+
+ /* construct python object */
+ PyObject *pychemistry = pytools_return(&chemistry, class_chemistry_global_data);
+
+ return pychemistry;
+}
diff --git a/src/chemistry_wrapper.h b/src/chemistry_wrapper.h
new file mode 100644
index 0000000000000000000000000000000000000000..31ab1593f820845104c044fb64a312c767bcdd6b
--- /dev/null
+++ b/src/chemistry_wrapper.h
@@ -0,0 +1,27 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
+#ifndef __PYSWIFTSIM_CHEMISTRY_H__
+#define __PYSWIFTSIM_CHEMISTRY_H__
+
+#include "pyswiftsim_tools.h"
+
+PyObject* pychemistry_init(PyObject* self, PyObject* args);
+
+#endif // __PYSWIFTSIM_CHEMISTRY_H__
+
diff --git a/src/config_wrapper.h b/src/config_wrapper.h
deleted file mode 100644
index 5a17aeb3c933dc9b36a7884bc11125c7a2b5c9ef..0000000000000000000000000000000000000000
--- a/src/config_wrapper.h
+++ /dev/null
@@ -1,36 +0,0 @@
-#ifndef __PYSWIFTSIM_CONFIG_H__
-#define __PYSWIFTSIM_CONFIG_H__
-
-#include "pyswiftsim_tools.h"
-
-/**
- * @brief give the cooling type name
- *
- * @return Cooling name (PyUnicode)
- */
-PyObject* config_get_cooling() {
- char *cooling_name;
- /* lambda */
-#ifdef COOLING_CONST_LAMBDA
- cooling_name = "const_lambda";
-
- /* grackle */
-#elif defined(COOLING_GRACKLE)
-#if COOLING_GRACKLE_MODE == 0
- cooling_name = "grackle";
-#elif COOLING_GRACKLE_MODE == 1
- cooling_name = "grackle1";
-#elif COOLING_GRACKLE_MODE == 2
- cooling_name = "grackle2";
-#elif COOLING_GRACKLE_MODE == 3
- cooling_name = "grackle3";
-#else
- error("Grackle mode unknown");
-#endif // COOLING_GRACKLE_MODE
-#endif // COOLING_GRACKLE
-
- return PyUnicode_FromString(cooling_name);
-};
-
-
-#endif // __PYSWIFTSIM_CONFIG_H__
diff --git a/src/cooling_wrapper.c b/src/cooling_wrapper.c
index 11f5bdb97a6d694e8f3c396553034d31bc723015..8def5d2d1e23507a0edd4f438816df7a02aff522 100644
--- a/src/cooling_wrapper.c
+++ b/src/cooling_wrapper.c
@@ -1,50 +1,27 @@
-#include "pyswiftsim_tools.h"
-#include "cooling_wrapper.h"
-#include <omp.h>
-
-
-/**
- * @brief Initialize the cooling
+/*******************************************************************************
+ * This file is part of PYSWIFTSIM.
+ * Copyright (c) 2018 loic hausammann (loic.hausammann@epfl.ch)
*
- * args is expecting pyswiftsim classes in the following order:
- * SwiftParams, UnitSystem and PhysConst.
+ * 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.
*
- * @param self calling object
- * @param args arguments
- * @return CoolingFunctionData
- */
-PyObject* pycooling_init(PyObject* self, PyObject* args) {
- PyObject *pyparams;
- PyObject *pyus;
- PyObject *pypconst;
-
- /* parse arguments */
- if (!PyArg_ParseTuple(args, "OOO", &pyparams, &pyus, &pypconst))
- return NULL;
-
- struct swift_params *params = (struct swift_params*) pytools_construct(pyparams, class_swift_params);
- if (params == NULL)
- return NULL;
-
- struct unit_system *us = (struct unit_system*) pytools_construct(pyus, class_unit_system);
- if (us == NULL)
- return NULL;
-
- struct phys_const *pconst = (struct phys_const*) pytools_construct(pypconst, class_phys_const);
- if (pconst == NULL)
- return NULL;
-
-
- struct cooling_function_data cooling;
+ * 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"
- /* init cooling */
- cooling_init_backend(params, us, pconst, &cooling);
+/* Local includes */
+#include "pyswiftsim_tools.h"
+#include "part_wrapper.h"
- /* construct python object */
- PyObject *pycooling = pytools_return(&cooling, class_cooling_function_data);
-
- return pycooling;
-}
/**
* @brief Set the cooling element fractions
@@ -54,6 +31,7 @@ PyObject* pycooling_init(PyObject* self, PyObject* args) {
* @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);
@@ -78,6 +56,9 @@ void pycooling_set_fractions(struct xpart *xp, PyArrayObject* frac, const int id
#endif // COOLING_GRACKLE_MODE
#endif // COOLING_GRACKLE
+#else
+ message("Not implemented");
+#endif
}
/**
@@ -92,41 +73,34 @@ void pycooling_set_fractions(struct xpart *xp, PyArrayObject* frac, const int id
* @param args arguments
* @return cooling rate
*/
-PyArrayObject* pycooling_rate(PyObject* self, PyObject* args) {
+PyObject* pycooling_rate(PyObject* self, PyObject* args) {
import_array();
- PyObject *pycooling;
- PyObject *pyus;
- PyObject *pypconst;
+
+ char *filename;
PyArrayObject *rho;
PyArrayObject *energy;
PyArrayObject *fractions = NULL;
float dt = 1e-3;
+ int with_cosmo = 0;
/* parse argument */
- if (!PyArg_ParseTuple(args,
- "OOOOO|fOO",
- &pypconst,
- &pyus,
- &pycooling,
- &rho,
- &energy,
- &dt,
- &fractions
- ))
+ if (!PyArg_ParseTuple(
+ args, "siOO|fO", &filename, &with_cosmo,
+ &rho, &energy, &dt, &fractions))
return NULL;
/* check numpy array */
- if (pytools_check_array(energy, 1, NPY_FLOAT) != SUCCESS)
+ if (pytools_check_array(energy, 1, NPY_FLOAT, "Energy") != SUCCESS)
return NULL;
- if (pytools_check_array(rho, 1, NPY_FLOAT) != SUCCESS)
+ if (pytools_check_array(rho, 1, NPY_FLOAT, "Density") != SUCCESS)
return NULL;
if (fractions != NULL &&
- pytools_check_array(fractions, 2, NPY_FLOAT) != SUCCESS)
+ pytools_check_array(fractions, 2, NPY_FLOAT, "Fractions") != SUCCESS)
return NULL;
if (PyArray_DIM(energy, 0) != PyArray_DIM(rho, 0))
@@ -138,36 +112,42 @@ PyArrayObject* pycooling_rate(PyObject* self, PyObject* args) {
size_t N = PyArray_DIM(energy, 0);
- /* transform PyObject to C struct */
- struct cooling_function_data *cooling = (struct cooling_function_data*) pytools_construct(pycooling, class_cooling_function_data);
- if (cooling == NULL)
- return NULL;
+ /* Initialize everything */
+ struct swift_params params;
+ parser_read_file(filename, ¶ms);
- struct unit_system *us = (struct unit_system*) pytools_construct(pyus, class_unit_system);
- if (us == NULL)
- return NULL;
+ struct unit_system us;
+ units_init_from_params(&us, ¶ms, "InternalUnitSystem");
- struct phys_const *pconst = (struct phys_const*) pytools_construct(pypconst, class_phys_const);
- if (pconst == NULL)
- return NULL;
+ struct phys_const pconst;
+ phys_const_init(&us, ¶ms, &pconst);
+
+ struct cooling_function_data cooling;
+
+ /* init cooling */
+ cooling_init_backend(¶ms, &us, &pconst, &cooling);
+
+ struct chemistry_global_data chemistry;
+ chemistry_init_backend(¶ms, &us, &pconst, &chemistry);
+
+ struct cosmology cosmo;
+
+ if (with_cosmo)
+ cosmology_init(¶ms, &us, &pconst, &cosmo);
+ else
+ cosmology_init_no_cosmo(&cosmo);
struct part p;
struct xpart xp;
-#ifdef COOLING_GRACKLE
- /* set grackle_data */
- grackle_data = &cooling->chemistry;
-#endif
-
+ pyswift_part_set_to_zero(&p, &xp);
+ hydro_first_init_part(&p, &xp);
+
/* return object */
- PyArrayObject *rate = PyArray_SimpleNew(PyArray_NDIM(energy), PyArray_DIMS(energy), NPY_FLOAT);
+ PyArrayObject *rate = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(energy), PyArray_DIMS(energy), NPY_FLOAT);
- /* Release GIL */
- Py_BEGIN_ALLOW_THREADS;
-
/* loop over all particles */
-#pragma omp for
for(size_t i = 0; i < N; i++)
{
/* set particle data */
@@ -178,20 +158,22 @@ PyArrayObject* pycooling_rate(PyObject* self, PyObject* args) {
if (fractions != NULL)
pycooling_set_fractions(&xp, fractions, i);
else
- cooling_first_init_part(&p, &xp, cooling);
+ cooling_first_init_part(&pconst, &us, &cosmo, &cooling, &p, &xp);
+
+ chemistry_first_init_part(&pconst, &us, &cosmo, &chemistry, &p, &xp);
/* compute cooling rate */
- float *tmp = PyArray_GETPTR1(rate, i);
#ifdef COOLING_GRACKLE
- *tmp = cooling_rate(pconst, us, cooling, &p, &xp, dt);
+ float *tmp = PyArray_GETPTR1(rate, i);
+ *tmp = cooling_new_energy(&pconst, &us, &cosmo, &cooling, &p, &xp, dt);
+ *tmp = *tmp - hydro_get_physical_internal_energy(&p, &xp, &cosmo);
+ *tmp /= dt;
#else
- *tmp = cooling_rate(pconst, us, cooling, &p, &xp);
+ message("Not implemented");
+ //*tmp = cooling_rate(pconst, us, cosmo, cooling, &p, &xp);
#endif
}
- /* Acquire GIL */
- Py_END_ALLOW_THREADS;
-
- return rate;
+ return (PyObject *) rate;
}
diff --git a/src/cooling_wrapper.h b/src/cooling_wrapper.h
index 95717e8ac5d2867b25ddcebe33548fe86ecfa6cc..a9eeeb4d6a5061637795536f76af71cc439c6747 100644
--- a/src/cooling_wrapper.h
+++ b/src/cooling_wrapper.h
@@ -1,11 +1,27 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
#ifndef __PYSWIFTSIM_COOLING_H__
#define __PYSWIFTSIM_COOLING_H__
#include "pyswiftsim_tools.h"
-PyObject* pycooling_init(PyObject* self, PyObject* args);
-
-PyArrayObject* pycooling_rate(PyObject* self, PyObject* args);
+PyObject* pycooling_rate(PyObject* self, PyObject* args);
#endif // __PYSWIFTSIM_COOLING_H__
diff --git a/src/cosmology_wrapper.c b/src/cosmology_wrapper.c
new file mode 100644
index 0000000000000000000000000000000000000000..d9789187106f40ef505d72e3d1a2456ebbb59c63
--- /dev/null
+++ b/src/cosmology_wrapper.c
@@ -0,0 +1,67 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * 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 "pyswiftsim_tools.h"
+#include "cosmology_wrapper.h"
+
+/**
+ * @brief Initialize the chemistry
+ *
+ * args is expecting pyswiftsim classes in the following order:
+ * SwiftParams, UnitSystem and PhysConst.
+ *
+ * @param self calling object
+ * @param args arguments
+ * @return Cosmology
+ */
+PyObject* pycosmology_init(PyObject* self, PyObject* args) {
+ PyObject *pyparams;
+ PyObject *pyus;
+ PyObject *pypconst;
+ int with_cosmo = 0;
+
+ /* parse arguments */
+ if (!PyArg_ParseTuple(args, "OOO|i", &pyparams, &pyus, &pypconst, &with_cosmo))
+ return NULL;
+
+ struct swift_params *params = (struct swift_params*) pytools_construct(pyparams, class_swift_params);
+ if (params == NULL)
+ return NULL;
+
+ struct unit_system *us = (struct unit_system*) pytools_construct(pyus, class_unit_system);
+ if (us == NULL)
+ return NULL;
+
+ struct phys_const *pconst = (struct phys_const*) pytools_construct(pypconst, class_phys_const);
+ if (pconst == NULL)
+ return NULL;
+
+ struct cosmology cosmo;
+
+ /* init cooling */
+ if (with_cosmo)
+ cosmology_init(params, us, pconst, &cosmo);
+ else
+ cosmology_init_no_cosmo(&cosmo);
+
+ /* construct python object */
+ PyObject *pycosmo = pytools_return(&cosmo, class_cosmology);
+
+ return pycosmo;
+}
diff --git a/src/cosmology_wrapper.h b/src/cosmology_wrapper.h
new file mode 100644
index 0000000000000000000000000000000000000000..8a674060e1c79f01090f441ac2ad7e8cdc1ff340
--- /dev/null
+++ b/src/cosmology_wrapper.h
@@ -0,0 +1,27 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
+#ifndef __PYSWIFTSIM_COSMOLOGY_H__
+#define __PYSWIFTSIM_COSMOLOGY_H__
+
+#include "pyswiftsim_tools.h"
+
+PyObject* pycosmology_init(PyObject* self, PyObject* args);
+
+#endif // __PYSWIFTSIM_COSMOLOGY_H__
+
diff --git a/src/parser_wrapper.c b/src/parser_wrapper.c
index 3bf033f208428b2480f796ccdf415f7a09f8a2b7..9f4920fa55529604c19f02c8c0d591af9cb50bac 100644
--- a/src/parser_wrapper.c
+++ b/src/parser_wrapper.c
@@ -1,3 +1,21 @@
+/*******************************************************************************
+ * 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 "parser_wrapper.h"
#include "pyswiftsim_tools.h"
diff --git a/src/parser_wrapper.h b/src/parser_wrapper.h
index 529d352a7698f9180688033ca913efd870ddd272..223e6663d1e042857d2552bc7e5fe27ee2014b0f 100644
--- a/src/parser_wrapper.h
+++ b/src/parser_wrapper.h
@@ -1,3 +1,21 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
#ifndef __PYSWIFTSIM_PARSER_H__
#define __PYSWIFTSIM_PARSER_H__
diff --git a/src/part_wrapper.c b/src/part_wrapper.c
index 76d4512ae235c880ecd52d4c993f58f7e8ac1dbd..d8cefd848401a1e5eedb5839e5f12188e2c7738c 100644
--- a/src/part_wrapper.c
+++ b/src/part_wrapper.c
@@ -1,37 +1,30 @@
+/*******************************************************************************
+ * 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 "pyswiftsim_tools.h"
#include <Python.h>
#include <stdlib.h>
#include <string.h>
-PyObject* pypart_test_struct(PyObject *self, PyObject *args)
-{
-
- size_t N = sizeof(struct part);
-
- /* initialize particle */
- struct part *p = malloc(N);
- p->id = 0;
- for(size_t k = 0; k < DIM; k++)
- {
- p->x[k] = 1. + k;
- p->v[k] = 4. + k;
- p->a_hydro[k] = 7. + k;
- }
-
- p->h = 10;
- p->mass = 11.;
- p->rho = 12.;
- p->entropy = 13;
- p->entropy_dt = 14;
-
- hydro_init_part(p, NULL);
-
- /* create python object */
- PyObject *object = pytools_return(p, class_part);
-
- free(p);
-
- return object;
+void pyswift_part_set_to_zero(struct part *p, struct xpart *xp) {
+ p->v[0] = 0;
+ p->v[1] = 0;
+ p->v[2] = 0;
+ p->entropy = 0;
}
-
diff --git a/src/part_wrapper.h b/src/part_wrapper.h
index 5875f59da06cc9aba5ca9605a3ac4b377cc1a673..bdc3d52b5807871a4c0e656150686032917f17d8 100644
--- a/src/part_wrapper.h
+++ b/src/part_wrapper.h
@@ -1,17 +1,27 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
#ifndef __PYSWIFTSIM_PART_H__
#define __PYSWIFTSIM_PART_H__
#include <Python.h>
+#include "pyswiftsim_tools.h"
-/**
- * @brief Test function for the struct
- *
- * args is expecting no argument.
- *
- * @param self calling object
- * @param args arguments
- * @return Part
- */
-PyObject* pypart_test_struct(PyObject *self, PyObject *args);
-
+void pyswift_part_set_to_zero(struct part *p, struct xpart *xp);
+
#endif // __PYSWIFTSIM_PART_H__
diff --git a/src/pyswiftsim.h b/src/pyswiftsim.h
new file mode 100644
index 0000000000000000000000000000000000000000..f215e81490bd56169a241e862d8ad6056c42cf38
--- /dev/null
+++ b/src/pyswiftsim.h
@@ -0,0 +1,32 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
+#ifndef __PYSWIFTSIM_H__
+#define __PYSWIFTSIM_H__
+
+#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
+
+#include <Python.h>
+#include <numpy/arrayobject.h>
+
+/* Assert is already defined in Python.h */
+#undef assert
+#include "swift.h"
+
+
+#endif // __PYSWIFTSIM_H__
diff --git a/src/pyswiftsim_tools.c b/src/pyswiftsim_tools.c
index b8673de8a40465ef1806f847c9381b04bbd8dcff..5b872ec5a283030da0420dbb0bf05a38b43e774d 100644
--- a/src/pyswiftsim_tools.c
+++ b/src/pyswiftsim_tools.c
@@ -1,3 +1,21 @@
+/*******************************************************************************
+ * 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 "pyswiftsim_tools.h"
@@ -6,7 +24,9 @@ const size_t class_size[class_count] = {
sizeof(struct part),
sizeof(struct swift_params),
sizeof(struct phys_const),
- sizeof(struct cooling_function_data)
+ sizeof(struct cooling_function_data),
+ sizeof(struct chemistry_global_data),
+ sizeof(struct cosmology)
};
char *class_name[class_count] = {
@@ -14,22 +34,28 @@ char *class_name[class_count] = {
"Part",
"SwiftParams",
"PhysConst",
- "CoolingFunctionData"
+ "CoolingFunctionData",
+ "ChemistryGlobalData",
+ "Cosmology"
};
+/**
+ * @brief Import a python library/function
+ *
+ * @param module Module name
+ * @param object_name object to import from module
+ * @return (PyObject) imported object
+ */
PyObject* pytools_import(char* module_name, char* object_name)
{
/* load module */
PyObject *module;
module = PyImport_ImportModule(module_name);
-
if (module == NULL)
- {
- pyerror("Failed to import module '%s'.", module_name);
- }
-
+ pyerror("Failed to import module '%s'.", module_name);
+
/* get module dictionary */
PyObject *dict;
@@ -37,21 +63,25 @@ PyObject* pytools_import(char* module_name, char* object_name)
Py_DECREF(module);
if (dict == NULL)
- {
- pyerror("Failed to get module '%s' dictionary", module_name);
- }
+ pyerror("Failed to get module '%s' dictionary", module_name);
/* get right object */
PyObject *python_obj = PyDict_GetItemString(dict, object_name);
- Py_DECREF(dict);
-
if (python_obj == NULL)
pyerror("Object %s does not exist in module %s", object_name, module_name);
+ Py_INCREF(python_obj);
return python_obj;
}
+/**
+ * @brief Construct a python object from a C struct
+ *
+ * @param p pointer to the C struct
+ * @param class #swift_class of the pointer
+ * @return PyObject of the required class
+ */
PyObject* pytools_return(void *p, int class)
{
@@ -70,14 +100,14 @@ PyObject* pytools_return(void *p, int class)
/* import python class */
python_class = pytools_import(module_name, class_pyname);
-
- if (python_class == NULL)
+ if (!python_class)
return NULL;
if (!PyCallable_Check(python_class))
{
Py_DECREF(python_class);
- pyerror("Unable to import class %s from %s", class_pyname, module_name);
+ message("Redo");
+ pyerror("Unable to import class %30s from %.30s", class_pyname, module_name);
}
/* create object */
@@ -94,21 +124,26 @@ PyObject* pytools_return(void *p, int class)
}
+/**
+ * @brief get type name in C
+ *
+ * @param obj PyObject from which to get the name0
+ * @return type name
+ */
char* pytools_get_type_name(PyObject *obj)
{
+ /* check input */
+ if (!obj)
+ pyerror("Requires a non null object");
+
/* get object type */
- PyObject *type = PyObject_Type(obj);
+ PyObject *type = (PyObject *)obj->ob_type;
if (type == NULL)
- {
- Py_DECREF(type);
pyerror("Unable to get type");
- }
/* get object name */
PyObject* recv = PyObject_Str(type);
- Py_DECREF(type);
-
- if (recv == NULL)
+ if (!recv)
{
Py_DECREF(recv);
pyerror("Unable to get string representation");
@@ -119,15 +154,20 @@ char* pytools_get_type_name(PyObject *obj)
char *name = PyUnicode_AsUTF8AndSize(recv, &size);
Py_DECREF(recv);
- if (name == NULL)
- {
+ if (!name)
pyerror("Unable to convert string to char");
- }
return name;
}
+/**
+ * @brief Construct a C struct from a python object
+ *
+ * @param obj pointer to the python object
+ * @param class #swift_class of the pointer
+ * @return pointer to the required struct
+ */
char* pytools_construct(PyObject* obj, int class)
{
char *module_name = "pyswiftsim.structure";
@@ -142,6 +182,8 @@ char* pytools_construct(PyObject* obj, int class)
/* import class */
PyObject *pyclass = pytools_import(module_name, class_pyname);
+ if (!pyclass)
+ return NULL;
/* check if classes correspond */
int test = PyObject_IsInstance(obj, pyclass);
@@ -149,16 +191,15 @@ char* pytools_construct(PyObject* obj, int class)
if (!test)
{
char *recv = pytools_get_type_name(obj);
- if (recv == NULL)
+ if (!recv)
return NULL;
pyerror("Expecting class %s, received %s", class_pyname, recv);
}
-
/* copy python class' data to C */
PyObject* data = PyObject_GetAttrString(obj, "data");
- if (data == NULL)
+ if (!data)
pyerror("Unable to get the attribute 'data'");
char *ret = PyBytes_AsString(data);
@@ -168,26 +209,35 @@ char* pytools_construct(PyObject* obj, int class)
}
-int pytools_check_array(PyArrayObject *obj, int dim, int type)
+/**
+ * @brief Check if object is the expected PyArray type
+ *
+ * @param obj PyArray to check
+ * @param dim required dimension
+ * @param int required type
+ * @param name Name to print for errors
+ * @return #error_code result of the test
+ */
+int pytools_check_array(PyArrayObject *obj, int dim, int type, const char* name)
{
IMPORT_ARRAY1(1);
/* check if array */
if (!PyArray_Check(obj))
{
- pyerror("Expecting a numpy array");
+ pyerror("Expecting a numpy array for %s", name);
}
/* check if required dim */
if (PyArray_NDIM(obj) != dim)
{
- pyerror("Array should be a %i dimensional object", dim);
+ pyerror("Array should be a %i dimensional object for %s", dim, name);
}
/* check data type */
if (PyArray_TYPE(obj) != type)
{
- pyerror("Wrong array type");
+ pyerror("Wrong array type for %s", name);
}
return SUCCESS;
diff --git a/src/pyswiftsim_tools.h b/src/pyswiftsim_tools.h
index 65cb3df581181723e2df87a3d5ff42f79a27059d..3cc461c5019251511dfc91d04b73da64d3b49d5b 100644
--- a/src/pyswiftsim_tools.h
+++ b/src/pyswiftsim_tools.h
@@ -1,22 +1,33 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
#ifndef __PYSWIFTSIM_TOOLS_H__
#define __PYSWIFTSIM_TOOLS_H__
-#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
-
-#include <Python.h>
-#include <numpy/arrayobject.h>
-
-#undef assert
-#include <swift.h>
+#include "pyswiftsim.h"
#define DIM 3
#define STRING_SIZE 200
-/* Set the error message for python (still need to return NULL) */
+/* Set the error message for python */
#define pyerror(s, ...) \
({ \
char error_msg[STRING_SIZE]; \
- PyErr_Print(); \
sprintf(error_msg, "\n%s:%s():%i: " s, __FILE__, \
__FUNCTION__, __LINE__, ##__VA_ARGS__); \
PyErr_SetString(PyExc_RuntimeError, error_msg); \
@@ -26,7 +37,8 @@
/* debugging function */
#define pydebug(s, ...) \
({ \
- PyErr_Print(); \
+ if (!PyErr_Occurred()) \
+ PyErr_Print(); \
fprintf(stderr, "\n%s:%s():%i: Test " s "\n", __FILE__, \
__FUNCTION__, __LINE__, ##__VA_ARGS__); \
})
@@ -49,6 +61,8 @@ enum swift_class {
class_swift_params,
class_phys_const,
class_cooling_function_data,
+ class_chemistry_global_data,
+ class_cosmology,
class_count /* should always be last! */
};
@@ -64,49 +78,14 @@ enum error_code {
};
-/**
- * @brief Construct a python object from a C struct
- *
- * @param p pointer to the C struct
- * @param class #swift_class of the pointer
- * @return PyObject of the required class
- */
PyObject* pytools_return(void* p, int class);
-/**
- * @brief Construct a C struct from a python object
- *
- * @param obj pointer to the python object
- * @param class #swift_class of the pointer
- * @return pointer to the required struct
- */
char* pytools_construct(PyObject* obj, int class);
-/**
- * @brief Import a python library/function
- *
- * @param module Module name
- * @param object_name object to import from module
- * @return (PyObject) imported object
- */
PyObject* pytools_import(char* module, char* object_name);
-/**
- * @brief get type name in C
- *
- * @param obj PyObject from which to get the name0
- * @return type name
- */
char* pytools_get_type_name(PyObject *obj);
-/**
- * @brief Check if object is the expected PyArray type
- *
- * @param obj PyArray to check
- * @param dim required dimension
- * @param int required type
- * @return #error_code result of the test
- */
-int pytools_check_array(PyArrayObject *obj, int dim, int type);
+int pytools_check_array(PyArrayObject *obj, int dim, int type, const char* name);
#endif // __PYSWIFTSIM_TOOLS_H__
diff --git a/src/units_wrapper.c b/src/units_wrapper.c
index fe3a23a0cac0ed011cbaf9417d77c0d933e4721c..5fc066f1e28e155e9cf268a50b564ab4ae8d9a42 100644
--- a/src/units_wrapper.c
+++ b/src/units_wrapper.c
@@ -1,3 +1,21 @@
+/*******************************************************************************
+ * 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 "pyswiftsim_tools.h"
#include <Python.h>
@@ -20,8 +38,6 @@ PyObject* pyunit_system_test_struct(PyObject *self, PyObject *args)
/* construct python object */
PyObject *object = pytools_return(us, class_unit_system);
free(us);
- if (object == NULL)
- return NULL;
return object;
}
@@ -29,9 +45,10 @@ PyObject* pyunit_system_test_struct(PyObject *self, PyObject *args)
PyObject* pyunit_system_init(PyObject *self, PyObject *args)
{
PyObject* parser;
+ char* section = "InternalUnitSystem";
/* parse arguement */
- if (!PyArg_ParseTuple(args, "O", &parser))
+ if (!PyArg_ParseTuple(args, "O|s", &parser, §ion))
return NULL;
struct swift_params *params = (struct swift_params*) pytools_construct(parser, class_swift_params);
@@ -41,11 +58,12 @@ PyObject* pyunit_system_init(PyObject *self, PyObject *args)
/* initialize units */
struct unit_system us;
- units_init(&us, params, "InternalUnitSystem");
+ units_init_from_params(&us, params, section);
/* initialize phys_const */
struct phys_const pconst;
- phys_const_init(&us, &pconst);
+
+ phys_const_init(&us, params, &pconst);
/* create python object to return */
PyObject *pyus = pytools_return(&us, class_unit_system);
@@ -58,5 +76,3 @@ PyObject* pyunit_system_init(PyObject *self, PyObject *args)
return PyTuple_Pack(2, pyus, pypconst);
}
-
-
diff --git a/src/units_wrapper.h b/src/units_wrapper.h
index 867e055863d24203dbfb41213ef6fd837252a28e..fa42a09dda6961ba1065dbd65bc785c3cacc04ea 100644
--- a/src/units_wrapper.h
+++ b/src/units_wrapper.h
@@ -1,3 +1,21 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
#ifndef __PYSWIFTSIM_UNITS_H__
#define __PYSWIFTSIM_UNITS_H__
diff --git a/src/wrapper.c b/src/wrapper.c
index 46076c5cb627eaed51654bb02107fcaeefcbbfb6..65dbad5592af49b1ca4900b3c99ad357edd56031 100644
--- a/src/wrapper.c
+++ b/src/wrapper.c
@@ -1,9 +1,28 @@
+/*******************************************************************************
+ * 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 "units_wrapper.h"
#include "part_wrapper.h"
#include "parser_wrapper.h"
#include "cooling_wrapper.h"
+#include "cosmology_wrapper.h"
+#include "chemistry_wrapper.h"
#include "pyswiftsim_tools.h"
-#include "config_wrapper.h"
#include <Python.h>
#include <math.h>
@@ -14,20 +33,17 @@
static PyMethodDef wrapper_methods[] = {
- {"partTestStruct", pypart_test_struct, METH_VARARGS,
- "Construct a part object and return it."},
-
- {"unitSystemTestStruct", pyunit_system_test_struct, METH_VARARGS,
- "Construct a unit_system object and return it."},
-
{"parserReadFile", pyparser_read_file, METH_VARARGS,
"Read a swift params file."},
{"unitSystemInit", pyunit_system_init, METH_VARARGS,
"Construct a unit_system object and return it."},
- {"coolingInit", pycooling_init, METH_VARARGS,
- "Initialize cooling."},
+ {"chemistryInit", pychemistry_init, METH_VARARGS,
+ "Initialize chemistry."},
+
+ {"cosmologyInit", pycosmology_init, METH_VARARGS,
+ "Initialize cosmology."},
{"coolingRate", pycooling_rate, METH_VARARGS,
"Compute the cooling rate.\n\n"
@@ -44,11 +60,13 @@ static PyMethodDef wrapper_methods[] = {
"\t Time step in pyus units\n"
"fractions: np.array, optional\n"
"\t Fraction of each cooling element (including metals)\n"
+ "Returns\n"
+ "-------\n\n"
+ "rate: np.array\n"
+ "\t Cooling rate\n"
},
- {"configGetCooling", config_get_cooling, METH_VARARGS,
- "Get the cooling type."},
-
+
{NULL, NULL, 0, NULL} /* Sentinel */
};
@@ -59,11 +77,15 @@ static struct PyModuleDef wrapper_cmodule = {
"wrapper",
"Wrapper around the SPH cosmological simulation code SWIFT",
-1,
- wrapper_methods
+ wrapper_methods,
+ NULL, /* m_slots */
+ NULL, /* m_traverse */
+ NULL, /* m_clear */
+ NULL, /* m_free */
};
-PyMODINIT_FUNC PyInit_wrapper(void)
+PyMODINIT_FUNC PyInit_libpyswiftsim(void)
{
PyObject *m;
@@ -74,8 +96,5 @@ PyMODINIT_FUNC PyInit_wrapper(void)
Py_Initialize();
m = PyModule_Create(&wrapper_cmodule);
- if (m == NULL)
- return NULL;
-
return m;
}
diff --git a/test/test_struct.py b/test/test_struct.py
index 7c228d849979d2229cfd9d741862f54535f63744..2d9df52e88b5a247619ce1df6065bf7c72b6076f 100644
--- a/test/test_struct.py
+++ b/test/test_struct.py
@@ -1,46 +1,29 @@
#!/usr/bin/env python3
+"""
+This file is part of PYSWIFT.
+Copyright (c) 2018 Loic Hausammann (loic.hausammann@epfl.ch)
-from pyswiftsim import wrapper
-from pyswiftsim import structure
-
-us = wrapper.unitSystemTestStruct()
-
-print(us)
-
-part = wrapper.partTestStruct()
-
-print(part)
-
-#us = structure.Units(t)
-
-#print(type(us), us.struct_name)
-#print(us.struct_size_format)
+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.
-#for i in us.struct_name:
-# print(i, getattr(us, i))
+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/>.
"""
-print(us.pos)
-us.pos = [12., 16., 11.]
-#us.pos[0] = 10.
-
-print(us)
-print(us.pos)
-print("setitem")
-us.pos[2] = 20.
+from pyswiftsim import wrapper
-print(us.pos)
-print(us)
+filename = "parameter_example.yml"
+params = wrapper.parserReadFile(filename)
-d = {}
-d["UnitMass_in_cgs"] = 1.
-d["UnitLength_in_cgs"] = 2.
-d["UnitTime_in_cgs"] = 3.
-d["UnitCurrent_in_cgs"] = 4.
-d["UnitTemperature_in_cgs"] = 5.
+print(params)
-us = structure.UnitSystem(d)
+us = wrapper.unitSystemInit(params)
print(us)
-"""