Commit 14aec30e authored by lhausamm's avatar lhausamm
Browse files

Remove python wrapper

parent 881ad540
from pyswiftsim cimport part_def
from pyswiftsim.part_def cimport part
from pyswiftsim.part_def cimport xpart
cdef extern from "parser.h":
struct swift_params:
pass
void parser_read_file(const char *file_name, swift_params *params);
void parser_print_params(const swift_params *params);
cdef extern from "units.h":
struct unit_system:
double UnitMass_in_cgs;
double UnitLength_in_cgs;
double UnitTime_in_cgs;
double UnitCurrent_in_cgs;
double UnitTemperature_in_cgs;
void units_init(unit_system* us, const swift_params* params,
const char* category)
void units_print_backend(const unit_system *us)
cdef extern from "adiabatic_index.h":
float hydro_gamma
cdef extern from "physical_constants.h":
struct phys_const:
double const_boltzmann_k
double const_proton_mass
pass
void phys_const_init(unit_system* us, phys_const* internal_const);
cdef extern from "cooling_struct.h":
struct cooling_function_data:
pass
cdef extern from "cooling.h":
void cooling_init(const swift_params* parameter_file,
const unit_system* us,
const phys_const* phys_const,
cooling_function_data* cooling);
void cooling_print_backend(
const cooling_function_data* cooling)
void cooling_cool_part(
const phys_const* phys_const,
const unit_system* us,
const cooling_function_data* cooling,
part* p, xpart* xp, float dt)
float cooling_rate(
const phys_const* const phys_const, const unit_system* us,
const cooling_function_data* cooling, const part* p)
cimport pyswiftsim
from pyswiftsim cimport pointer
from pyswiftsim.pointer cimport Pointer
from pyswiftsim.pointer cimport pointer_check
from pyswiftsim.pointer cimport pointer_create
from pyswiftsim cimport part_def
from pyswiftsim.part_def cimport hydro_space
from pyswiftsim cimport equation_of_state
from pyswiftsim.equation_of_state cimport gas_entropy_from_internal_energy
from pyswiftsim.equation_of_state cimport gas_soundspeed_from_internal_energy
from libc.stdlib cimport malloc, free
from libc.math cimport M_PI
cimport numpy as np
import numpy as np
ctypedef np.float_t DTYPE_t
_hydro_gamma = hydro_gamma
def read_params(str filename):
# transform python to C
tmp = filename.encode(u"ascii")
cdef char *cfilename = tmp
if cfilename == NULL:
raise Exception("Unable to convert filename to char*")
# init params
cdef swift_params *params = <swift_params*> malloc(sizeof(swift_params));
parser_read_file(cfilename, params)
#parser_print_params(params)
return pointer_create(
params,
pointer.type_swift_params
)
def init_units(Pointer p_params):
# deal with inputs
pointer_check(p_params, pointer.type_swift_params)
cdef swift_params *params = <swift_params*> p_params._data
# init units
cdef unit_system *units = <unit_system*> malloc(sizeof(unit_system));
cdef const char* category = "InternalUnitSystem"
units_init(units, params, category)
#units_print_backend(units)
# init constants
cdef phys_const *constants = <phys_const*> malloc(sizeof(phys_const));
phys_const_init(units, constants)
# return
d = {}
d["constants"] = pointer_create(constants, pointer.type_phys_const)
d["units"] = pointer_create(units, pointer.type_unit_system)
return d
def pycooling_init(Pointer p_params,
Pointer p_units,
Pointer p_constants):
# deal with inputs
pointer_check(p_params, pointer.type_swift_params)
cdef swift_params *params = <swift_params*> p_params._data
pointer_check(p_units, pointer.type_unit_system)
cdef unit_system *units = <unit_system*> p_units._data
pointer_check(p_constants, pointer.type_phys_const)
cdef phys_const *constants = <phys_const*> p_constants._data
#units_print_backend(units)
# init cooling
cdef cooling_function_data *cooling = <cooling_function_data*> malloc(sizeof(cooling_function_data));
cooling_init(params, units, constants, cooling)
# print results
#cooling_print_backend(cooling)
# store and return results
return pointer_create(cooling, pointer.type_cooling_function_data)
def pycooling_rate(Pointer p_units,
Pointer p_cooling,
Pointer p_constants,
np.ndarray [DTYPE_t, ndim=1] rho,
np.ndarray[DTYPE_t, ndim=1] u):
if rho.shape[0] != u.shape[0]:
raise Exception("Rho and u should be of the same size!")
cdef int N = rho.shape[0]
cdef np.ndarray[DTYPE_t, ndim=1] rate = np.empty_like(rho)
# deal with inputs
# p_units
pointer_check(p_units, pointer.type_unit_system)
cdef unit_system *units = <unit_system*> p_units._data
# p_cooling
pointer_check(p_cooling, pointer.type_cooling_function_data)
cdef cooling_function_data *cooling = <cooling_function_data*> p_cooling._data
# p_constants
pointer_check(p_constants, pointer.type_phys_const)
cdef phys_const *constants = <phys_const*> p_constants._data
# initialize particles
cdef part *p = <part*> malloc(sizeof(part) * N)
cdef xpart *xp = <xpart*> malloc(sizeof(xpart) * N)
cdef hydro_space hs # shadowfax stuff => no need to initialize
for i in range(N):
# init
part_def.hydro_init_part(&p[i], &hs)
p[i].id = i
p[i].mass = 100.
p[i].h = 1.
p[i].rho = rho[i]
p[i].entropy = gas_entropy_from_internal_energy(rho[i], u[i])
for j in range(3):
p[i].x[j] = 0
p[i].v[j] = 0
# compute cooling
for i in range(N):
# compute cooling
rate[i] = cooling_rate(constants, units, cooling, &p[i])
return rate
def soundspeed_from_internal_energy(np.ndarray [DTYPE_t, ndim=1] rho,
np.ndarray [DTYPE_t, ndim=1] u):
if rho.shape[0] != u.shape[0]:
raise Exception("Rho and u should be of the same size!")
cdef int N = rho.shape[0]
cdef np.ndarray[DTYPE_t, ndim=1] cs = np.empty_like(rho)
for i in range(N):
cs[i] = gas_soundspeed_from_internal_energy(rho[i], u[i])
return cs
cdef extern from "equation_of_state.h":
cdef float gas_entropy_from_internal_energy(float density, float u)
cdef float gas_soundspeed_from_entropy(float density, float entropy)
cdef float gas_soundspeed_from_internal_energy(float density, float u)
cdef extern from "part.h":
struct part:
float h
float mass
float rho
long long id
double x[3]
float v[3]
float entropy
struct xpart:
pass
cdef extern from "hydro.h":
void hydro_init_part(
part *p, const hydro_space *hs)
cdef extern from "hydro_space.h":
struct hydro_space:
pass
cdef extern from "debug.h":
void printParticle(const part *parts, const xpart *xparts,
long long int id, size_t N);
cdef extern from "clocks.h":
void clocks_set_cpufreq(unsigned long long freq);
cdef enum pointer_type:
type_swift_params
type_unit_system
type_phys_const
type_cooling_function_data
type_count
cdef const char **pointer_name
cdef class Pointer:
cdef void *_data
cdef int _data_type
cdef _setup(self, void *params, int data_type)
cdef pointer_create(void *params, int data_type)
cdef pointer_check(Pointer ptr, int data_type)
from libc.stdlib cimport free
from cooling cimport unit_system, phys_const
clocks_set_cpufreq(0); # estimate automatically cpufreq and init time
pointer_name = [
"swift_params",
"unit_system",
"phys_const",
"cooling_function_data",
"count"
]
# wrapper around C pointers
cdef class Pointer:
def __cinit__(self):
self._data = NULL
cdef _setup(self, void *params, int data_type):
self._data = params
self._data_type = data_type
return self
def __str__(self):
txt = "%s" % pointer_name[self._data_type]
return txt
def __del__(self):
free(self._data)
def __getattr__(self, value):
data_type = self._data_type
if data_type == type_unit_system:
return pointer_unit_system_getattr(self, value)
if data_type == type_phys_const:
return pointer_phys_const_getattr(self, value)
else:
raise Exception("Data type not implemented")
# method to call when creating a pointer
cdef pointer_create(void *params, int data_type):
return Pointer()._setup(params, data_type)
# method to call when receiving a pointer
cdef pointer_check(Pointer ptr, int data_type):
if (ptr._data_type != data_type):
exp = pointer_name[data_type]
rec = pointer_name[ptr._data_type]
raise Exception("Expecting pointer of type %s, got %s" % (exp, rec))
if (ptr._data == NULL):
data_name = pointer_name[ptr._data_type]
raise Exception("Value not set for pointer type %s" % data_name)
cdef pointer_unit_system_getattr(Pointer ptr, str attr):
pointer_check(ptr, type_unit_system)
cdef unit_system us = (<unit_system*> ptr._data)[0]
if attr == "UnitMass_in_cgs":
return us.UnitMass_in_cgs
elif attr == "UnitLength_in_cgs":
return us.UnitLength_in_cgs
elif attr == "UnitTime_in_cgs":
return us.UnitTime_in_cgs
elif attr == "UnitCurrent_in_cgs":
return us.UnitCurrent_in_cgs
elif attr == "UnitTemperature_in_cgs":
return us.UnitTemperature_in_cgs
else:
raise Exception("unit_system does not contain %s" % attr)
cdef pointer_phys_const_getattr(Pointer ptr, str attr):
pointer_check(ptr, type_phys_const)
cdef phys_const consts = (<phys_const*> ptr._data)[0]
if attr == "const_boltzmann_k":
return consts.const_boltzmann_k
if attr == "const_proton_mass":
return consts.const_proton_mass
else:
raise Exception("phys_const does not contain %s" % attr)
#!/usr/bin/env python3
__author__ = "loic hausammann"
__copyright__ = "GPLv3"
descr = "Swift is a cosmological hydrodynamical code: SPH With Inter-dependent Fine-grained Tasking"
# import
from setuptools import setup, find_packages, Extension
from Cython.Build import cythonize
import numpy
import os
# read makefile in order to get libraries
makefile = "../src/Makefile"
# save environment
old = {}
def saveEnviron(key, old):
""" Save environment variable before modifying them
Key is the environment variable and old is a dictionary
containing all the modified values.
"""
old[key] = ""
if key in os.environ:
old[key] = os.environ[key]
return old
old = saveEnviron("CC", old)
old = saveEnviron("LDFLAGS", old)
def getValueFromMakefile(f, key):
"""
Read a file and return the value of a variable.
f is an open file and value if the requested key.
The variable should be matching the following patern:
KEY = VALUE
"""
test = key + " = "
ret = None
for line in f:
start = line[:len(test)]
if test == start:
ret = line[len(test):]
f.seek(0)
if ret is None:
raise Exception("Unable to find key %s" % key)
else:
# remove last character "\n"
return ret[:-1]
with open(makefile, "r") as f:
hdf5_root = getValueFromMakefile(f, "H5CC")
cc = getValueFromMakefile(f, "CC")
# python requirement
install_requires = []
# include
include = [
numpy.get_include(),
"../src",
hdf5_root + "/include"
]
# libraries
lib = ["m", "hdf5", "swiftsim"]
# Extension object required by setup
ext = []
# pointer
tmp = Extension("pyswiftsim.pointer",
["pyswiftsim/pointer.pyx"],
include_dirs=include,
libraries=lib)
tmp = cythonize(tmp)
ext.extend(tmp)
# cooling wrapper
tmp = Extension("pyswiftsim.cooling",
["pyswiftsim/cooling.pyx"],
include_dirs=include,
libraries=lib)
tmp = cythonize(tmp)
ext.extend(tmp)
# scripts
scripts = []
# data to copy
data = []
# set environment variables, please save the environment variable before modifying them
# path to libswiftsim
ldflags = "LDFLAGS"
if ldflags not in os.environ:
os.environ[ldflags] = ""
os.environ[ldflags] += " -L" + hdf5_root + "/lib"
# compiler
os.environ["CC"] = cc
setup(
name = "pyswiftsim",
version = "0.6.0",
author = "Hausammann Loic",
author_email = "loic.hausammann@epfl.ch",
description = descr,
license = "GPLv3",
keywords = "swift hpc cosmology",
url = "https://gitlab.cosma.dur.ac.uk/swift/swiftsim",
packages = find_packages(),
data_files = data,
scripts = scripts,
install_requires = install_requires,
ext_modules = ext,
)
def restoreEnviron(old):
for key, value in old.items():
os.environ[key] = value
restoreEnviron(old)
#!/usr/bin/env python3
from pyswiftsim import pointer
from pyswiftsim import cooling
import astropy.units
import matplotlib.pyplot as plt
import numpy as np
# number of particles
N = 100
# params file
params = "/home/loikki/swift_test/cooling_box/coolingBox.yml"
# read params
params = cooling.read_params(params)
#print(params)
# init units
d = cooling.init_units(params)
units = d["units"]
consts = d["constants"]
# init cooling
cooling_data = cooling.pycooling_init(params, units, consts)
# Init variables
rho = np.logspace(-6, 4, N).astype(np.float) # in cm-3
rho *= consts.const_proton_mass * units.UnitLength_in_cgs**3
T = np.logspace(1,9, N).astype(np.float) # in K
rho, T = np.meshgrid(rho, T)
shape = rho.shape
rho = rho.reshape(np.prod(shape))
T = T.reshape(np.prod(shape))
# u = k_b T / (gamma - 1) mu m_p
mu = 2
u = consts.const_boltzmann_k * T / ((cooling._hydro_gamma - 1.) * mu * consts.const_proton_mass)
# compute cooling
rate = cooling.pycooling_rate(units, cooling_data, consts, rho, u)
cs = cooling.soundspeed_from_internal_energy(rho, u)
rho = rho.reshape(shape) / (units.UnitLength_in_cgs**3 * consts.const_proton_mass)
T = T.reshape(shape)
# du / dt [e / s] = [kg m2 / s3]
rate = rate.reshape(shape) * units.UnitMass_in_cgs * units.UnitLength_in_cgs**2 / units.UnitTime_in_cgs**3
cs = cs.reshape(shape) * units.UnitLength_in_cgs / units.UnitTime_in_cgs
L_cool = -rate * cs / astropy.units.kpc.to("cm")
plt.title("Cooling Length")
plt.contourf(rho, T, np.log(L_cool))
plt.xlabel("Density [$cm^{-3}$]")
plt.ylabel("Temperature [K]")
plt.colorbar()
ax = plt.gca()
ax.set_xscale('log')
ax.set_yscale('log')
plt.show()
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment