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..4daa85e84e708c56e89190d458e47d6af9f0edcb
--- /dev/null
+++ b/Makefile.am
@@ -0,0 +1,22 @@
+# This file is part of SWIFT.
+# Copyright (c) 2012 pedro.gonnet@durham.ac.uk
+# 2015 matthieu.schaller@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/>.
+
+# 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/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_rate/cooling.yml b/examples/cooling_rate/cooling.yml
index 0dd7d8b65ab0d4f54be66ce743182350b889b230..ba04594f9b4abd193f53f2438f299aec701e16d1 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:
diff --git a/examples/cooling_rate/plot_cooling.py b/examples/cooling_rate/plot_cooling.py
index 4471dea1043201c857e68f7f0b9cc91315da4583..a8096db72053d0b529449808a3685d769da2f87b 100644
--- a/examples/cooling_rate/plot_cooling.py
+++ b/examples/cooling_rate/plot_cooling.py
@@ -17,21 +17,31 @@ 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
+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 (for comparison)
-primordial_chemistry = 3
+primordial_chemistry = 0
+
+# grackle simulation
+grackle_simulation = 0
+grackle_simulation = "Simulation {}".format(grackle_simulation)
+
+# name of this simulation
+simulation_name = "Grackle"
# reference data
grackle_filename = "grackle.hdf5"
@@ -98,7 +108,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
@@ -108,88 +118,50 @@ 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 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,
- }
-
- 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, primordial_chemistry, cooling, chemistry):
+def read_grackle_data(filename, us):
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]]
+ data = data[grackle_simulation]
# 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
@@ -215,32 +187,10 @@ def read_grackle_data(filename, us, primordial_chemistry, cooling, 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 cosmo
- 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 getParser(filename):
+ with open(filename, "r") as stream:
+ stream = yaml.load(stream)
+ return stream
def plot_solution(rate, data, us):
@@ -255,16 +205,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
@@ -274,13 +227,10 @@ def plot_solution(rate, data, us):
# do plot
if N_rho == 1:
- # save data
- write_solution(T, rho, lam)
-
# 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),
@@ -334,81 +284,25 @@ def plot_solution(rate, data, us):
plt.show()
-def write_solution(T, rho, lam):
- f = File(swift_output_filename, "a")
-
- # create data set for this cooling
- data_name = wrapper.configGetCooling()
-
- if data_name not in f:
- print("Creating Dataset")
- data = f.create_group(data_name)
- else:
- data = f[data_name]
-
- i = getSimulationIndice(data, cooling, chemistry, output=True)
- if i != -1:
- print("Updating Simulation")
- keys = list(data.keys())
- data.pop(keys[i])
-
- data = data.create_group("Simulation %i" % len(data))
-
- tmp = data.create_group("Units")
- tmp.attrs["UnitSystem"] = "cgs"
-
- tmp = data.create_group("Params")
- fields = get_io_fields(cooling, chemistry, output=True)
- for key in fields.keys():
- tmp.attrs[key] = fields[key]
-
- # temperature
- dset = data.create_dataset(
- "Temperature", T.shape, dtype=T.dtype)
- dset[:] = T
-
- # density
- dset = data.create_dataset(
- "Density", rho.shape, dtype=rho.dtype)
- dset[:] = rho
-
- # cooling rate
- dset = data.create_dataset(
- "Lambda", lam.shape, dtype=lam.dtype)
- dset[:] = lam
-
- f.close()
-
-
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"]
+ parser = getParser(filename)
+ us = parser["InternalUnitSystem"]
if isfile(grackle_filename) and N_rho == 1:
- d = read_grackle_data(grackle_filename, us,
- primordial_chemistry, cooling,
- chemistry)
+ 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,
- cosmo,
- cooling,
- chemistry,
- 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"])
else:
fields = get_fields(primordial_chemistry)
N = len(fields)
@@ -416,13 +310,10 @@ if __name__ == "__main__":
for i in range(N):
frac[:, i] = d[fields[i]]
- rate = wrapper.coolingRate(pconst, us,
- cosmo,
- cooling,
- chemistry,
- 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/src/Makefile.am b/src/Makefile.am
index 40c61f96f1dbd77e1caebebedd95ffbe8d61ac05..a77473fd9cf45fb642c0175a3d22490b9b7f9084 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -15,7 +15,7 @@
# 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)
+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
diff --git a/src/cooling_wrapper.c b/src/cooling_wrapper.c
index 242dc7f76930d1bcafc734a0e4aff3830d73380b..8def5d2d1e23507a0edd4f438816df7a02aff522 100644
--- a/src/cooling_wrapper.c
+++ b/src/cooling_wrapper.c
@@ -16,53 +16,12 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
-#include "pyswiftsim_tools.h"
#include "cooling_wrapper.h"
-#include <omp.h>
-
-/**
- * @brief Initialize the cooling
- *
- * args is expecting pyswiftsim classes in the following order:
- * SwiftParams, UnitSystem and PhysConst.
- *
- * @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;
-
- /* init cooling */
- cooling_init_backend(params, us, pconst, &cooling);
-
- /* construct python object */
- PyObject *pycooling = pytools_return(&cooling, class_cooling_function_data);
+/* Local includes */
+#include "pyswiftsim_tools.h"
+#include "part_wrapper.h"
- return pycooling;
-}
/**
* @brief Set the cooling element fractions
@@ -117,31 +76,20 @@ void pycooling_set_fractions(struct xpart *xp, PyArrayObject* frac, const int id
PyObject* pycooling_rate(PyObject* self, PyObject* args) {
import_array();
- PyObject *pycooling;
- PyObject *pyus;
- PyObject *pycosmo;
- PyObject *pypconst;
- PyObject *pychemistry;
+
+ char *filename;
PyArrayObject *rho;
PyArrayObject *energy;
PyArrayObject *fractions = NULL;
float dt = 1e-3;
+ int with_cosmo = 0;
/* parse argument */
- if (!PyArg_ParseTuple(args,
- "OOOOOOO|fO",
- &pypconst,
- &pyus,
- &pycosmo,
- &pycooling,
- &pychemistry,
- &rho,
- &energy,
- &dt,
- &fractions
- ))
+ if (!PyArg_ParseTuple(
+ args, "siOO|fO", &filename, &with_cosmo,
+ &rho, &energy, &dt, &fractions))
return NULL;
/* check numpy array */
@@ -164,37 +112,41 @@ PyObject* 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 chemistry_global_data *chemistry = (struct chemistry_global_data*) pytools_construct(pychemistry, class_chemistry_global_data);
- if (chemistry == NULL)
- return NULL;
+ struct cooling_function_data cooling;
- struct cosmology *cosmo = (struct cosmology*) pytools_construct(pycosmo, class_cosmology);
- if (cosmo == NULL)
- return NULL;
+ /* 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;
+ pyswift_part_set_to_zero(&p, &xp);
+ hydro_first_init_part(&p, &xp);
+
/* return object */
PyArrayObject *rate = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(energy), PyArray_DIMS(energy), NPY_FLOAT);
- /* Release GIL */
- Py_BEGIN_ALLOW_THREADS;
-
/* loop over all particles */
for(size_t i = 0; i < N; i++)
{
@@ -206,159 +158,22 @@ PyObject* pycooling_rate(PyObject* self, PyObject* args) {
if (fractions != NULL)
pycooling_set_fractions(&xp, fractions, i);
else
- cooling_first_init_part(pconst, us, cosmo, cooling, &p, &xp);
+ cooling_first_init_part(&pconst, &us, &cosmo, &cooling, &p, &xp);
- chemistry_first_init_part(pconst, us, cosmo, chemistry, &p, &xp);
+ chemistry_first_init_part(&pconst, &us, &cosmo, &chemistry, &p, &xp);
/* compute cooling rate */
#ifdef COOLING_GRACKLE
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 = cooling_new_energy(&pconst, &us, &cosmo, &cooling, &p, &xp, dt);
+ *tmp = *tmp - hydro_get_physical_internal_energy(&p, &xp, &cosmo);
*tmp /= dt;
- message("%g", *tmp);
#else
message("Not implemented");
//*tmp = cooling_rate(pconst, us, cosmo, cooling, &p, &xp);
#endif
}
- /* Acquire GIL */
- Py_END_ALLOW_THREADS;
-
return (PyObject *) rate;
}
-
-
-/**
- * @brief Compute cooling
- *
- * args is expecting pyswiftsim classes in the following order:
- * PhysConst, UnitSystem, Cosmology and CoolingFunctionData.
- * Then two numpy arrays (density and specific energy) and a
- * float for the time step
- *
- * @param self calling object
- * @param args arguments
- * @return New energy
- */
-PyObject* pycooling_do_cooling(PyObject* self, PyObject* args) {
- import_array();
-
- PyObject *pycooling;
- PyObject *pyus;
- PyObject *pycosmo;
- PyObject *pypconst;
- PyObject *pychemistry;
-
- PyArrayObject *rho;
- PyArrayObject *energy;
- PyArrayObject *fractions = NULL;
-
- float dt = 1e-3;
-
- /* parse argument */
- if (!PyArg_ParseTuple(args,
- "OOOOOOOf|O",
- &pypconst,
- &pyus,
- &pycosmo,
- &pycooling,
- &pychemistry,
- &rho,
- &energy,
- &dt,
- &fractions
- ))
- return NULL;
-
- /* check numpy array */
- if (pytools_check_array(energy, 1, NPY_FLOAT, "Energy") != SUCCESS)
- return NULL;
-
- if (pytools_check_array(rho, 1, NPY_FLOAT, "Density") != SUCCESS)
- return NULL;
-
- if (fractions != NULL &&
- pytools_check_array(fractions, 2, NPY_FLOAT, "Fractions") != SUCCESS)
- return NULL;
-
- if (PyArray_DIM(energy, 0) != PyArray_DIM(rho, 0))
- pyerror("Density and energy should have the same dimension");
-
- if (fractions != NULL &&
- PyArray_DIM(fractions, 0) != PyArray_DIM(rho,0))
- pyerror("Fractions should have the same first dimension than the others");
-
- 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;
-
- 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 = (struct chemistry_global_data*) pytools_construct(pychemistry, class_chemistry_global_data);
- if (chemistry == NULL)
- return NULL;
-
- struct cosmology *cosmo = (struct cosmology*) pytools_construct(pycosmo, class_cosmology);
- if (cosmo == NULL)
- return NULL;
-
- struct part p;
- struct xpart xp;
-
-#ifdef COOLING_GRACKLE
- /* set grackle_data */
- grackle_data = &cooling->chemistry;
-#endif
-
- /* return object */
- PyArrayObject *new_energy = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(energy), PyArray_DIMS(energy), NPY_FLOAT);
-
- /* Release GIL */
- Py_BEGIN_ALLOW_THREADS;
-
- /* loop over all particles */
- for(size_t i = 0; i < N; i++)
- {
- /* set particle data */
- p.rho = *(float*) PyArray_GETPTR1(rho, i);
- float u = *(float*) PyArray_GETPTR1(energy, i);
- p.entropy = gas_entropy_from_internal_energy(p.rho, u);
-
- /* initialize fields */
- //hydro_set_internal_energy_dt(&p, 0.);
-
- chemistry_first_init_part(pconst, us, cosmo, chemistry, &p, &xp);
-
- if (fractions != NULL)
- pycooling_set_fractions(&xp, fractions, i);
- else
- cooling_first_init_part(pconst, us, cosmo, cooling, &p, &xp);
-
- /* compute cooling */
- cooling_cool_part(pconst, us, cosmo, NULL, cooling, &p, &xp, dt, dt);
- message("Need to redo");
-
- /* update energy */
- float *tmp = PyArray_GETPTR1(new_energy, i);
- *tmp = u + dt * hydro_get_physical_internal_energy_dt(&p, NULL);
- }
-
- /* Acquire GIL */
- Py_END_ALLOW_THREADS;
-
- return (PyObject *) new_energy;
-
-}
diff --git a/src/cooling_wrapper.h b/src/cooling_wrapper.h
index e80cb1c38f7d97d2b00ac7161d8c98f9cbcaf578..a9eeeb4d6a5061637795536f76af71cc439c6747 100644
--- a/src/cooling_wrapper.h
+++ b/src/cooling_wrapper.h
@@ -21,11 +21,7 @@
#include "pyswiftsim_tools.h"
-PyObject* pycooling_init(PyObject* self, PyObject* args);
-
PyObject* pycooling_rate(PyObject* self, PyObject* args);
-PyObject* pycooling_do_cooling(PyObject* self, PyObject* args);
-
#endif // __PYSWIFTSIM_COOLING_H__
diff --git a/src/part_wrapper.c b/src/part_wrapper.c
index 65ceec2618a250e4faeae64a9243179c83539fc0..d8cefd848401a1e5eedb5839e5f12188e2c7738c 100644
--- a/src/part_wrapper.c
+++ b/src/part_wrapper.c
@@ -22,34 +22,9 @@
#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 05a2fd8b4b57d1f8847a0fc4757e287cc00ae740..bdc3d52b5807871a4c0e656150686032917f17d8 100644
--- a/src/part_wrapper.h
+++ b/src/part_wrapper.h
@@ -20,16 +20,8 @@
#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.h b/src/pyswiftsim_tools.h
index ca2b0d03129c88da287f71e7ccb78342fb35c094..3cc461c5019251511dfc91d04b73da64d3b49d5b 100644
--- a/src/pyswiftsim_tools.h
+++ b/src/pyswiftsim_tools.h
@@ -19,14 +19,7 @@
#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 "../config.h"
-#include "swift.h"
+#include "pyswiftsim.h"
#define DIM 3
#define STRING_SIZE 200
diff --git a/src/wrapper.c b/src/wrapper.c
index 2bc1f5f830bce9fd56773c3017f069394b0bae92..65dbad5592af49b1ca4900b3c99ad357edd56031 100644
--- a/src/wrapper.c
+++ b/src/wrapper.c
@@ -23,7 +23,6 @@
#include "cosmology_wrapper.h"
#include "chemistry_wrapper.h"
#include "pyswiftsim_tools.h"
-#include "swift.h"
#include <Python.h>
#include <math.h>
@@ -40,9 +39,6 @@ static PyMethodDef wrapper_methods[] = {
{"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."},
@@ -70,26 +66,6 @@ static PyMethodDef wrapper_methods[] = {
"\t Cooling rate\n"
},
- {"doCooling", pycooling_do_cooling, METH_VARARGS,
- "Compute the cooling.\n\n"
- "Parameters\n"
- "----------\n\n"
- "pyconst: swift physical constant\n"
- "pyus: swift unit system\n"
- "cooling: swift cooling structure\n"
- "rho: np.array\n"
- "\t Mass density in pyus units\n"
- "energy: np.array\n"
- "\t Internal energy in pyus units\n"
- "dt: float\n"
- "\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"
- "new_energy: np.array\n"
- "\t Energy of the particle after dt\n"
- },
{NULL, NULL, 0, NULL} /* Sentinel */
};