diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
new file mode 100644
index 0000000000000000000000000000000000000000..6251ddd24e369f246e64c15b92c399bcf5dca2a0
--- /dev/null
+++ b/.gitlab-ci.yml
@@ -0,0 +1,11 @@
+before_script:
+ # compile swift
+ - git clone https://gitlab.cosma.dur.ac.uk/swift/swiftsim.git
+ - cd swiftsim
+ - ./autogen
+ - ./configure --with-subgrid=GEAR
+ - make -j 4
+
+test_build:
+ script:
+ - python setup.py build
\ No newline at end of file
diff --git a/examples/cooling_box/CloudyData_UVB=HM2012.h5 b/examples/cooling_box/CloudyData_UVB=HM2012.h5
deleted file mode 100644
index 1e35a914e33cb6465b48babba8730ba47a92f6f2..0000000000000000000000000000000000000000
Binary files a/examples/cooling_box/CloudyData_UVB=HM2012.h5 and /dev/null differ
diff --git a/examples/cooling_box/cooling.yml b/examples/cooling_box/cooling.yml
deleted file mode 100644
index cbcb4c555033d6204cd1fbffcd1a70e613a3d47c..0000000000000000000000000000000000000000
--- a/examples/cooling_box/cooling.yml
+++ /dev/null
@@ -1,60 +0,0 @@
-# Define the system of units to use internally.
-InternalUnitSystem:
- UnitMass_in_cgs: 2.e33 # Grams
- UnitLength_in_cgs: 3.085e21 # Centimeters
- UnitVelocity_in_cgs: 1e5 # Centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
-
-# Parameters governing the time integration
-TimeIntegration:
- time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 5e-2 # The end time of the simulation (in internal units).
- dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
- dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
-
-# Parameters governing the snapshots
-Snapshots:
- basename: sedov # Common part of the name of output files
- time_first: 0. # Time of the first output (in internal units)
- delta_time: 1e-2 # Time difference between consecutive outputs (in internal units)
-
-# Parameters governing the conserved quantities statistics
-Statistics:
- delta_time: 1e-3 # Time between statistics output
- tasks_per_cell: 1000
-
-# Parameters for the hydrodynamics scheme
-SPH:
- resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
- CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
-
-# Parameters related to the initial conditions
-InitialConditions:
- file_name: ./sedov.hdf5 # The file to read
- h_scaling: 3.33
-
-
-# Cooling with Grackle 2.0
-GrackleCooling:
- CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table
- WithUVbackground: 1 # Enable or not the UV background
- Redshift: 0 # Redshift to use (-1 means time based redshift)
- WithMetalCooling: 1 # Enable or not the metal cooling
- ProvideVolumetricHeatingRates: 0 # User provide volumetric heating rates
- 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.
- theta: 0.7 # Opening angle (Multipole acceptance criterion)
- epsilon: 0.1 # Softening length (in internal units).
- a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
- r_cut_max: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
- r_cut_min: 0.1 # (Optional) Cut-off in number of top-level cells below which no truncation of FMM forces are performed (this is the default value).
-
-GearChemistry:
- InitialMetallicity: 0.01295
\ No newline at end of file
diff --git a/examples/cooling_box/generate_grackle_data.py b/examples/cooling_box/generate_grackle_data.py
deleted file mode 100644
index e9ed1bb77a4e05a2534773fdbfb8d40d51f8e152..0000000000000000000000000000000000000000
--- a/examples/cooling_box/generate_grackle_data.py
+++ /dev/null
@@ -1,262 +0,0 @@
-#!/usr/bin/env python
-########################################################################
-#
-# Cooling rate data generation. This code is a modified version
-# of a code developed by the Grackle Team.
-#
-#
-# Copyright (c) 2013-2016, Grackle Development Team.
-#
-# Distributed under the terms of the Enzo Public Licence.
-#
-# The full license is in the file LICENSE, distributed with this
-# software.
-########################################################################
-
-"""
-Generate (or update) a hdf5 file containing
-the grackle data for each primordial chemistry
-"""
-
-import numpy as np
-from h5py import File
-
-from pygrackle import \
- chemistry_data, \
- setup_fluid_container
-
-from pygrackle.utilities.physical_constants import \
- mass_hydrogen_cgs
-
-filename = "grackle.hdf5"
-
-# metal_mass_fraction = 1e-10
-metal_mass_fraction = 0.01295
-
-tolerance = 1e-1
-max_iter = 1e8
-
-t_end = 1. # in cooling time
-N_time = 1000
-
-part_density = 1.0 # in atom/cc
-part_temperature = 1e6 # in K
-
-
-def generate_data(primordial_chemistry):
- """
- Compute the cooling a single particles and save its evolution
- to an hdf5 file
- """
- current_redshift = 0.
-
- # Set solver parameters
- my_chemistry = chemistry_data()
- my_chemistry.use_grackle = 1
- my_chemistry.with_radiative_cooling = 1
- my_chemistry.primordial_chemistry = primordial_chemistry
- my_chemistry.metal_cooling = 1
- my_chemistry.UVbackground = 1
- my_chemistry.self_shielding_method = 0
- my_chemistry.H2_self_shielding = 0
- my_chemistry.grackle_data_file = "CloudyData_UVB=HM2012.h5"
-
- # Set units
- my_chemistry.comoving_coordinates = 0 # proper units
- my_chemistry.a_units = 1.0
- my_chemistry.a_value = 1.0 / (1.0 + current_redshift) / \
- my_chemistry.a_units
- my_chemistry.length_units = 3.085e21
- my_chemistry.density_units = 2.0e33 / my_chemistry.length_units**3
- my_chemistry.velocity_units = 1e5
- my_chemistry.time_units = my_chemistry.length_units / \
- my_chemistry.velocity_units
-
- density = np.array([part_density * mass_hydrogen_cgs])
- # Call convenience function for setting up a fluid container.
- # This container holds the solver parameters, units, and fields.
- temperature = np.array([part_temperature])
- fc = setup_fluid_container(my_chemistry,
- metal_mass_fraction=metal_mass_fraction,
- temperature=temperature,
- density=density,
- converge=True,
- tolerance=tolerance,
- max_iterations=max_iter)
-
- # compute time steps
- fc.calculate_cooling_time()
- t_max = t_end * abs(fc["cooling_time"])
-
- t = np.linspace(0, t_max, N_time)
-
- f, data = write_input(
- fc, my_chemistry, primordial_chemistry)
-
- dt = t[1] - t[0]
-
- # compute evolution
- E = np.zeros(t.shape)
- E[0] = fc["energy"]
- for i in range(N_time-1):
- fc.solve_chemistry(dt)
- E[i+1] = fc["energy"]
-
- write_output(f, data, fc, t, E)
-
-
-def write_output(f, data, fc, t, E):
- """
- Write the output of the simulation
- """
- data.create_group("Output")
- tmp = data["Output"]
-
- # Energy
- dset = tmp.create_dataset("Energy", E.shape,
- dtype=E.dtype)
- dset[:] = E
-
- # Time
- dset = tmp.create_dataset("Time", t.shape,
- dtype=t.dtype)
- dset[:] = t
-
- f.close()
-
-
-def simulationAreEqual(fields, data):
- """
- Check if current simulation has already
- been run
- """
- test = True
- for key in fields.keys():
- if data.attrs[key] != fields[key]:
- test = False
- break
-
- return test
-
-
-def getSameSimulation(data, density, my_chemistry):
- """
- Get indice of the same simulation.
- returns -1 if not found
- """
- fields = {
- "MetalCooling": my_chemistry.metal_cooling,
- "UVBackground": my_chemistry.UVbackground,
- "SelfShieldingMethod": my_chemistry.self_shielding_method,
- "MetalMassFraction": metal_mass_fraction,
- "ScaleFactor": my_chemistry.a_value,
- "Temperature": part_temperature,
- "Density": density,
- }
- ind = 0
- for d in data:
- if "Params" in data[d]:
- tmp = data[d]["Params"]
- if simulationAreEqual(fields, tmp):
- return ind
- ind += 1
-
- return -1
-
-
-def write_input(fc, my_chemistry, primordial_chemistry):
- """
- Write the simulation input in the HDF5 file
- """
- f = File(filename, "a")
- data_name = "PrimordialChemistry%i" % primordial_chemistry
- if data_name in f:
- data = f[data_name]
- else:
- print("Creating Dataset")
- data = f.create_group(data_name)
-
- i = getSameSimulation(data, fc["density"], my_chemistry)
- if i != -1:
- print("Updating Simulation")
- data.pop(data.keys()[i])
- data = data.create_group("Simulation %i" % len(data))
-
- # Units
- data.create_group("Units")
- tmp = data["Units"]
- tmp.attrs["Length"] = my_chemistry.length_units
- tmp.attrs["Density"] = my_chemistry.density_units
- tmp.attrs["Velocity"] = my_chemistry.velocity_units
- tmp.attrs["Time"] = my_chemistry.time_units
-
- # Parameters
- data.create_group("Params")
- tmp = data["Params"]
- tmp.attrs["MetalCooling"] = my_chemistry.metal_cooling
- tmp.attrs["UVBackground"] = my_chemistry.UVbackground
- tmp.attrs["SelfShieldingMethod"] = my_chemistry.self_shielding_method
- tmp.attrs["MetalMassFraction"] = metal_mass_fraction
- tmp.attrs["ScaleFactor"] = my_chemistry.a_value
- tmp.attrs["Tolerance"] = tolerance
- tmp.attrs["MaxIter"] = max_iter
- tmp.attrs["Density"] = fc["density"]
- tmp.attrs["Temperature"] = part_temperature
- tmp.attrs["Energy"] = fc["energy"][0]
-
- # Inputs
- data.create_group("Input")
- tmp = data["Input"]
-
- write_fractions(tmp, fc, primordial_chemistry)
-
- return f, data
-
-
-def write_fractions(tmp, fc, primordial_chemistry):
- """
- Write the element fractions in the HDF5 file
- """
- fields = get_fields(primordial_chemistry)
- for i in fields:
- dset = tmp.create_dataset(i, fc[i].shape,
- dtype=fc[i].dtype)
- dset[:] = fc[i] / fc["density"]
-
-
-def get_fields(primordial_chemistry):
- """
- Return the fields present in grackle
- """
- fields = [
- "metal"
- ]
- if primordial_chemistry > 0:
- fields.extend([
- "HI",
- "HII",
- "HeI",
- "HeII",
- "HeIII",
- "de"])
-
- if primordial_chemistry > 1:
- fields.extend([
- "HM",
- "H2I",
- "H2II"
- ])
-
- if primordial_chemistry > 2:
- fields.extend([
- "DI",
- "DII",
- "HDI"
- ])
- return fields
-
-
-if __name__ == "__main__":
- for i in range(4):
- print("Computing Primordial Chemistry %i" % i)
- generate_data(i)
diff --git a/examples/cooling_box/getCoolingTable.sh b/examples/cooling_box/getCoolingTable.sh
deleted file mode 100644
index 809958bfc236e5ab0f7be924c62789fa13b3ac29..0000000000000000000000000000000000000000
--- a/examples/cooling_box/getCoolingTable.sh
+++ /dev/null
@@ -1,2 +0,0 @@
-#!/bin/bash
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5
diff --git a/examples/cooling_box/plot_cooling.py b/examples/cooling_box/plot_cooling.py
deleted file mode 100644
index 917ea82d48a324c1a5fbceb076b89f865129a5dd..0000000000000000000000000000000000000000
--- a/examples/cooling_box/plot_cooling.py
+++ /dev/null
@@ -1,335 +0,0 @@
-#!/usr/bin/env python3
-"""
-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/>.
-"""
-
-from pyswiftsim import wrapper
-
-import numpy as np
-import matplotlib.pyplot as plt
-from os.path import isfile
-from h5py import File
-from glob import glob
-
-plt.rc('text', usetex=True)
-
-# PARAMETERS
-
-# grackle primordial chemistry (for comparison)
-primordial_chemistry = 0
-
-# reference data
-grackle_filename = "grackle.hdf5"
-compute_equilibrium = False
-
-# swift params filename
-filename = "cooling.yml"
-
-# particle data
-# in atom/cc if generating
-# if reading => same than in grackle.hdf5
-part_density = 2.4570967948472597e7
-part_temperature = 1e6 # in K
-gamma = 5./3.
-
-# default parameters
-t_end = 1.
-N_time = 100000
-
-# cooling output
-swift_output_filename = "swift.hdf5"
-
-# simulations
-simulations = {
- "SWIFT": "/home/loikki/data/swiftsim/examples/CoolingBox/coolingBox_*.hdf5"
-}
-# SCRIPT
-
-
-def get_fields(primordial_chemistry):
- fields = [
- "metal"
- ]
- if primordial_chemistry > 0:
- fields.extend([
- "HI",
- "HII",
- "HeI",
- "HeII",
- "HeIII",
- "de"])
-
- if primordial_chemistry > 1:
- fields.extend([
- "HM",
- "H2I",
- "H2II"
- ])
-
- if primordial_chemistry > 2:
- fields.extend([
- "DI",
- "DII",
- "HDI"
- ])
- return fields
-
-
-def generate_default_initial_condition(us, pconst):
- print("Generating default initial conditions")
- d = {}
- # generate grid
- d["temperature"] = part_temperature
-
- # Deal with units
- rho = part_density * us.UnitLength_in_cgs**3 * pconst.const_proton_mass
- d["density"] = [rho]
-
- energy = pconst.const_boltzmann_k * part_temperature \
- / us.UnitTemperature_in_cgs
- energy /= (gamma - 1.) * pconst.const_proton_mass
- d["energy"] = energy
-
- d["time"] = np.linspace(0, t_end, N_time)
-
- 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,
- "Temperature": part_temperature,
- "Density": part_density,
- }
-
- 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, pconst, primordial_chemistry,
- cooling, chemistry):
- 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]]
-
- # read units
- tmp = data["Units"].attrs
-
- u_len = tmp["Length"] / us.UnitLength_in_cgs
- u_time = tmp["Time"] / us.UnitTime_in_cgs
- u_den = tmp["Density"] * us.UnitLength_in_cgs**3 / us.UnitMass_in_cgs
-
- # read input
- tmp = data["Params"]
- energy = tmp.attrs["Energy"] * u_len**2 / u_time**2
- d["energy"] = energy
-
- density = tmp.attrs["Density"] * u_den
- d["density"] = density
-
- # # read fractions
- # tmp = data["Input"]
- # for i in get_fields(primordial_chemistry):
- # d[i] = tmp[i][:]
-
- # read output
- tmp = data["Output"]
-
- energy = tmp["Energy"][:] * u_len**2 / u_time**2
- d["out_energy"] = energy
-
- t = tmp["Time"][:] * u_time
- d["time"] = t
-
- f.close()
- 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
-
- 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(energy, data, us):
- print("Plotting solution")
- rho = data["density"]
- time = data["time"]
-
- ref = False
- if "out_energy" in data:
- ref = True
- grackle_energy = data["out_energy"]
-
- # change units => cgs
- rho *= us.UnitMass_in_cgs / us.UnitLength_in_cgs**3
-
- energy *= us.UnitLength_in_cgs**2 / us.UnitTime_in_cgs**2
-
- if ref:
- grackle_energy *= us.UnitLength_in_cgs**2 / us.UnitTime_in_cgs**2
-
- time *= us.UnitTime_in_cgs
-
- # do plot
-
- plt.figure()
- plt.plot(time, energy, label="PySWIFTsim, %s" % wrapper.configGetCooling())
- if ref:
- plt.plot(time, grackle_energy,
- label="Grackle, Prim. Chem. %i" % primordial_chemistry)
-
- for i in simulations:
- files = glob(simulations[i])
- files.sort()
- E = np.zeros([len(files)])
- t = np.zeros([len(files)])
- j = 0
- for f in files:
- f = File(f)
- tmp = f["Units"]
- UL = tmp.attrs["Unit length in cgs (U_L)"]
- UT = tmp.attrs["Unit time in cgs (U_t)"]
- UE = UL**2 / UT**2
- tmp = f["PartType0"]
- E[j] = np.mean(tmp["InternalEnergy"]) * UE
- t[j] = f["Header"].attrs["Time"] * UT
- print(t[j])
- j += 1
- cooling_type = f["SubgridScheme"].attrs["Cooling Model"]
- plt.plot(t, E, label=i + ": " + cooling_type.decode("utf8"))
- plt.xlabel("Time [s]")
- plt.ylabel("Energy [erg]")
- plt.legend()
- plt.show()
-
-
-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"]
-
- if isfile(grackle_filename) and False:
- d = read_grackle_data(grackle_filename, us, pconst,
- primordial_chemistry, cooling,
- chemistry)
- else:
- d = generate_default_initial_condition(us, pconst)
-
- t = d["time"]
- dt = t[1] - t[0]
-
- # du / dt
- print("Computing cooling...")
-
- N = len(t)
- energy = np.zeros(t.shape)
- energy[0] = d["energy"]
- rho = np.array(d["density"], dtype=np.float32)
- 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,
- cosmo,
- cooling,
- chemistry,
- rho,
- ene,
- dt)
- else:
- fields = get_fields(primordial_chemistry)
- N = len(fields)
- frac = np.zeros([1, N])
- for j in range(N):
- frac[:, j] = d[fields[j]]
- energy[i+1] = wrapper.doCooling(pconst, us,
- cosmo,
- cooling,
- chemistry,
- rho,
- ene,
- dt,
- frac.astype(np.float32))
-
- plot_solution(energy, d, us)
diff --git a/examples/cooling_box/run.sh b/examples/cooling_box/run.sh
deleted file mode 100644
index d2f01dae43b010839cddca2d11d20d474dd76f0f..0000000000000000000000000000000000000000
--- a/examples/cooling_box/run.sh
+++ /dev/null
@@ -1,17 +0,0 @@
-#!/bin/bash
-
-# Get the Grackle cooling table
-if [ ! -e CloudyData_UVB=HM2012.h5 ]
-then
- echo "Fetching the Cloudy tables required by Grackle..."
- ./getCoolingTable.sh
-fi
-
-# generate grackle data
-if [ ! -e grackle.hdf5 ]
-then
- echo "Generating Grackle data..."
- ./generate_grackle_data.py
-fi
-
-./plot_cooling.py
diff --git a/examples/cooling_rate/cooling.yml b/examples/cooling_rate/cooling.yml
deleted file mode 100644
index 10608d34b8e2fcdb2ad2db441e5b0ae655195541..0000000000000000000000000000000000000000
--- a/examples/cooling_rate/cooling.yml
+++ /dev/null
@@ -1,71 +0,0 @@
-# Define the system of units to use internally.
-InternalUnitSystem:
- UnitMass_in_cgs: 1.989e43 # Grams
- UnitLength_in_cgs: 3.085e21 # Centimeters
- UnitVelocity_in_cgs: 20725573.785998672 # Centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
-
-# Parameters governing the time integration
-TimeIntegration:
- time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 5e-2 # The end time of the simulation (in internal units).
- dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
- dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
-
-# Parameters governing the snapshots
-Snapshots:
- basename: sedov # Common part of the name of output files
- time_first: 0. # Time of the first output (in internal units)
- delta_time: 1e-2 # Time difference between consecutive outputs (in internal units)
-
-# Parameters governing the conserved quantities statistics
-Statistics:
- delta_time: 1e-3 # Time between statistics output
- tasks_per_cell: 1000
-
-# Parameters for the hydrodynamics scheme
-SPH:
- resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
- CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
-
-# Parameters related to the initial conditions
-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:
- CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table
- WithUVbackground: 1 # Enable or not the UV background
- Redshift: 0 # Redshift to use (-1 means time based redshift)
- WithMetalCooling: 1 # Enable or not the metal cooling
- ProvideVolumetricHeatingRates: 0 # User provide volumetric heating rates
- ProvideSpecificHeatingRates: 0 # User provide specific heating rates
- SelfShieldingMethod: 0 # Grackle (<= 3) or Gear self shielding method
- ConvergenceLimit: 1e-4
- MaxSteps: 20000
- Omega: 0.8
-
-Gravity:
- eta: 0.025 # Constant dimensionless multiplier for time integration.
- theta: 0.7 # Opening angle (Multipole acceptance criterion)
- epsilon: 0.1 # Softening length (in internal units).
- a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
- r_cut_max: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
- r_cut_min: 0.1 # (Optional) Cut-off in number of top-level cells below which no truncation of FMM forces are performed (this is the default value).
-
-GearChemistry:
- InitialMetallicity: 0.01295
\ No newline at end of file
diff --git a/examples/cooling_rate/generate_grackle_data.py b/examples/cooling_rate/generate_grackle_data.py
deleted file mode 100644
index 290c6611dd5a454242de02f94d8a6efa97291cca..0000000000000000000000000000000000000000
--- a/examples/cooling_rate/generate_grackle_data.py
+++ /dev/null
@@ -1,256 +0,0 @@
-#!/usr/bin/env python
-########################################################################
-#
-# Cooling rate data generation. This code is a modified version
-# of a code developed by the Grackle Team.
-#
-#
-# Copyright (c) 2013-2016, Grackle Development Team.
-#
-# Distributed under the terms of the Enzo Public Licence.
-#
-# The full license is in the file LICENSE, distributed with this
-# software.
-########################################################################
-
-"""
-Generate (or update) a hdf5 file containing
-the grackle data for each primordial chemistry
-"""
-
-import numpy as np
-from h5py import File
-from copy import deepcopy
-
-from pygrackle import \
- chemistry_data, \
- setup_fluid_container
-
-from pygrackle.utilities.physical_constants import \
- mass_hydrogen_cgs, \
- sec_per_Myr
-
-filename = "grackle.hdf5"
-
-debug = 1
-
-# metal_mass_fraction = 1e-10
-metal_mass_fraction = 0.01295
-
-tolerance = 1e-5
-max_iter = 1e8
-
-
-def generate_data(primordial_chemistry):
- """
- Compute the cooling rate of a bunch of particles and save
- them to an hdf5 file
- """
- current_redshift = 0.
-
- # Set solver parameters
- my_chemistry = chemistry_data()
- my_chemistry.use_grackle = 1
- my_chemistry.with_radiative_cooling = 1
- my_chemistry.primordial_chemistry = primordial_chemistry
- my_chemistry.metal_cooling = 1
- my_chemistry.UVbackground = 1
- my_chemistry.self_shielding_method = 0
- my_chemistry.H2_self_shielding = 0
- my_chemistry.grackle_data_file = "CloudyData_UVB=HM2012.h5"
-
- # Set units
- my_chemistry.comoving_coordinates = 0 # proper units
- my_chemistry.a_units = 1.0
- my_chemistry.a_value = 1.0 / (1.0 + current_redshift) / \
- my_chemistry.a_units
- my_chemistry.length_units = 3.085e21
- my_chemistry.density_units = 1.989e43 / my_chemistry.length_units**3
- my_chemistry.velocity_units = 20725573.785998672
- my_chemistry.time_units = my_chemistry.length_units / \
- my_chemistry.velocity_units
-
- density = mass_hydrogen_cgs
- dt = 1e-8 * sec_per_Myr / my_chemistry.time_units
- # Call convenience function for setting up a fluid container.
- # This container holds the solver parameters, units, and fields.
- temperature = np.logspace(1, 9, 1000)
- fc = setup_fluid_container(my_chemistry,
- metal_mass_fraction=metal_mass_fraction,
- temperature=temperature,
- density=density,
- converge=True,
- tolerance=tolerance,
- max_iterations=max_iter)
-
- f, data = write_input(
- fc, my_chemistry, temperature, dt, primordial_chemistry)
-
- old = deepcopy(fc)
- fc.solve_chemistry(dt)
- rate = (fc["energy"] - old["energy"]) / dt
-
- write_output(f, data, fc, rate)
-
-
-def write_output(f, data, fc, rate):
- """
- Write the output of the simulation
- """
- data.create_group("Output")
- tmp = data["Output"]
- # Rate
- dset = tmp.create_dataset("Rate", rate.shape,
- dtype=rate.dtype)
- dset[:] = rate
- # Energy
- dset = tmp.create_dataset("Energy", fc["energy"].shape,
- dtype=fc["energy"].dtype)
- dset[:] = fc["energy"]
-
- f.close()
-
-
-def simulationAreEqual(fields, data):
- """
- Check if current simulation has already
- been run
- """
- test = True
- for key in fields.keys():
- if data.attrs[key] != fields[key]:
- test = False
- break
-
- return test
-
-
-def getSameSimulation(data, my_chemistry):
- """
- Get indice of the same simulation.
- returns -1 if not found
- """
- fields = {
- "MetalCooling": my_chemistry.metal_cooling,
- "UVBackground": my_chemistry.UVbackground,
- "SelfShieldingMethod": my_chemistry.self_shielding_method,
- "MetalMassFraction": metal_mass_fraction,
- "ScaleFactor": my_chemistry.a_value
- }
- ind = 0
- for d in data:
- if "Params" in data[d]:
- tmp = data[d]["Params"]
- if simulationAreEqual(fields, tmp):
- return ind
- ind += 1
-
- return -1
-
-
-def write_input(fc, my_chemistry, temperature, dt, primordial_chemistry):
- """
- Write the simulation input in the HDF5 file
- """
- f = File(filename, "a")
- data_name = "PrimordialChemistry%i" % primordial_chemistry
- if data_name in f:
- data = f[data_name]
- else:
- print("Creating Dataset")
- data = f.create_group(data_name)
-
- i = getSameSimulation(data, my_chemistry)
- if i != -1:
- print("Updating Simulation")
- data.pop(data.keys()[i])
- data = data.create_group("Simulation %i" % len(data))
-
- # Units
- data.create_group("Units")
- tmp = data["Units"]
- tmp.attrs["Length"] = my_chemistry.length_units
- tmp.attrs["Density"] = my_chemistry.density_units
- tmp.attrs["Velocity"] = my_chemistry.velocity_units
- tmp.attrs["Time"] = my_chemistry.time_units
-
- # Parameters
- data.create_group("Params")
- tmp = data["Params"]
- tmp.attrs["MetalCooling"] = my_chemistry.metal_cooling
- tmp.attrs["UVBackground"] = my_chemistry.UVbackground
- tmp.attrs["SelfShieldingMethod"] = my_chemistry.self_shielding_method
- tmp.attrs["MetalMassFraction"] = metal_mass_fraction
- tmp.attrs["ScaleFactor"] = my_chemistry.a_value
- tmp.attrs["TimeStep"] = dt
- tmp.attrs["Tolerance"] = tolerance
- tmp.attrs["MaxIter"] = max_iter
-
- # Inputs
- data.create_group("Input")
- tmp = data["Input"]
- # energy
- dset = tmp.create_dataset("Energy", fc["energy"].shape,
- dtype=fc["energy"].dtype)
- dset[:] = fc["energy"]
- # density
- dset = tmp.create_dataset("Density", fc["density"].shape,
- dtype=fc["density"].dtype)
- dset[:] = fc["density"]
- # Temperature
- dset = tmp.create_dataset("Temperature", temperature.shape,
- dtype=temperature.dtype)
- dset[:] = temperature
-
- write_fractions(tmp, fc, primordial_chemistry)
-
- return f, data
-
-
-def write_fractions(tmp, fc, primordial_chemistry):
- """
- Write the element fractions in the HDF5 file
- """
- fields = get_fields(primordial_chemistry)
- for i in fields:
- dset = tmp.create_dataset(i, fc[i].shape,
- dtype=fc[i].dtype)
- dset[:] = fc[i] / fc["density"]
-
-
-def get_fields(primordial_chemistry):
- """
- Return the fields present in grackle
- """
- fields = [
- "metal"
- ]
- if primordial_chemistry > 0:
- fields.extend([
- "HI",
- "HII",
- "HeI",
- "HeII",
- "HeIII",
- "de"])
-
- if primordial_chemistry > 1:
- fields.extend([
- "HM",
- "H2I",
- "H2II"
- ])
-
- if primordial_chemistry > 2:
- fields.extend([
- "DI",
- "DII",
- "HDI"
- ])
- return fields
-
-
-if __name__ == "__main__":
- for i in range(4):
- print("Computing Primordial Chemistry %i" % i)
- generate_data(i)
diff --git a/examples/cooling_rate/plot_cooling.py b/examples/cooling_rate/plot_cooling.py
deleted file mode 100644
index 62bc5377c5de9357a543ac34faeb92c4de7ef063..0000000000000000000000000000000000000000
--- a/examples/cooling_rate/plot_cooling.py
+++ /dev/null
@@ -1,330 +0,0 @@
-#!/usr/bin/env python3
-"""
-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/>.
-"""
-
-from libpyswiftsim import coolingRate
-
-from copy import deepcopy
-import numpy as np
-import matplotlib.pyplot as plt
-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 = 0
-
-# reference data
-grackle_filename = "grackle.hdf5"
-compute_equilibrium = False
-
-# swift params filename
-filename = "cooling.yml"
-
-# if grackle_filename does not exist
-# use following values
-
-# density in atom / cm3
-N_rho = 1
-# with N_rho > 1, the code is not implemented to deal
-# with the reference
-if N_rho == 1:
- default_density = np.array([1.])
-else:
- default_density = np.logspace(-3, 1, N_rho)
-
-# temperature in K
-N_T = 1000
-default_temperature = np.logspace(1, 9, N_T)
-
-# time step in s
-default_dt = units.Myr * 1e-8
-default_dt = default_dt.to("s") / units.s
-
-# adiabatic index
-gamma = 5. / 3.
-
-# cooling output
-swift_output_filename = "swift.hdf5"
-
-# SCRIPT
-
-
-def get_fields(primordial_chemistry):
- fields = [
- "metal"
- ]
- if primordial_chemistry > 0:
- fields.extend([
- "HI",
- "HII",
- "HeI",
- "HeII",
- "HeIII",
- "de"])
-
- if primordial_chemistry > 1:
- fields.extend([
- "HM",
- "H2I",
- "H2II"
- ])
-
- if primordial_chemistry > 2:
- fields.extend([
- "DI",
- "DII",
- "HDI"
- ])
- return fields
-
-
-def generate_default_initial_condition(us):
- print("Generating default initial conditions")
- d = {}
- # generate grid
- rho, T = np.meshgrid(default_density, default_temperature)
- rho = deepcopy(rho.reshape(-1))
- T = T.reshape(-1)
- d["temperature"] = T
-
- # Deal with units
- m_p = constants.m_p.to("{} g".format(us["UnitLength_in_cgs"]))
- rho *= us["UnitLength_in_cgs"]**3 * m_p
- d["density"] = rho
-
- 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 / unit_time
- d["time_step"] = dt
-
- return d
-
-
-def read_grackle_data(filename, us):
- print("Reading initial conditions")
- f = File(filename, "r")
- data = f["PrimordialChemistry%i" % primordial_chemistry]
-
-<<<<<<< HEAD
- print(cooling)
- 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]
->>>>>>> 083548de07e9dd7ccb17d962bb3e253d656c0a83
-
- # read units
- tmp = data["Units"].attrs
-
- 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"][:] / float(us["UnitTemp_in_cgs"])
- d["temperature"] = T
-
- density = tmp["Density"][:] * u_den
- d["density"] = density
-
- dt = data["Params"].attrs["TimeStep"] * u_time
- d["time_step"] = dt
-
- # read fractions
- for i in get_fields(primordial_chemistry):
- d[i] = tmp[i][:]
-
- # read output
- tmp = data["Output"]
-
- energy = tmp["Energy"][:] * u_len**2 / u_time**2
- d["out_energy"] = energy
-
- rate = tmp["Rate"][:] * u_len**2 / (u_time**3)
- d["rate"] = rate
-
- f.close()
- return d
-
-
-def getParser(filename):
- with open(filename, "r") as stream:
- stream = yaml.load(stream)
- return stream
-
-
-def plot_solution(rate, data, us):
- print("Plotting solution")
- energy = data["energy"]
- rho = data["density"]
- T = data["temperature"]
-
- ref = False
- if "rate" in data:
- ref = True
- grackle_rate = data["rate"]
-
- # change units => cgs
- 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 *= float(us["UnitTemp_in_cgs"])
-
- energy *= u_l**2 / u_t**2
-
- rate *= u_l**2 / u_t**3
-
- if ref:
- grackle_rate *= u_l**2 / u_t**3
-
- # lambda cooling
- lam = rate * rho
-
- if ref:
- grackle_lam = grackle_rate * rho
-
- # do plot
- if N_rho == 1:
- # plot Lambda vs T
- plt.figure()
- plt.loglog(T, np.abs(lam),
- label="SWIFT")
- if ref:
- # plot reference
- plt.loglog(T, np.abs(grackle_lam),
- label="Grackle, Prim. Chem. %i" % primordial_chemistry)
- plt.legend()
- plt.xlabel("Temperature [K]")
- plt.ylabel("$\\Lambda$ [erg s$^{-1}$ cm$^{3}$]")
-
- # plot error vs T
- plt.figure()
- plt.plot(T, (lam - grackle_lam) / grackle_lam)
- plt.gca().set_xscale("log")
- plt.xlabel("Temperature [K]")
- plt.ylabel(
- r"$(\Lambda - \Lambda_\textrm{ref}) / \Lambda_\textrm{ref}$")
-
- plt.show()
-
- else:
- shape = [N_rho, N_T]
- cooling_time = energy / rate
- cooling_length = np.sqrt(gamma * (gamma-1.) * energy) * cooling_time
-
- cooling_length = np.log10(np.abs(cooling_length) / units.kpc.to('cm'))
-
- # reshape
- rho = rho.reshape(shape)
- T = T.reshape(shape)
- energy = energy.reshape(shape)
- cooling_length = cooling_length.reshape(shape)
-
- _min = -7
- _max = 7
- N_levels = 100
- levels = np.linspace(_min, _max, N_levels)
- plt.figure()
- plt.contourf(rho, T, cooling_length, levels)
- plt.xlabel("Density [atom/cm3]")
- plt.ylabel("Temperature [K]")
-
- ax = plt.gca()
- ax.set_xscale("log")
- ax.set_yscale("log")
-
- cbar = plt.colorbar()
- tc = np.arange(_min, _max, 1.5)
- cbar.set_ticks(tc)
- cbar.set_ticklabels(tc)
- cbar.ax.set_ylabel("Log( Cooling Length / kpc )")
-
- plt.show()
-
-
-if __name__ == "__main__":
-
- parser = getParser(filename)
- us = parser["InternalUnitSystem"]
-
- if isfile(grackle_filename) and N_rho == 1:
- d = read_grackle_data(grackle_filename, us)
- else:
- d = generate_default_initial_condition(us)
-
- # du / dt
- print("Computing cooling...")
- rate = np.zeros(d["density"].shape)
-
- if compute_equilibrium:
-<<<<<<< HEAD
- rate = wrapper.coolingRate(
- filename,
- 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"])
->>>>>>> 083548de07e9dd7ccb17d962bb3e253d656c0a83
- else:
- fields = get_fields(primordial_chemistry)
- N = len(fields)
- frac = np.zeros([len(d["density"]), N])
- for i in range(N):
- frac[:, i] = d[fields[i]]
-
- 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/examples/cooling_rate/run.sh b/examples/cooling_rate/run.sh
deleted file mode 100644
index 3564de07f4f0751a51f025fc018bf5ed5c07e64c..0000000000000000000000000000000000000000
--- a/examples/cooling_rate/run.sh
+++ /dev/null
@@ -1,17 +0,0 @@
-#!/bin/bash
-
-# Get the Grackle cooling table
-if [ ! -e CloudyData_UVB=HM2012.h5 ]
-then
- echo "Fetching the Cloudy tables required by Grackle..."
- ./getCoolingTable.sh
-fi
-
-# Generate Grackle data if not present
-if [ ! -e grackle.hdf5 ]
-then
- echo "Generating Grackle Data..."
- ./generate_grackle_data.py
-fi
-
-./plot_cooling.py
diff --git a/pyswiftsim/dev.py b/pyswiftsim/dev.py
deleted file mode 100644
index 64561b3bcb36adc66df701352e5705d95afead36..0000000000000000000000000000000000000000
--- a/pyswiftsim/dev.py
+++ /dev/null
@@ -1,88 +0,0 @@
-#!/usr/bin/env python3
-"""
-This file is part of PYSWIFT.
-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/>.
-"""
-
-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
- from it. (for dev, no guarantee of working).
- """
-
- struct_format = ""
- struct_attribute = []
-
- with open(filename, "r") as f:
- line = ""
- # skip until reaching struct
- while ("struct" not in line):
- line = f.readline()
-
- count = 0
- in_struct = False
-
- while (not in_struct or count != 0):
- # count { and }
- count += line.count("{")
- if (count > 0):
- in_struct = True
- count -= line.count("}")
- line = f.readline()
-
- for k in symbol_table.keys():
- if k not in line:
- continue
-
- i = line.index(k) + len(k)
-
- if (line[i] == "*"):
- struct_format += "p"
- i += 1
- else:
- struct_format += symbol_table[k]
-
- j = line.index(";")
- struct_attribute.append(line[i:j])
-
- return struct_format, struct_attribute
-
-
-if __name__ == "__main__":
-
- 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}
- """.format(attr=attr,
- form=struct_format)
- print(txt)
diff --git a/pyswiftsim/structure.py b/pyswiftsim/structure.py
deleted file mode 100644
index 184dff7107432e44cf71836d197c732db18a9ef0..0000000000000000000000000000000000000000
--- a/pyswiftsim/structure.py
+++ /dev/null
@@ -1,563 +0,0 @@
-"""
-This file is part of PYSWIFT.
-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/>.
-"""
-from pyswiftsim import wrapper
-
-import struct
-from ctypes import *
-
-PARSER_MAX_LINE_SIZE = 256
-PARSER_MAX_NO_OF_PARAMS = 256
-PARSER_MAX_NO_OF_SECTIONS = 64
-
-
-######################################################################
-# #
-# SwiftStruct #
-# #
-######################################################################
-class SwiftStruct(struct.Struct):
- """
- Abstract class constructing a wrapper around a C structure
-
- Parameters
- ----------
-
- struct_format: str
- Text defining the data type in the structure (see struct module for
- more informations)
-
- data: str, list
- Bytes string containing all the data of a structure or
- dictionary of attribute (all the attribute should be present)
-
- parent: :class:`SwiftStruct`
- Parent containing the data (usefull for :class:`ArrayStruct`)
- """
- def __init__(self, struct_format, data, parent):
- super().__init__(struct_format)
- self.parent = parent # parent for ArrayStruct
-
- if isinstance(data, bytes):
- self.data = data # bytes string data
- elif isinstance(data, dict):
- tmp = []
- for i in self.struct_name:
- tmp.append(data[i])
-
- self.data = self.pack(*tmp)
- else:
- raise ValueError("Data should be either bytes or dict, "
- "received ", type(data))
-
- @property
- def struct_name(self):
- """
- List of the structure attribute names.
-
- Should be implemented for each sub class
- """
- return self.__class__._name
-
- @property
- def struct_format(self):
- """
- String containing the data type of each attribute
- """
- return self.__class__._format
-
- @property
- def struct_substruct(self):
- """
- Dictionary containing the class of each substructure.
- See for example SwiftParams
- """
- return {}
-
- def _getInfoFromName(self, name):
- """
- Compute index, format and size from an attribute name.
-
- Parameters
- ----------
-
- name: str
- Name of an attribute
-
- Returns
- -------
-
- i: int, slice
- index of the attribute
-
- form: char
- data type (see module struct for more information)
-
- n: int
- number of element in attribute
- """
- struct_size, form = self.struct_size_format
- i = self.struct_name.index(name)
- form = form[i]
- n = struct_size[i]
- i = sum(struct_size[:i])
-
- if n > 1:
- i = slice(i, i+n)
-
- return i, form, n
-
- @property
- def struct_size_format(self):
- """
- Compute size and format for each field
-
- Returns
- -------
-
- out_nber: list
- number of element for each attribute
-
- out_form: list
- format for each attribute
- """
- out_nber = []
- out_form = []
- form = self.struct_format
-
- N = len(form)
- ii = 0
- while ii < N:
- v = form[ii]
- count = ""
- while v.isdigit():
- count += v
- ii += 1
- v = form[ii]
-
- if count == "":
- count = 1
-
- if v == "s":
- count = 1
-
- count = int(count)
- out_nber.append(count)
- out_form.append(v)
- ii += 1
-
- return out_nber, out_form
-
- def __str__(self, tab=""):
- txt = tab + "%s:\n" % type(self)
- for name in self.struct_name:
- d = getattr(self, name)
- txt += tab + "\t%s: " % name
- if isinstance(d, SwiftStruct):
- txt += d.__str__(tab+"\t")
- else:
- txt += "%s\n" % d
- return txt
-
- def __getattr__(self, name):
- # case where the attribute is not in the structure
- if name not in self.struct_name:
- return object.__getattribute__(self, name)
-
- # case where the attribute is in the structure
- else:
- i, form, n = self._getInfoFromName(name)
-
- data = self.unpack(self.data)
- if n == 1:
- # if substruct
- if name in self.struct_substruct:
- d = self.struct_substruct[name]
- cl = d["class"]
- tmp = []
- size = struct.calcsize(cl._format)
- for j in range(d["size"]):
- data_tmp = data[i][size*j:(j+1)*size]
- tmp.append(cl(data_tmp, parent=self))
-
- if d["size"] == 1:
- return tmp[0]
- else:
- return tmp
- # other case => array
- else:
- return data[i]
-
- else:
- # transform scalar -> vector
- nform = str(n) + form
- i = slice(i.start, i.start+n)
- data = data[i]
-
- # compress data and create return struct
- data = struct.pack(nform, *data)
- return ArrayStruct(nform, data, self, name)
-
- def __setattr__(self, name, value):
- # case where the attribute is not in the structure
- if name not in self.struct_name:
- object.__setattr__(self, name, value)
- # case where the attribute is in the structure
- else:
- data = list(self.unpack(self.data))
-
- i, form, n = self._getInfoFromName(name)
-
- if isinstance(value, ArrayStruct):
- value = value.getArray()
-
- data[i] = value
- self.data = self.pack(*data)
-
-
-######################################################################
-# #
-# ArrayStruct #
-# #
-######################################################################
-class ArrayStruct(SwiftStruct):
- _name = [
- "array_data"
- ]
-
- def __init__(self, struct_format, data, parent, name):
- super().__init__(struct_format, data, parent)
- self._format = struct_format
- self._name = name
-
- def __getitem__(self, ii):
- data = self.unpack(self.data)
- data = self._clean(data)
- return data[ii]
-
- def __setitem__(self, ii, value):
- data = list(self.unpack(self.data))
- data[ii] = value
- setattr(self.parent, self._name, data)
-
- def _clean(self, data):
- if self._format[-1] == "c":
- data = b"".join(data)
- # find end of txt character
- n = data.index(b"\x00")
- data = data[:n].decode("utf-8")
-
- return data
-
- def __str__(self, tab=""):
- data = self.unpack(self.data)
- data = self._clean(data)
- return tab + str(data) + "\n"
-
- def getArray(self):
- return self.unpack(self.data)
-
- @property
- def struct_format(self):
- return self._format
-
-
-######################################################################
-# #
-# UnitSystem #
-# #
-######################################################################
-class UnitSystem(SwiftStruct):
- _format = "ddddd"
- _name = [
- "UnitMass_in_cgs",
- "UnitLength_in_cgs",
- "UnitTime_in_cgs",
- "UnitCurrent_in_cgs",
- "UnitTemperature_in_cgs",
- ]
-
- def __init__(self, data, parent=None):
- super().__init__(self.struct_format, data, parent)
-
-
-######################################################################
-# #
-# chemistry_part_data #
-# #
-######################################################################
-class ChemistryPartData(SwiftStruct):
- _format = "f"
- _name = [
- "initial_metallicity"
- ]
-
- def __init__(self, data, parent=None):
- super().__init__(self.struct_format, data, parent)
-
-
-######################################################################
-# #
-# Part #
-# #
-######################################################################
-class Part(SwiftStruct):
- _format = "qP3d3f3ffffffN7f4{chemistry}sc".format(
- chemistry=struct.calcsize(ChemistryPartData._format)
- )
- _name = [
- "id",
- "gpart",
- "pos",
- "vel",
- "a_hydro",
- "h",
- "mass",
- "rho",
- "entropy",
- "entropy_dt",
- "last_offset",
- "density",
- "chemistry_data",
- "time_bin"
- ]
-
- def __init__(self, data, parent=None):
- super().__init__(self.struct_format, data, parent)
- print("ERROR, need to fix density/time_bin")
-
-
-######################################################################
-# #
-# Parameter #
-# #
-######################################################################
-class Parameter(SwiftStruct):
- _format = "{line_size}c{line_size}cii".format(
- line_size=PARSER_MAX_LINE_SIZE
- )
-
- _name = [
- "name",
- "value",
- "used",
- "is_default"
- ]
-
- def __init__(self, data, parent=None):
- super().__init__(self.struct_format, data, parent)
-
-
-######################################################################
-# #
-# Section #
-# #
-######################################################################
-class Section(SwiftStruct):
- _format = "{line_size}c".format(
- line_size=PARSER_MAX_LINE_SIZE
- )
- _name = [
- "name"
- ]
-
- def __init__(self, data, parent=None):
- super().__init__(self.struct_format, data, parent)
-
-
-######################################################################
-# #
-# SwiftParams #
-# #
-######################################################################
-class SwiftParams(SwiftStruct):
- _format = "{sec}s{data}sii{line_size}c".format(
- sec=struct.calcsize(Section._format)*PARSER_MAX_NO_OF_SECTIONS,
- data=struct.calcsize(Parameter._format)*PARSER_MAX_NO_OF_PARAMS,
- line_size=PARSER_MAX_LINE_SIZE
- )
-
- _name = [
- "section",
- "data_params",
- "sectionCount",
- "paramCount",
- "filename"
- ]
-
- def __init__(self, data, parent=None):
- super().__init__(self.struct_format, data, parent)
-
- @property
- def struct_substruct(self):
- sec = {
- "class": Section,
- "size": PARSER_MAX_NO_OF_SECTIONS
- }
-
- param = {
- "class": Parameter,
- "size": PARSER_MAX_NO_OF_PARAMS
- }
- return {
- "section": sec,
- "data_params": param
- }
-
-
-######################################################################
-# #
-# PhysConst #
-# #
-######################################################################
-class PhysConst(SwiftStruct):
- _format = "ddddddddddddddddddddd"
- _name = [
- "const_newton_G",
- "const_speed_light_c",
- "const_planck_h",
- "const_planck_hbar",
- "const_boltzmann_k",
- "const_avogadro_number",
- "const_thomson_cross_section",
- "const_stefan_boltzmann",
- "const_electron_charge",
- "const_electron_volt",
- "const_electron_mass",
- "const_proton_mass",
- "const_year",
- "const_astronomical_unit",
- "const_parsec",
- "const_light_year",
- "const_solar_mass",
- "const_earth_mass",
- "const_T_CMB_0",
- "const_primordial_He_fraction",
- "const_reduced_hubble"
- ]
-
- def __init__(self, data, parent=None):
- super().__init__(self.struct_format, data, parent)
-
-
-class CoolingFunctionData(SwiftStruct):
- cooling_type = wrapper.configGetCooling()
- if cooling_type == "const_lambda":
- _format = "dffff"
- _name = [
- "lambda",
- "hydrogen_mass_abundance",
- "mean_molecular_weight",
- "min_energy",
- "cooling_tstep_mult"
- ]
- elif "grackle" in cooling_type:
- _format = "200cid{code_units}c{chemistry}ciiiifif".format(
- code_units=56,
- chemistry=248
- )
- _name = [
- "cloudy_table",
- "with_uv_background",
- "redshift",
- "code_units",
- "chemistry",
- "with_metal_cooling",
- "provide_volumetric_heating_rate",
- "provide_specificy_heating_rate",
- "self_shielding_method",
- "convergence_limit",
- "max_step",
- "omega"
- ]
-
- else:
- raise ValueError(
- "Cooling Type %s not implemented" % cooling_type)
-
- def __init__(self, data, parent=None):
- super().__init__(self.struct_format, data, parent)
-
-
-class ChemistryGlobalData(SwiftStruct):
- _format = "f"
- _name = [
- "initial_metallicity",
- ]
-
- 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_Balsara_eps',
- 'a_factor_grav_accel',
- 'a_factor_hydro_accel',
- 'z',
- 'H',
- 'critical_density',
- 'critical_density_0',
- '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',
- 'hydro_kick_corr_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)
diff --git a/src/chemistry_wrapper.c b/src/chemistry_wrapper.c
deleted file mode 100644
index 9c6cdc09f6f3919de4d69aa5ab1bab516808ed79..0000000000000000000000000000000000000000
--- a/src/chemistry_wrapper.c
+++ /dev/null
@@ -1,63 +0,0 @@
-/*******************************************************************************
- * 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/>.
- *
- ******************************************************************************/
-
-#include "pyswiftsim_tools.h"
-#include "chemistry_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 ChemistryGlobalData
- */
-PyObject* pychemistry_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 chemistry_global_data chemistry;
-
- /* init cooling */
- chemistry_init_backend(params, us, pconst, &chemistry);
-
- /* construct python object */
- PyObject *pychemistry = pytools_return(&chemistry, class_chemistry_global_data);
-
- return pychemistry;
-}
diff --git a/src/chemistry_wrapper.h b/src/chemistry_wrapper.h
deleted file mode 100644
index 31ab1593f820845104c044fb64a312c767bcdd6b..0000000000000000000000000000000000000000
--- a/src/chemistry_wrapper.h
+++ /dev/null
@@ -1,27 +0,0 @@
-/*******************************************************************************
- * 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_CHEMISTRY_H__
-#define __PYSWIFTSIM_CHEMISTRY_H__
-
-#include "pyswiftsim_tools.h"
-
-PyObject* pychemistry_init(PyObject* self, PyObject* args);
-
-#endif // __PYSWIFTSIM_CHEMISTRY_H__
-
diff --git a/src/cooling_wrapper.c b/src/cooling_wrapper.c
index 8def5d2d1e23507a0edd4f438816df7a02aff522..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 100644
--- a/src/cooling_wrapper.c
+++ b/src/cooling_wrapper.c
@@ -1,179 +0,0 @@
-/*******************************************************************************
- * 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/>.
- *
- ******************************************************************************/
-#include "cooling_wrapper.h"
-
-/* Local includes */
-#include "pyswiftsim_tools.h"
-#include "part_wrapper.h"
-
-
-/**
- * @brief Set the cooling element fractions
- *
- * @param xp The #xpart to set
- * @param frac The numpy array containing the fractions (id, element)
- * @param idx The id (in frac) of the particle to set
- */
-void pycooling_set_fractions(struct xpart *xp, PyArrayObject* frac, const int idx) {
-#ifdef COOLING_GRACKLE
- struct cooling_xpart_data *data = &xp->cooling_data;
- data->metal_frac = *(float*)PyArray_GETPTR2(frac, idx, 0);
-
-#ifdef COOLING_GRACKLE
-#if COOLING_GRACKLE_MODE > 0
- data->HI_frac = *(float*)PyArray_GETPTR2(frac, idx, 1);
- data->HII_frac = *(float*)PyArray_GETPTR2(frac, idx, 2);
- data->HeI_frac = *(float*)PyArray_GETPTR2(frac, idx, 3);
- data->HeII_frac = *(float*)PyArray_GETPTR2(frac, idx, 4);
- data->HeIII_frac = *(float*)PyArray_GETPTR2(frac, idx, 5);
- data->e_frac = *(float*)PyArray_GETPTR2(frac, idx, 6);
-#endif // COOLING_GRACKLE_MODE
-#if COOLING_GRACKLE_MODE > 1
- data->HM_frac = *(float*)PyArray_GETPTR2(frac, idx, 7);
- data->H2I_frac = *(float*)PyArray_GETPTR2(frac, idx, 8);
- data->H2II_frac = *(float*)PyArray_GETPTR2(frac, idx, 9);
-#endif // COOLING_GRACKLE_MODE
-#if COOLING_GRACKLE_MODE > 2
- data->DI_frac = *(float*)PyArray_GETPTR2(frac, idx, 10);
- data->DII_frac = *(float*)PyArray_GETPTR2(frac, idx, 11);
- data->HDI_frac = *(float*)PyArray_GETPTR2(frac, idx, 12);
-#endif // COOLING_GRACKLE_MODE
-#endif // COOLING_GRACKLE
-
-#else
- message("Not implemented");
-#endif
-}
-
-/**
- * @brief Compute cooling rate
- *
- * args is expecting pyswiftsim classes in the following order:
- * PhysConst, UnitSystem and CoolingFunctionData.
- * Then two numpy arrays (density and specific energy) and an optional
- * float for the time step
- *
- * @param self calling object
- * @param args arguments
- * @return cooling rate
- */
-PyObject* pycooling_rate(PyObject* self, PyObject* args) {
- import_array();
-
-
- char *filename;
-
- PyArrayObject *rho;
- PyArrayObject *energy;
- PyArrayObject *fractions = NULL;
-
- float dt = 1e-3;
- int with_cosmo = 0;
-
- /* parse argument */
- if (!PyArg_ParseTuple(
- args, "siOO|fO", &filename, &with_cosmo,
- &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);
-
-
- /* Initialize everything */
- struct swift_params params;
- parser_read_file(filename, ¶ms);
-
- struct unit_system us;
- units_init_from_params(&us, ¶ms, "InternalUnitSystem");
-
- struct phys_const pconst;
- phys_const_init(&us, ¶ms, &pconst);
-
- struct cooling_function_data cooling;
-
- /* 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);
-
- /* 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);
-
- if (fractions != NULL)
- pycooling_set_fractions(&xp, fractions, i);
- else
- cooling_first_init_part(&pconst, &us, &cosmo, &cooling, &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 /= dt;
-#else
- message("Not implemented");
- //*tmp = cooling_rate(pconst, us, cosmo, cooling, &p, &xp);
-#endif
- }
-
- return (PyObject *) rate;
-
-}
diff --git a/src/cooling_wrapper.h b/src/cooling_wrapper.h
index a9eeeb4d6a5061637795536f76af71cc439c6747..86841cff6fcfed69a26a4fa89e42bd0ab9d1ab0e 100644
--- a/src/cooling_wrapper.h
+++ b/src/cooling_wrapper.h
@@ -19,7 +19,7 @@
#ifndef __PYSWIFTSIM_COOLING_H__
#define __PYSWIFTSIM_COOLING_H__
-#include "pyswiftsim_tools.h"
+#include "pyswiftsim_includes.h"
PyObject* pycooling_rate(PyObject* self, PyObject* args);
diff --git a/src/cosmology_wrapper.c b/src/cosmology_wrapper.c
deleted file mode 100644
index d9789187106f40ef505d72e3d1a2456ebbb59c63..0000000000000000000000000000000000000000
--- a/src/cosmology_wrapper.c
+++ /dev/null
@@ -1,67 +0,0 @@
-/*******************************************************************************
- * 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;
-}
diff --git a/src/cosmology_wrapper.h b/src/cosmology_wrapper.h
deleted file mode 100644
index 8a674060e1c79f01090f441ac2ad7e8cdc1ff340..0000000000000000000000000000000000000000
--- a/src/cosmology_wrapper.h
+++ /dev/null
@@ -1,27 +0,0 @@
-/*******************************************************************************
- * 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__
-
diff --git a/src/parser_wrapper.c b/src/parser_wrapper.c
deleted file mode 100644
index 9f4920fa55529604c19f02c8c0d591af9cb50bac..0000000000000000000000000000000000000000
--- a/src/parser_wrapper.c
+++ /dev/null
@@ -1,41 +0,0 @@
-/*******************************************************************************
- * 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/>.
- *
- ******************************************************************************/
-#include "parser_wrapper.h"
-#include "pyswiftsim_tools.h"
-
-
-PyObject* pyparser_read_file(PyObject *self, PyObject *args)
-{
-
- char *filename;
-
- /* parse argument */
- if (!PyArg_ParseTuple(args, "s", &filename))
- return NULL;
-
- struct swift_params params;
-
- /* parse file */
- parser_read_file(filename, ¶ms);
-
- /* create return python object */
- PyObject* obj = pytools_return(¶ms, class_swift_params);
-
- return obj;
-}
diff --git a/src/parser_wrapper.h b/src/parser_wrapper.h
deleted file mode 100644
index 223e6663d1e042857d2552bc7e5fe27ee2014b0f..0000000000000000000000000000000000000000
--- a/src/parser_wrapper.h
+++ /dev/null
@@ -1,36 +0,0 @@
-/*******************************************************************************
- * 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_PARSER_H__
-#define __PYSWIFTSIM_PARSER_H__
-
-#include <Python.h>
-
-/**
- * @brief Read swift parameters
- *
- * args is expecting a string.
- *
- * @param self calling object
- * @param args arguments
- * @return SwiftParams
- */
-PyObject* pyparser_read_file(PyObject *self, PyObject *args);
-
-
-#endif // __PYSWIFTSIM_PARSER_H__
diff --git a/src/part_wrapper.c b/src/part_wrapper.c
deleted file mode 100644
index d8cefd848401a1e5eedb5839e5f12188e2c7738c..0000000000000000000000000000000000000000
--- a/src/part_wrapper.c
+++ /dev/null
@@ -1,30 +0,0 @@
-/*******************************************************************************
- * 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/>.
- *
- ******************************************************************************/
-#include "pyswiftsim_tools.h"
-
-#include <Python.h>
-#include <stdlib.h>
-#include <string.h>
-
-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
deleted file mode 100644
index bdc3d52b5807871a4c0e656150686032917f17d8..0000000000000000000000000000000000000000
--- a/src/part_wrapper.h
+++ /dev/null
@@ -1,27 +0,0 @@
-/*******************************************************************************
- * 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_PART_H__
-#define __PYSWIFTSIM_PART_H__
-
-#include <Python.h>
-#include "pyswiftsim_tools.h"
-
-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_includes.h
similarity index 100%
rename from src/pyswiftsim.h
rename to src/pyswiftsim_includes.h
diff --git a/src/pyswiftsim_tools.c b/src/pyswiftsim_tools.c
index 5b872ec5a283030da0420dbb0bf05a38b43e774d..4fd232542387b913e7604f99cb95eb18ae3f8e70 100644
--- a/src/pyswiftsim_tools.c
+++ b/src/pyswiftsim_tools.c
@@ -19,196 +19,6 @@
#include "pyswiftsim_tools.h"
-const size_t class_size[class_count] = {
- sizeof(struct unit_system),
- sizeof(struct part),
- sizeof(struct swift_params),
- sizeof(struct phys_const),
- sizeof(struct cooling_function_data),
- sizeof(struct chemistry_global_data),
- sizeof(struct cosmology)
-};
-
-char *class_name[class_count] = {
- "UnitSystem",
- "Part",
- "SwiftParams",
- "PhysConst",
- "CoolingFunctionData",
- "ChemistryGlobalData",
- "Cosmology"
-};
-
-
-/**
- * @brief Import a python library/function
- *
- * @param module Module name
- * @param object_name object to import from module
- * @return (PyObject) imported object
- */
-PyObject* pytools_import(char* module_name, char* object_name)
-{
- /* load module */
- PyObject *module;
-
- module = PyImport_ImportModule(module_name);
- if (module == NULL)
- pyerror("Failed to import module '%s'.", module_name);
-
- /* get module dictionary */
- PyObject *dict;
-
- dict = PyModule_GetDict(module);
- Py_DECREF(module);
-
- if (dict == NULL)
- pyerror("Failed to get module '%s' dictionary", module_name);
-
- /* get right object */
- PyObject *python_obj = PyDict_GetItemString(dict, object_name);
- if (python_obj == NULL)
- pyerror("Object %s does not exist in module %s", object_name, module_name);
-
- Py_INCREF(python_obj);
- return python_obj;
-}
-
-
-/**
- * @brief Construct a python object from a C struct
- *
- * @param p pointer to the C struct
- * @param class #swift_class of the pointer
- * @return PyObject of the required class
- */
-PyObject* pytools_return(void *p, int class)
-{
-
- PyObject *python_class;
- size_t nber_bytes;
-
- char module_name[STRING_SIZE] = "pyswiftsim.structure";
-
- /* check class */
- if (class >= class_count)
- pyerror("Class %i does not exists", class);
-
- /* get class information */
- nber_bytes = class_size[class];
- char *class_pyname = class_name[class];
-
- /* import python class */
- python_class = pytools_import(module_name, class_pyname);
- if (!python_class)
- return NULL;
-
- if (!PyCallable_Check(python_class))
- {
- Py_DECREF(python_class);
- message("Redo");
- pyerror("Unable to import class %30s from %.30s", class_pyname, module_name);
- }
-
- /* create object */
- PyObject *object, *args;
-
- args = PyTuple_Pack(1, PyBytes_FromStringAndSize((char *) p, nber_bytes));
-
- object = PyObject_CallObject(python_class, args);
-
- Py_DECREF(args);
- Py_DECREF(python_class);
-
- return object;
-
-}
-
-/**
- * @brief get type name in C
- *
- * @param obj PyObject from which to get the name0
- * @return type name
- */
-char* pytools_get_type_name(PyObject *obj)
-{
- /* check input */
- if (!obj)
- pyerror("Requires a non null object");
-
- /* get object type */
- PyObject *type = (PyObject *)obj->ob_type;
- if (type == NULL)
- pyerror("Unable to get type");
-
- /* get object name */
- PyObject* recv = PyObject_Str(type);
- if (!recv)
- {
- Py_DECREF(recv);
- pyerror("Unable to get string representation");
- }
-
- /* transform to C */
- Py_ssize_t size;
- char *name = PyUnicode_AsUTF8AndSize(recv, &size);
- Py_DECREF(recv);
-
- if (!name)
- pyerror("Unable to convert string to char");
-
- return name;
-}
-
-
-/**
- * @brief Construct a C struct from a python object
- *
- * @param obj pointer to the python object
- * @param class #swift_class of the pointer
- * @return pointer to the required struct
- */
-char* pytools_construct(PyObject* obj, int class)
-{
- char *module_name = "pyswiftsim.structure";
- char *class_pyname;
-
- /* check python class */
- if (class >= class_count)
- pyerror("Class %i does not exists", class);
-
- /* get class information */
- class_pyname = class_name[class];
-
- /* import class */
- PyObject *pyclass = pytools_import(module_name, class_pyname);
- if (!pyclass)
- return NULL;
-
- /* check if classes correspond */
- int test = PyObject_IsInstance(obj, pyclass);
- Py_DECREF(pyclass);
- if (!test)
- {
- char *recv = pytools_get_type_name(obj);
- if (!recv)
- return NULL;
- pyerror("Expecting class %s, received %s", class_pyname, recv);
- }
-
- /* copy python class' data to C */
- PyObject* data = PyObject_GetAttrString(obj, "data");
-
- if (!data)
- pyerror("Unable to get the attribute 'data'");
-
- char *ret = PyBytes_AsString(data);
-
- Py_DECREF(data);
- return ret;
-}
-
-
/**
* @brief Check if object is the expected PyArray type
*
diff --git a/src/pyswiftsim_tools.h b/src/pyswiftsim_tools.h
index 3cc461c5019251511dfc91d04b73da64d3b49d5b..b7dd07b2b5972cb8b77d1afe0ffdb0b1dc7855fc 100644
--- a/src/pyswiftsim_tools.h
+++ b/src/pyswiftsim_tools.h
@@ -24,25 +24,6 @@
#define DIM 3
#define STRING_SIZE 200
-/* Set the error message for python */
-#define pyerror(s, ...) \
- ({ \
- char error_msg[STRING_SIZE]; \
- sprintf(error_msg, "\n%s:%s():%i: " s, __FILE__, \
- __FUNCTION__, __LINE__, ##__VA_ARGS__); \
- PyErr_SetString(PyExc_RuntimeError, error_msg); \
- return FAIL; \
- })
-
-/* debugging function */
-#define pydebug(s, ...) \
- ({ \
- if (!PyErr_Occurred()) \
- PyErr_Print(); \
- fprintf(stderr, "\n%s:%s():%i: Test " s "\n", __FILE__, \
- __FUNCTION__, __LINE__, ##__VA_ARGS__); \
- })
-
#define IMPORT_ARRAY1(ret) \
{ \
if (PyArray_API == NULL) \
@@ -54,38 +35,6 @@
#define IMPORT_ARRAY() IMPORT_ARRAY1(NUMPY_IMPORT_ARRAY_RETVAL);
-/* Enum swift classes */
-enum swift_class {
- class_unit_system,
- class_part,
- class_swift_params,
- class_phys_const,
- class_cooling_function_data,
- class_chemistry_global_data,
- class_cosmology,
- class_count /* should always be last! */
-};
-
-/* size of each structure in enum swift_class */
-extern const size_t class_size[];
-/* name of each Python class representing a swift class */
-extern char *class_name[];
-
-/* error code in pyswiftsim */
-enum error_code {
- FAIL = 0, // ensure NULL == FAIL
- SUCCESS,
-};
-
-
-PyObject* pytools_return(void* p, int class);
-
-char* pytools_construct(PyObject* obj, int class);
-
-PyObject* pytools_import(char* module, char* object_name);
-
-char* pytools_get_type_name(PyObject *obj);
-
int pytools_check_array(PyArrayObject *obj, int dim, int type, const char* name);
#endif // __PYSWIFTSIM_TOOLS_H__
diff --git a/src/units_wrapper.c b/src/units_wrapper.c
deleted file mode 100644
index 5fc066f1e28e155e9cf268a50b564ab4ae8d9a42..0000000000000000000000000000000000000000
--- a/src/units_wrapper.c
+++ /dev/null
@@ -1,78 +0,0 @@
-/*******************************************************************************
- * 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/>.
- *
- ******************************************************************************/
-#include "pyswiftsim_tools.h"
-
-#include <Python.h>
-#include <stdlib.h>
-#include <string.h>
-
-PyObject* pyunit_system_test_struct(PyObject *self, PyObject *args)
-{
-
- size_t N = sizeof(struct unit_system);
-
- /* initialize array */
- struct unit_system *us = malloc(N);
- us->UnitMass_in_cgs = 1.;
- us->UnitLength_in_cgs = 2.;
- us->UnitTime_in_cgs = 3.;
- us->UnitCurrent_in_cgs = 4.;
- us->UnitTemperature_in_cgs = 5.;
-
- /* construct python object */
- PyObject *object = pytools_return(us, class_unit_system);
- free(us);
-
- return object;
-}
-
-PyObject* pyunit_system_init(PyObject *self, PyObject *args)
-{
- PyObject* parser;
- char* section = "InternalUnitSystem";
-
- /* parse arguement */
- if (!PyArg_ParseTuple(args, "O|s", &parser, §ion))
- return NULL;
-
- struct swift_params *params = (struct swift_params*) pytools_construct(parser, class_swift_params);
-
- if (params == NULL)
- return NULL;
-
- /* initialize units */
- struct unit_system us;
- units_init_from_params(&us, params, section);
-
- /* initialize phys_const */
- struct phys_const pconst;
-
- phys_const_init(&us, params, &pconst);
-
- /* create python object to return */
- PyObject *pyus = pytools_return(&us, class_unit_system);
- if (pyus == NULL)
- return NULL;
-
- PyObject *pypconst = pytools_return(&pconst, class_phys_const);
- if (pypconst == NULL)
- return NULL;
-
- return PyTuple_Pack(2, pyus, pypconst);
-}
diff --git a/src/units_wrapper.h b/src/units_wrapper.h
deleted file mode 100644
index fa42a09dda6961ba1065dbd65bc785c3cacc04ea..0000000000000000000000000000000000000000
--- a/src/units_wrapper.h
+++ /dev/null
@@ -1,46 +0,0 @@
-/*******************************************************************************
- * 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_UNITS_H__
-#define __PYSWIFTSIM_UNITS_H__
-
-#include <Python.h>
-
-/**
- * @brief Test function for the struct
- *
- * args is expecting no argument.
- *
- * @param self calling object
- * @param args arguments
- * @return UnitSystem
- */
-PyObject* pyunit_system_test_struct(PyObject *self, PyObject *args);
-
-/**
- * @brief Initialize the unit system
- *
- * args is expecting a pyswift SwiftParms object.
- *
- * @param self calling object
- * @param args arguments
- * @return UnitSystem
- */
-PyObject* pyunit_system_init(PyObject *self, PyObject *args);
-
-#endif // __PYSWIFTSIM_UNITS_H__
diff --git a/src/wrapper.c b/src/wrapper.c
index 65dbad5592af49b1ca4900b3c99ad357edd56031..af15b99499ce9bd9ca7ea085da4ad418cf9bc856 100644
--- a/src/wrapper.c
+++ b/src/wrapper.c
@@ -16,13 +16,7 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
-#include "units_wrapper.h"
-#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 <Python.h>
#include <math.h>
@@ -33,18 +27,6 @@
static PyMethodDef wrapper_methods[] = {
- {"parserReadFile", pyparser_read_file, METH_VARARGS,
- "Read a swift params file."},
-
- {"unitSystemInit", pyunit_system_init, METH_VARARGS,
- "Construct a unit_system object and return it."},
-
- {"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"
diff --git a/test/test_struct.py b/test/test_struct.py
deleted file mode 100644
index 2d9df52e88b5a247619ce1df6065bf7c72b6076f..0000000000000000000000000000000000000000
--- a/test/test_struct.py
+++ /dev/null
@@ -1,29 +0,0 @@
-#!/usr/bin/env python3
-"""
-This file is part of PYSWIFT.
-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/>.
-"""
-
-from pyswiftsim import wrapper
-
-filename = "parameter_example.yml"
-params = wrapper.parserReadFile(filename)
-
-print(params)
-
-us = wrapper.unitSystemInit(params)
-
-print(us)