Skip to content
Snippets Groups Projects
Commit 6d8f7237 authored by lhausamm's avatar lhausamm
Browse files

Add cosmology

parent f78bb16c
No related branches found
No related tags found
1 merge request!5Last commit before huge rework
......@@ -45,6 +45,8 @@ GrackleCooling:
ProvideSpecificHeatingRates: 0 # User provide specific heating rates
SelfShieldingMethod: 0 # Grackle (<= 3) or Gear self shielding method
OutputMode: 0 # Write in output corresponding primordial chemistry mode
MaxSteps: 1000
ConvergenceLimit: 1e-2
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
......
......@@ -30,7 +30,7 @@ plt.rc('text', usetex=True)
# PARAMETERS
# grackle primordial chemistry (for comparison)
primordial_chemistry = 1
primordial_chemistry = 3
# reference data
grackle_filename = "grackle.hdf5"
......@@ -213,6 +213,9 @@ def initialize_swift(filename):
d["us"] = us
d["pconst"] = pconst
cosmo = wrapper.cosmologyInit(params, us, pconst)
d["cosmo"] = cosmo
# init cooling
cooling = wrapper.coolingInit(params, us, pconst)
d["cooling"] = cooling
......@@ -281,6 +284,7 @@ if __name__ == "__main__":
pconst = d["pconst"]
us = d["us"]
params = d["params"]
cosmo = d["cosmo"]
cooling = d["cooling"]
chemistry = d["chemistry"]
......@@ -304,7 +308,9 @@ if __name__ == "__main__":
for i in range(N-1):
ene = np.array([energy[i]], dtype=np.float32)
if not compute_equilibrium:
energy[i+1] = wrapper.doCooling(pconst, us, cooling,
energy[i+1] = wrapper.doCooling(pconst, us,
cosmo,
cooling,
chemistry,
rho,
ene,
......@@ -315,7 +321,9 @@ if __name__ == "__main__":
frac = np.zeros([1, N])
for j in range(N):
frac[:, j] = d[fields[j]]
energy[i+1] = wrapper.doCooling(pconst, us, cooling,
energy[i+1] = wrapper.doCooling(pconst, us,
cosmo,
cooling,
chemistry,
rho,
ene,
......
......@@ -31,7 +31,7 @@ plt.rc('text', usetex=True)
# PARAMETERS
# grackle primordial chemistry (for comparison)
primordial_chemistry = 1
primordial_chemistry = 3
# reference data
grackle_filename = "grackle.hdf5"
......@@ -228,6 +228,10 @@ def initialize_swift(filename):
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
......@@ -382,6 +386,7 @@ if __name__ == "__main__":
pconst = d["pconst"]
us = d["us"]
params = d["params"]
cosmo = d["cosmo"]
cooling = d["cooling"]
chemistry = d["chemistry"]
......@@ -397,7 +402,9 @@ if __name__ == "__main__":
rate = np.zeros(d["density"].shape)
if compute_equilibrium:
rate = wrapper.coolingRate(pconst, us, cooling,
rate = wrapper.coolingRate(pconst, us,
cosmo,
cooling,
chemistry,
d["density"].astype(np.float32),
d["energy"].astype(np.float32),
......@@ -409,7 +416,9 @@ if __name__ == "__main__":
for i in range(N):
frac[:, i] = d[fields[i]]
rate = wrapper.coolingRate(pconst, us, cooling,
rate = wrapper.coolingRate(pconst, us,
cosmo,
cooling,
chemistry,
d["density"].astype(np.float32),
d["energy"].astype(np.float32),
......
......@@ -17,12 +17,15 @@ You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
"""
from sys import argv
symbol_table = {
" int ": "i",
" double ": "d",
" char ": "c",
}
def _generateFormatString(filename):
"""
Open a file containing a single struct and generate a format string
......@@ -31,7 +34,7 @@ def _generateFormatString(filename):
struct_format = ""
struct_attribute = []
with open(filename, "r") as f:
line = ""
# skip until reaching struct
......@@ -63,21 +66,20 @@ def _generateFormatString(filename):
j = line.index(";")
struct_attribute.append(line[i:j])
return struct_format, struct_attribute
if __name__ == "__main__":
filename = "chemistry_data.h"
filename = argv[-1]
struct_format, struct_attribute = _generateFormatString(filename)
attr = "[\n"
for i in struct_attribute:
attr += "\t'%s',\n" % i
attr += "]"
txt = """
_format = "{form}"
_name = {attr}
......
......@@ -495,3 +495,60 @@ class ChemistryGlobalData(SwiftStruct):
def __init__(self, data, parent=None):
super().__init__(self.struct_format, data, parent)
######################################################################
# #
# PhysConst #
# #
######################################################################
class Cosmology(SwiftStruct):
_format = "ddddddddddddddddddddddddddddddddddddppppdd"
_name = [
'a',
'a_inv',
'a2_inv',
'a3_inv',
'a_factor_internal_energy',
'a_factor_pressure',
'a_factor_sound_speed',
'a_factor_mu',
'a_factor_grav_accel',
'a_factor_hydro_accel',
'z',
'H',
'critical_density',
'time_step_factor',
'a_dot',
'time',
'lookback_time',
'w',
'a_begin',
'a_end',
'time_begin',
'time_end',
'time_base',
'time_base_inv',
'h',
'H0',
'Hubble_time',
'Omega_m',
'Omega_b',
'Omega_lambda',
'Omega_r',
'Omega_k',
'w_0',
'w_a',
'log_a_begin',
'log_a_end',
'drift_fac_interp_table',
'grav_kick_fac_interp_table',
'hydro_kick_fac_interp_table',
'time_interp_table',
'time_interp_table_offset',
'universe_age_at_present_day',
]
def __init__(self, data, parent=None):
super().__init__(self.struct_format, data, parent)
......@@ -115,6 +115,7 @@ PyArrayObject* pycooling_rate(PyObject* self, PyObject* args) {
PyObject *pycooling;
PyObject *pyus;
PyObject *pycosmo;
PyObject *pypconst;
PyObject *pychemistry;
......@@ -126,9 +127,10 @@ PyArrayObject* pycooling_rate(PyObject* self, PyObject* args) {
/* parse argument */
if (!PyArg_ParseTuple(args,
"OOOOOO|fO",
"OOOOOOO|fO",
&pypconst,
&pyus,
&pycosmo,
&pycooling,
&pychemistry,
&rho,
......@@ -176,14 +178,13 @@ PyArrayObject* pycooling_rate(PyObject* self, PyObject* args) {
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 *rate = PyArray_SimpleNew(PyArray_NDIM(energy), PyArray_DIMS(energy), NPY_FLOAT);
......@@ -202,16 +203,16 @@ PyArrayObject* pycooling_rate(PyObject* self, PyObject* args) {
if (fractions != NULL)
pycooling_set_fractions(&xp, fractions, i);
else
cooling_first_init_part(&p, &xp, cooling);
cooling_first_init_part(pconst, us, cosmo, cooling, &p, &xp);
chemistry_first_init_part(&p, &xp, chemistry);
chemistry_first_init_part(pconst, us, cosmo, chemistry, &p, &xp);
/* compute cooling rate */
float *tmp = PyArray_GETPTR1(rate, i);
#ifdef COOLING_GRACKLE
*tmp = cooling_rate(pconst, us, cooling, &p, &xp, dt);
*tmp = cooling_rate(pconst, us, cosmo, cooling, &p, &xp, dt);
#else
*tmp = cooling_rate(pconst, us, cooling, &p, &xp);
*tmp = cooling_rate(pconst, us, cosmo, cooling, &p, &xp);
#endif
}
......@@ -227,7 +228,7 @@ PyArrayObject* pycooling_rate(PyObject* self, PyObject* args) {
* @brief Compute cooling
*
* args is expecting pyswiftsim classes in the following order:
* PhysConst, UnitSystem and CoolingFunctionData.
* PhysConst, UnitSystem, Cosmology and CoolingFunctionData.
* Then two numpy arrays (density and specific energy) and a
* float for the time step
*
......@@ -240,6 +241,7 @@ PyArrayObject* pycooling_do_cooling(PyObject* self, PyObject* args) {
PyObject *pycooling;
PyObject *pyus;
PyObject *pycosmo;
PyObject *pypconst;
PyObject *pychemistry;
......@@ -251,9 +253,10 @@ PyArrayObject* pycooling_do_cooling(PyObject* self, PyObject* args) {
/* parse argument */
if (!PyArg_ParseTuple(args,
"OOOOOOf|O",
"OOOOOOOf|O",
&pypconst,
&pyus,
&pycosmo,
&pycooling,
&pychemistry,
&rho,
......@@ -301,6 +304,10 @@ PyArrayObject* pycooling_do_cooling(PyObject* self, PyObject* args) {
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;
......@@ -327,15 +334,15 @@ PyArrayObject* pycooling_do_cooling(PyObject* self, PyObject* args) {
/* initialize fields */
hydro_set_internal_energy_dt(&p, 0.);
chemistry_first_init_part(&p, &xp, chemistry);
chemistry_first_init_part(pconst, us, cosmo, chemistry, &p, &xp);
if (fractions != NULL)
pycooling_set_fractions(&xp, fractions, i);
else
cooling_first_init_part(&p, &xp, cooling);
cooling_first_init_part(pconst, us, cosmo, cooling, &p, &xp);
/* compute cooling */
cooling_cool_part(pconst, us, cooling, &p, &xp, dt);
cooling_cool_part(pconst, us, cosmo, cooling, &p, &xp, dt);
/* update energy */
float *tmp = PyArray_GETPTR1(new_energy, i);
......
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#include "pyswiftsim_tools.h"
#include "cosmology_wrapper.h"
/**
* @brief Initialize the chemistry
*
* args is expecting pyswiftsim classes in the following order:
* SwiftParams, UnitSystem and PhysConst.
*
* @param self calling object
* @param args arguments
* @return Cosmology
*/
PyObject* pycosmology_init(PyObject* self, PyObject* args) {
PyObject *pyparams;
PyObject *pyus;
PyObject *pypconst;
int with_cosmo = 0;
/* parse arguments */
if (!PyArg_ParseTuple(args, "OOO|i", &pyparams, &pyus, &pypconst, &with_cosmo))
return NULL;
struct swift_params *params = (struct swift_params*) pytools_construct(pyparams, class_swift_params);
if (params == NULL)
return NULL;
struct unit_system *us = (struct unit_system*) pytools_construct(pyus, class_unit_system);
if (us == NULL)
return NULL;
struct phys_const *pconst = (struct phys_const*) pytools_construct(pypconst, class_phys_const);
if (pconst == NULL)
return NULL;
struct cosmology cosmo;
/* init cooling */
if (with_cosmo)
cosmology_init(params, us, pconst, &cosmo);
else
cosmology_init_no_cosmo(&cosmo);
/* construct python object */
PyObject *pycosmo = pytools_return(&cosmo, class_cosmology);
return pycosmo;
}
/*******************************************************************************
* This file is part of PYSWIFTSIM.
* Copyright (c) 2018 loic hausammann (loic.hausammann@epfl.ch)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef __PYSWIFTSIM_COSMOLOGY_H__
#define __PYSWIFTSIM_COSMOLOGY_H__
#include "pyswiftsim_tools.h"
PyObject* pycosmology_init(PyObject* self, PyObject* args);
#endif // __PYSWIFTSIM_COSMOLOGY_H__
......@@ -25,7 +25,8 @@ const size_t class_size[class_count] = {
sizeof(struct swift_params),
sizeof(struct phys_const),
sizeof(struct cooling_function_data),
sizeof(struct chemistry_global_data)
sizeof(struct chemistry_global_data),
sizeof(struct cosmology)
};
char *class_name[class_count] = {
......@@ -34,7 +35,8 @@ char *class_name[class_count] = {
"SwiftParams",
"PhysConst",
"CoolingFunctionData",
"ChemistryGlobalData"
"ChemistryGlobalData",
"Cosmology"
};
......
......@@ -68,6 +68,7 @@ enum swift_class {
class_phys_const,
class_cooling_function_data,
class_chemistry_global_data,
class_cosmology,
class_count /* should always be last! */
};
......
......@@ -45,9 +45,10 @@ PyObject* pyunit_system_test_struct(PyObject *self, PyObject *args)
PyObject* pyunit_system_init(PyObject *self, PyObject *args)
{
PyObject* parser;
char* section = "InternalUnitSystem";
/* parse arguement */
if (!PyArg_ParseTuple(args, "O", &parser))
if (!PyArg_ParseTuple(args, "O|s", &parser, &section))
return NULL;
struct swift_params *params = (struct swift_params*) pytools_construct(parser, class_swift_params);
......@@ -57,11 +58,12 @@ PyObject* pyunit_system_init(PyObject *self, PyObject *args)
/* initialize units */
struct unit_system us;
units_init(&us, params, "InternalUnitSystem");
units_init(&us, params, section);
/* initialize phys_const */
struct phys_const pconst;
phys_const_init(&us, &pconst);
phys_const_init(&us, params, &pconst);
/* create python object to return */
PyObject *pyus = pytools_return(&us, class_unit_system);
......
......@@ -20,6 +20,7 @@
#include "part_wrapper.h"
#include "parser_wrapper.h"
#include "cooling_wrapper.h"
#include "cosmology_wrapper.h"
#include "chemistry_wrapper.h"
#include "pyswiftsim_tools.h"
#include "config_wrapper.h"
......@@ -45,6 +46,9 @@ static PyMethodDef wrapper_methods[] = {
{"chemistryInit", pychemistry_init, METH_VARARGS,
"Initialize chemistry."},
{"cosmologyInit", pycosmology_init, METH_VARARGS,
"Initialize cosmology."},
{"coolingRate", pycooling_rate, METH_VARARGS,
"Compute the cooling rate.\n\n"
"Parameters\n"
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment