From 083548de07e9dd7ccb17d962bb3e253d656c0a83 Mon Sep 17 00:00:00 2001
From: loikki <loic.hausammann@protonmail.ch>
Date: Wed, 13 Feb 2019 09:00:55 +0100
Subject: [PATCH] Start implementing PySWIFTsim as a subproject of SWIFT
---
AUTHORS | 2 +
COPYING | 1 +
ChangeLog | 0
INSTALL | 1 +
Makefile.am | 22 +++
NEWS | 0
README | 1 +
configure.ac | 145 +++++++++++++++
examples/cooling_rate/cooling.yml | 11 ++
examples/cooling_rate/plot_cooling.py | 229 ++++++-----------------
src/Makefile.am | 2 +-
src/cooling_wrapper.c | 257 ++++----------------------
src/cooling_wrapper.h | 4 -
src/part_wrapper.c | 35 +---
src/part_wrapper.h | 14 +-
src/pyswiftsim.h | 32 ++++
src/pyswiftsim_tools.h | 9 +-
src/wrapper.c | 24 ---
18 files changed, 321 insertions(+), 468 deletions(-)
create mode 100644 AUTHORS
create mode 120000 COPYING
create mode 100644 ChangeLog
create mode 120000 INSTALL
create mode 100644 Makefile.am
create mode 100644 NEWS
create mode 120000 README
create mode 100644 configure.ac
create mode 100644 src/pyswiftsim.h
diff --git a/AUTHORS b/AUTHORS
new file mode 100644
index 0000000..d202675
--- /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 0000000..88798ab
--- /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 0000000..e69de29
diff --git a/INSTALL b/INSTALL
new file mode 120000
index 0000000..ddcdb76
--- /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 0000000..4daa85e
--- /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 0000000..e69de29
diff --git a/README b/README
new file mode 120000
index 0000000..602d87d
--- /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 0000000..f9353c5
--- /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 0dd7d8b..ba04594 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 4471dea..a8096db 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 40c61f9..a77473f 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 242dc7f..8def5d2 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 e80cb1c..a9eeeb4 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 65ceec2..d8cefd8 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 05a2fd8..bdc3d52 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 0000000..f215e81
--- /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 ca2b0d0..3cc461c 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 2bc1f5f..65dbad5 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 */
};
--
GitLab