Skip to content
Snippets Groups Projects
Commit 743713e3 authored by lhausamm's avatar lhausamm
Browse files

Read config.h in makeInput.py

parent ba230049
No related branches found
No related tags found
1 merge request!486Add chemistry in part
......@@ -19,6 +19,7 @@
import h5py
from numpy import *
import re as regex
# Generates a swift IC file containing a cartesian distribution of particles
# at a constant density and pressure in a cubic box
......@@ -32,12 +33,24 @@ P = 1. # Pressure
gamma = 5./3. # Gas adiabatic index
fileName = "input.hdf5"
# read flags from config.h
flags = {}
config_file = "../config.h"
with open(config_file, "r") as f:
reg = "#define (\w*) (.*)\n"
for line in f:
m = regex.match(reg, line)
if m:
flags[m.group(1)] = m.group(2)
cooling_grackle = "COOLING_GRACKLE" in flags
#---------------------------------------------------
numPart = L**3
mass = boxSize**3 * rho / numPart
internalEnergy = P / ((gamma - 1.)*rho)
he_density = rho * 0.24
if cooling_grackle:
he_density = rho * 0.24
#Generate particles
coords = zeros((numPart, 3))
......@@ -46,7 +59,8 @@ m = zeros((numPart, 1))
h = zeros((numPart, 1))
u = zeros((numPart, 1))
ids = zeros((numPart, 1), dtype='L')
he = zeros((numPart, 1))
if cooling_grackle:
he = zeros((numPart, 1))
for i in range(L):
for j in range(L):
......@@ -65,7 +79,8 @@ for i in range(L):
h[index] = 2.251 * boxSize / L
u[index] = internalEnergy
ids[index] = index
he[index] = he_density
if cooling_grackle:
he[index] = he_density
#--------------------------------------------------
......@@ -112,7 +127,8 @@ ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
ds[()] = u
ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
ds[()] = ids
ds = grp.create_dataset('HeDensity', (numPart, 1), 'f')
ds[()] = he
if cooling_grackle:
ds = grp.create_dataset('HeDensity', (numPart, 1), 'f')
ds[()] = he
file.close()
......@@ -42,7 +42,9 @@ int main() {
const double boxSize = 1.;
const size_t L = 4;
const double rho = 2.;
#ifdef COOLING_GRACKLE
const float He_density = rho * 0.24;
#endif
/* Read data */
read_ic_single("input.hdf5", &us, dim, &parts, &gparts, &sparts, &Ngas,
......@@ -95,9 +97,6 @@ int main() {
assert(parts[n].a_hydro[2] == 0.);
#ifdef COOLING_GRACKLE
printf("Test %g, %g, %g\n", parts[n].cooling_data.He_density,
He_density,
parts[n].cooling_data.He_density-He_density);
assert(parts[n].cooling_data.He_density == He_density);
#endif
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment