Skip to content
Snippets Groups Projects
Commit aa642cd6 authored by Loic Hausammann's avatar Loic Hausammann
Browse files

Merge branch 'test' of gitlab.cosma.dur.ac.uk:swift/PySWIFTsim into test

parents 5b92dd29 083548de
No related branches found
No related tags found
1 merge request!5Last commit before huge rework
Loic Hausammann loic.hausammann@epfl.ch
Yves Revaz yves.revaz@epfl.ch
/usr/share/automake-1.15/COPYING
\ No newline at end of file
/usr/share/automake-1.15/INSTALL
\ No newline at end of file
...@@ -19,6 +19,3 @@ ACLOCAL_AMFLAGS = -I m4 ...@@ -19,6 +19,3 @@ ACLOCAL_AMFLAGS = -I m4
# Show the way... # Show the way...
SUBDIRS = src SUBDIRS = src
# Non-standard files that should be part of the distribution.
# EXTRA_DIST = INSTALL.swift .clang-format format.sh
NEWS 0 → 100644
README 0 → 120000
Readme.md
\ No newline at end of file
# 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
...@@ -34,6 +34,17 @@ InitialConditions: ...@@ -34,6 +34,17 @@ InitialConditions:
file_name: ./sedov.hdf5 # The file to read file_name: ./sedov.hdf5 # The file to read
h_scaling: 3.33 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 # Cooling with Grackle 2.0
GrackleCooling: GrackleCooling:
......
...@@ -17,18 +17,21 @@ You should have received a copy of the GNU Lesser General Public License ...@@ -17,18 +17,21 @@ 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/>. along with this program. If not, see <http://www.gnu.org/licenses/>.
""" """
from pyswiftsim import wrapper from libpyswiftsim import coolingRate
from copy import deepcopy from copy import deepcopy
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from astropy import units from astropy import units, constants
from os.path import isfile from os.path import isfile
from h5py import File from h5py import File
import yaml
plt.rc('text', usetex=True) plt.rc('text', usetex=True)
# PARAMETERS # PARAMETERS
# cosmology
with_cosmo = 0
# grackle primordial chemistry (for comparison) # grackle primordial chemistry (for comparison)
primordial_chemistry = 0 primordial_chemistry = 0
...@@ -98,7 +101,7 @@ def get_fields(primordial_chemistry): ...@@ -98,7 +101,7 @@ def get_fields(primordial_chemistry):
return fields return fields
def generate_default_initial_condition(us, pconst): def generate_default_initial_condition(us):
print("Generating default initial conditions") print("Generating default initial conditions")
d = {} d = {}
# generate grid # generate grid
...@@ -108,67 +111,32 @@ def generate_default_initial_condition(us, pconst): ...@@ -108,67 +111,32 @@ def generate_default_initial_condition(us, pconst):
d["temperature"] = T d["temperature"] = T
# Deal with units # 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 d["density"] = rho
energy = pconst.const_boltzmann_k * T / us.UnitTemperature_in_cgs unit_time = us["UnitLength_in_cgs"] / us["UnitVelocity_in_cgs"]
energy /= (gamma - 1.) * pconst.const_proton_mass 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 d["energy"] = energy
dt = default_dt / us.UnitTime_in_cgs dt = default_dt / unit_time
d["time_step"] = dt d["time_step"] = dt
return d return d
def get_io_fields(cooling, chemistry, output): def read_grackle_data(filename, us):
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):
print("Reading initial conditions") print("Reading initial conditions")
f = File(filename, "r") f = File(filename, "r")
data = f["PrimordialChemistry%i" % primordial_chemistry] data = f["PrimordialChemistry%i" % primordial_chemistry]
<<<<<<< HEAD
print(cooling) print(cooling)
i = getSimulationIndice(data, cooling, chemistry, output=False) i = getSimulationIndice(data, cooling, chemistry, output=False)
if i < 0: if i < 0:
...@@ -176,21 +144,27 @@ def read_grackle_data(filename, us, primordial_chemistry, cooling, chemistry): ...@@ -176,21 +144,27 @@ def read_grackle_data(filename, us, primordial_chemistry, cooling, chemistry):
else: else:
tmp = list(data.keys()) tmp = list(data.keys())
data = data[tmp[i]] data = data[tmp[i]]
=======
data = data[grackle_simulation]
>>>>>>> 083548de07e9dd7ccb17d962bb3e253d656c0a83
# read units # read units
tmp = data["Units"].attrs tmp = data["Units"].attrs
u_len = tmp["Length"] / us.UnitLength_in_cgs u_len = tmp["Length"] / float(us["UnitLength_in_cgs"])
u_den = tmp["Density"] * us.UnitLength_in_cgs**3 / us.UnitMass_in_cgs u_den = tmp["Density"] * float(us["UnitLength_in_cgs"])**3
u_time = tmp["Time"] / us.UnitTime_in_cgs 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 # read input
tmp = data["Input"] tmp = data["Input"]
d = {}
energy = tmp["Energy"][:] * u_len**2 / u_time**2 energy = tmp["Energy"][:] * u_len**2 / u_time**2
d["energy"] = energy d["energy"] = energy
T = tmp["Temperature"][:] / us.UnitTemperature_in_cgs T = tmp["Temperature"][:] / float(us["UnitTemp_in_cgs"])
d["temperature"] = T d["temperature"] = T
density = tmp["Density"][:] * u_den density = tmp["Density"][:] * u_den
...@@ -216,32 +190,10 @@ def read_grackle_data(filename, us, primordial_chemistry, cooling, chemistry): ...@@ -216,32 +190,10 @@ def read_grackle_data(filename, us, primordial_chemistry, cooling, chemistry):
return d return d
def initialize_swift(filename): def getParser(filename):
print("Initialization of SWIFT") with open(filename, "r") as stream:
d = {} stream = yaml.load(stream)
return stream
# 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 plot_solution(rate, data, us): def plot_solution(rate, data, us):
...@@ -256,16 +208,19 @@ def plot_solution(rate, data, us): ...@@ -256,16 +208,19 @@ def plot_solution(rate, data, us):
grackle_rate = data["rate"] grackle_rate = data["rate"]
# change units => cgs # 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: if ref:
grackle_rate *= us.UnitLength_in_cgs**2 / us.UnitTime_in_cgs**3 grackle_rate *= u_l**2 / u_t**3
# lambda cooling # lambda cooling
lam = rate * rho lam = rate * rho
...@@ -275,13 +230,10 @@ def plot_solution(rate, data, us): ...@@ -275,13 +230,10 @@ def plot_solution(rate, data, us):
# do plot # do plot
if N_rho == 1: if N_rho == 1:
# save data
write_solution(T, rho, lam)
# plot Lambda vs T # plot Lambda vs T
plt.figure() plt.figure()
plt.loglog(T, np.abs(lam), plt.loglog(T, np.abs(lam),
label="SWIFT, %s" % wrapper.configGetCooling()) label="SWIFT")
if ref: if ref:
# plot reference # plot reference
plt.loglog(T, np.abs(grackle_lam), plt.loglog(T, np.abs(grackle_lam),
...@@ -335,79 +287,33 @@ def plot_solution(rate, data, us): ...@@ -335,79 +287,33 @@ def plot_solution(rate, data, us):
plt.show() 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__": if __name__ == "__main__":
d = initialize_swift(filename) parser = getParser(filename)
pconst = d["pconst"] us = parser["InternalUnitSystem"]
us = d["us"]
params = d["params"]
cosmo = d["cosmo"]
cooling = d["cooling"]
chemistry = d["chemistry"]
if isfile(grackle_filename) and N_rho == 1: if isfile(grackle_filename) and N_rho == 1:
d = read_grackle_data(grackle_filename, us, d = read_grackle_data(grackle_filename, us)
primordial_chemistry, cooling,
chemistry)
else: else:
d = generate_default_initial_condition(us, pconst) d = generate_default_initial_condition(us)
# du / dt # du / dt
print("Computing cooling...") print("Computing cooling...")
rate = np.zeros(d["density"].shape) rate = np.zeros(d["density"].shape)
if compute_equilibrium: if compute_equilibrium:
<<<<<<< HEAD
rate = wrapper.coolingRate( rate = wrapper.coolingRate(
filename, filename,
d["density"].astype(np.float32), d["density"].astype(np.float32),
d["energy"].astype(np.float32), d["energy"].astype(np.float32),
d["time_step"]) d["time_step"])
=======
rate = coolingRate(filename, with_cosmo,
d["density"].astype(np.float32),
d["energy"].astype(np.float32),
d["time_step"])
>>>>>>> 083548de07e9dd7ccb17d962bb3e253d656c0a83
else: else:
fields = get_fields(primordial_chemistry) fields = get_fields(primordial_chemistry)
N = len(fields) N = len(fields)
...@@ -415,13 +321,10 @@ if __name__ == "__main__": ...@@ -415,13 +321,10 @@ if __name__ == "__main__":
for i in range(N): for i in range(N):
frac[:, i] = d[fields[i]] frac[:, i] = d[fields[i]]
rate = wrapper.coolingRate(pconst, us, rate = coolingRate(filename, with_cosmo,
cosmo, d["density"].astype(np.float32),
cooling, d["energy"].astype(np.float32),
chemistry, d["time_step"],
d["density"].astype(np.float32), frac.astype(np.float32))
d["energy"].astype(np.float32),
d["time_step"],
frac.astype(np.float32))
plot_solution(rate, d, us) plot_solution(rate, d, us)
...@@ -15,7 +15,7 @@ ...@@ -15,7 +15,7 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>. # along with this program. If not, see <http://www.gnu.org/licenses/>.
# Add the non-standard paths to the included library headers # 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 # Assign a "safe" version number
AM_LDFLAGS = $(HDF5_LDFLAGS) $(FFTW_LIBS) -version-info 0:0:0 AM_LDFLAGS = $(HDF5_LDFLAGS) $(FFTW_LIBS) -version-info 0:0:0
......
...@@ -16,53 +16,12 @@ ...@@ -16,53 +16,12 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. * along with this program. If not, see <http://www.gnu.org/licenses/>.
* *
******************************************************************************/ ******************************************************************************/
#include "pyswiftsim_tools.h"
#include "cooling_wrapper.h" #include "cooling_wrapper.h"
#include <omp.h>
/** /* Local includes */
* @brief Initialize the cooling #include "pyswiftsim_tools.h"
* #include "part_wrapper.h"
* 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);
return pycooling;
}
/** /**
* @brief Set the cooling element fractions * @brief Set the cooling element fractions
...@@ -117,31 +76,20 @@ void pycooling_set_fractions(struct xpart *xp, PyArrayObject* frac, const int id ...@@ -117,31 +76,20 @@ void pycooling_set_fractions(struct xpart *xp, PyArrayObject* frac, const int id
PyObject* pycooling_rate(PyObject* self, PyObject* args) { PyObject* pycooling_rate(PyObject* self, PyObject* args) {
import_array(); import_array();
PyObject *pycooling;
PyObject *pyus; char *filename;
PyObject *pycosmo;
PyObject *pypconst;
PyObject *pychemistry;
PyArrayObject *rho; PyArrayObject *rho;
PyArrayObject *energy; PyArrayObject *energy;
PyArrayObject *fractions = NULL; PyArrayObject *fractions = NULL;
float dt = 1e-3; float dt = 1e-3;
int with_cosmo = 0;
/* parse argument */ /* parse argument */
if (!PyArg_ParseTuple(args, if (!PyArg_ParseTuple(
"OOOOOOO|fO", args, "siOO|fO", &filename, &with_cosmo,
&pypconst, &rho, &energy, &dt, &fractions))
&pyus,
&pycosmo,
&pycooling,
&pychemistry,
&rho,
&energy,
&dt,
&fractions
))
return NULL; return NULL;
/* check numpy array */ /* check numpy array */
...@@ -164,37 +112,41 @@ PyObject* pycooling_rate(PyObject* self, PyObject* args) { ...@@ -164,37 +112,41 @@ PyObject* pycooling_rate(PyObject* self, PyObject* args) {
size_t N = PyArray_DIM(energy, 0); 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); /* Initialize everything */
if (cooling == NULL) struct swift_params params;
return NULL; parser_read_file(filename, &params);
struct unit_system *us = (struct unit_system*) pytools_construct(pyus, class_unit_system); struct unit_system us;
if (us == NULL) units_init_from_params(&us, &params, "InternalUnitSystem");
return NULL;
struct phys_const *pconst = (struct phys_const*) pytools_construct(pypconst, class_phys_const); struct phys_const pconst;
if (pconst == NULL) phys_const_init(&us, &params, &pconst);
return NULL;
struct chemistry_global_data *chemistry = (struct chemistry_global_data*) pytools_construct(pychemistry, class_chemistry_global_data); struct cooling_function_data cooling;
if (chemistry == NULL)
return NULL;
struct cosmology *cosmo = (struct cosmology*) pytools_construct(pycosmo, class_cosmology); /* init cooling */
if (cosmo == NULL) cooling_init_backend(&params, &us, &pconst, &cooling);
return NULL;
struct chemistry_global_data chemistry;
chemistry_init_backend(&params, &us, &pconst, &chemistry);
struct cosmology cosmo;
if (with_cosmo)
cosmology_init(&params, &us, &pconst, &cosmo);
else
cosmology_init_no_cosmo(&cosmo);
struct part p; struct part p;
struct xpart xp; struct xpart xp;
pyswift_part_set_to_zero(&p, &xp);
hydro_first_init_part(&p, &xp);
/* return object */ /* return object */
PyArrayObject *rate = (PyArrayObject *)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 */ /* loop over all particles */
for(size_t i = 0; i < N; i++) for(size_t i = 0; i < N; i++)
{ {
...@@ -206,159 +158,22 @@ PyObject* pycooling_rate(PyObject* self, PyObject* args) { ...@@ -206,159 +158,22 @@ PyObject* pycooling_rate(PyObject* self, PyObject* args) {
if (fractions != NULL) if (fractions != NULL)
pycooling_set_fractions(&xp, fractions, i); pycooling_set_fractions(&xp, fractions, i);
else 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 */ /* compute cooling rate */
#ifdef COOLING_GRACKLE #ifdef COOLING_GRACKLE
float *tmp = PyArray_GETPTR1(rate, i); float *tmp = PyArray_GETPTR1(rate, i);
*tmp = cooling_new_energy(pconst, us, cosmo, cooling, &p, &xp, dt); *tmp = cooling_new_energy(&pconst, &us, &cosmo, &cooling, &p, &xp, dt);
*tmp = *tmp - hydro_get_physical_internal_energy(&p, &xp, cosmo); *tmp = *tmp - hydro_get_physical_internal_energy(&p, &xp, &cosmo);
*tmp /= dt; *tmp /= dt;
message("%g", *tmp);
#else #else
message("Not implemented"); message("Not implemented");
//*tmp = cooling_rate(pconst, us, cosmo, cooling, &p, &xp); //*tmp = cooling_rate(pconst, us, cosmo, cooling, &p, &xp);
#endif #endif
} }
/* Acquire GIL */
Py_END_ALLOW_THREADS;
return (PyObject *) rate; 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;
}
...@@ -21,11 +21,7 @@ ...@@ -21,11 +21,7 @@
#include "pyswiftsim_tools.h" #include "pyswiftsim_tools.h"
PyObject* pycooling_init(PyObject* self, PyObject* args);
PyObject* pycooling_rate(PyObject* self, PyObject* args); PyObject* pycooling_rate(PyObject* self, PyObject* args);
PyObject* pycooling_do_cooling(PyObject* self, PyObject* args);
#endif // __PYSWIFTSIM_COOLING_H__ #endif // __PYSWIFTSIM_COOLING_H__
...@@ -22,34 +22,9 @@ ...@@ -22,34 +22,9 @@
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
PyObject* pypart_test_struct(PyObject *self, PyObject *args) void pyswift_part_set_to_zero(struct part *p, struct xpart *xp) {
{ p->v[0] = 0;
p->v[1] = 0;
size_t N = sizeof(struct part); p->v[2] = 0;
p->entropy = 0;
/* 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;
} }
...@@ -20,16 +20,8 @@ ...@@ -20,16 +20,8 @@
#define __PYSWIFTSIM_PART_H__ #define __PYSWIFTSIM_PART_H__
#include <Python.h> #include <Python.h>
#include "pyswiftsim_tools.h"
/** void pyswift_part_set_to_zero(struct part *p, struct xpart *xp);
* @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);
#endif // __PYSWIFTSIM_PART_H__ #endif // __PYSWIFTSIM_PART_H__
/*******************************************************************************
* 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__
...@@ -19,14 +19,7 @@ ...@@ -19,14 +19,7 @@
#ifndef __PYSWIFTSIM_TOOLS_H__ #ifndef __PYSWIFTSIM_TOOLS_H__
#define __PYSWIFTSIM_TOOLS_H__ #define __PYSWIFTSIM_TOOLS_H__
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION #include "pyswiftsim.h"
#include <Python.h>
#include <numpy/arrayobject.h>
#undef assert
#include "../config.h"
#include "swift.h"
#define DIM 3 #define DIM 3
#define STRING_SIZE 200 #define STRING_SIZE 200
......
...@@ -23,7 +23,6 @@ ...@@ -23,7 +23,6 @@
#include "cosmology_wrapper.h" #include "cosmology_wrapper.h"
#include "chemistry_wrapper.h" #include "chemistry_wrapper.h"
#include "pyswiftsim_tools.h" #include "pyswiftsim_tools.h"
#include "swift.h"
#include <Python.h> #include <Python.h>
#include <math.h> #include <math.h>
...@@ -40,9 +39,6 @@ static PyMethodDef wrapper_methods[] = { ...@@ -40,9 +39,6 @@ static PyMethodDef wrapper_methods[] = {
{"unitSystemInit", pyunit_system_init, METH_VARARGS, {"unitSystemInit", pyunit_system_init, METH_VARARGS,
"Construct a unit_system object and return it."}, "Construct a unit_system object and return it."},
{"coolingInit", pycooling_init, METH_VARARGS,
"Initialize cooling."},
{"chemistryInit", pychemistry_init, METH_VARARGS, {"chemistryInit", pychemistry_init, METH_VARARGS,
"Initialize chemistry."}, "Initialize chemistry."},
...@@ -70,26 +66,6 @@ static PyMethodDef wrapper_methods[] = { ...@@ -70,26 +66,6 @@ static PyMethodDef wrapper_methods[] = {
"\t Cooling rate\n" "\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 */ {NULL, NULL, 0, NULL} /* Sentinel */
}; };
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment