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

Remove everything to start again

parent a05c094e
No related branches found
No related tags found
1 merge request!6Create CI
Show whitespace changes
Inline Side-by-side
Showing
with 0 additions and 2468 deletions
examples/cooling_box/CloudyData_UVB=HM2012.h5 deleted 100644 → 0
+ 0
− 0
View file @ a05c094e
File deleted
examples/cooling_box/cooling.yml deleted 100644 → 0
+ 0
− 60
View file @ a05c094e
# 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
examples/cooling_box/generate_grackle_data.py deleted 100644 → 0
+ 0
− 262
View file @ a05c094e
#!/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)
examples/cooling_box/getCoolingTable.sh deleted 100644 → 0
+ 0
− 2
View file @ a05c094e
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5
examples/cooling_box/plot_cooling.py deleted 100644 → 0
+ 0
− 335
View file @ a05c094e
#!/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)
examples/cooling_box/run.sh deleted 100644 → 0
+ 0
− 17
View file @ a05c094e
#!/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
examples/cooling_rate/cooling.yml deleted 100644 → 0
+ 0
− 71
View file @ a05c094e
# 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
examples/cooling_rate/generate_grackle_data.py deleted 100644 → 0
+ 0
− 256
View file @ a05c094e
#!/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)
examples/cooling_rate/plot_cooling.py deleted 100644 → 0
+ 0
− 330
View file @ a05c094e
#!/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)
examples/cooling_rate/run.sh deleted 100644 → 0
+ 0
− 17
View file @ a05c094e
#!/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
pyswiftsim/dev.py deleted 100644 → 0
+ 0
− 88
View file @ a05c094e
#!/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)
pyswiftsim/structure.py deleted 100644 → 0
+ 0
− 563
View file @ a05c094e
"""
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)
src/chemistry_wrapper.c deleted 100644 → 0
+ 0
− 63
View file @ a05c094e
/*******************************************************************************
* 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;
}
src/chemistry_wrapper.h deleted 100644 → 0
+ 0
− 27
View file @ a05c094e
/*******************************************************************************
* 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__
src/cooling_wrapper.c deleted 100644 → 0
+ 0
− 179
View file @ a05c094e
/*******************************************************************************
* 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, &params);
struct unit_system us;
units_init_from_params(&us, &params, "InternalUnitSystem");
struct phys_const pconst;
phys_const_init(&us, &params, &pconst);
struct cooling_function_data cooling;
/* init cooling */
cooling_init_backend(&params, &us, &pconst, &cooling);
struct chemistry_global_data chemistry;
chemistry_init_backend(&params, &us, &pconst, &chemistry);
struct cosmology cosmo;
if (with_cosmo)
cosmology_init(&params, &us, &pconst, &cosmo);
else
cosmology_init_no_cosmo(&cosmo);
struct part p;
struct 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;
}
src/cooling_wrapper.h deleted 100644 → 0
+ 0
− 27
View file @ a05c094e
/*******************************************************************************
* 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_COOLING_H__
#define __PYSWIFTSIM_COOLING_H__
#include "pyswiftsim_tools.h"
PyObject* pycooling_rate(PyObject* self, PyObject* args);
#endif // __PYSWIFTSIM_COOLING_H__
src/cosmology_wrapper.c deleted 100644 → 0
+ 0
− 67
View file @ a05c094e
/*******************************************************************************
* 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;
}
src/cosmology_wrapper.h deleted 100644 → 0
+ 0
− 27
View file @ a05c094e
/*******************************************************************************
* 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__
src/parser_wrapper.c deleted 100644 → 0
+ 0
− 41
View file @ a05c094e
/*******************************************************************************
* 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, &params);
/* create return python object */
PyObject* obj = pytools_return(&params, class_swift_params);
return obj;
}
src/parser_wrapper.h deleted 100644 → 0
+ 0
− 36
View file @ a05c094e
/*******************************************************************************
* 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__
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment