Commit 6b94cbd5 authored by Stefan Arridge's avatar Stefan Arridge
Browse files

Added Cooling Box example

parent d76dbba0
import numpy as np
import matplotlib.pyplot as plt
import h5py as h5
import sys
filename = "./energy.txt"
plot_dir = "./"
cooling_lambda = 1.0e-22
n_H = 0.1
T_init = 1.0e5
T_floor = 1.0e4
plot_title = r"$\Lambda \, = \, %1.1g \mathrm{erg}\mathrm{cm^3}\mathrm{s^{-1}} \, \, T_{init} = %1.1g\mathrm{K} \, \, T_{floor} = %1.1g\mathrm{K} \, \, n_H = %1.1g\mathrm{cm^{-3}}$" %(cooling_lambda,T_init,T_floor,n_H)
plot_filename = "energy_plot_creasey_no_cooling_T_init_1p0e5_n_H_0p1.png"
#analytic_solution = np.zeros(n_snaps-1)
array = np.genfromtxt(filename,skip_header = 1)
time = array[:,0]
total_energy = array[:,2]
# plt.plot(e_kin_array,label = "Kinetic Energy")
# plt.plot(e_therm_array,label = "Internal Energy")
plt.plot(time,total_energy,'r',label = "Total Energy")
#plt.plot(time,analytic_solution,'--k',label = "Analytic Solution")
plt.xlabel("Time (code units)")
plt.ylabel("Energy (code units)")
#plt.ylim((0,1.0e8))
plt.legend(loc = "upper right")
#plt.title(plot_title)
full_plot_filename = "%s/%s" %(plot_dir,plot_filename)
if (int(sys.argv[1])==0):
plt.show()
else:
plt.savefig(full_plot_filename,format = "png")
plt.close()
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# 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/>.
#
##############################################################################
import h5py
import sys
from numpy import *
# Generates a swift IC file containing a cartesian distribution of particles
# at a constant density and pressure in a cubic box
# Parameters
periodic= 1 # 1 For periodic box
boxSize = 1 #1 kiloparsec
L = int(sys.argv[1]) # Number of particles along one axis
rho = 3.2e3 # Density in code units (0.1 hydrogen atoms per cm^3)
P = 4.5e7 # Pressure in code units (at 10^5K)
gamma = 5./3. # Gas adiabatic index
eta = 1.2349 # 48 ngbs with cubic spline kernel
fileName = "uniformBox.hdf5"
#---------------------------------------------------
numPart = L**3
mass = boxSize**3 * rho / numPart
print mass
internalEnergy = P / ((gamma - 1.)*rho)
#--------------------------------------------------
#File
file = h5py.File(fileName, 'w')
# Header
grp = file.create_group("/Header")
grp.attrs["BoxSize"] = boxSize
grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
grp.attrs["Time"] = 0.0
grp.attrs["NumFilesPerSnapshot"] = 1
grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
grp.attrs["Flag_Entropy_ICs"] = 0
#Runtime parameters
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
#Units
grp = file.create_group("/Units")
grp.attrs["Unit length in cgs (U_L)"] = 3.08e21
grp.attrs["Unit mass in cgs (U_M)"] = 2.0e33
grp.attrs["Unit time in cgs (U_t)"] = 3.08e16
grp.attrs["Unit current in cgs (U_I)"] = 1.
grp.attrs["Unit temperature in cgs (U_T)"] = 1.
#Particle group
grp = file.create_group("/PartType0")
v = zeros((numPart, 3))
ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
ds[()] = v
v = zeros(1)
m = full((numPart, 1), mass)
ds = grp.create_dataset('Masses', (numPart,1), 'f')
ds[()] = m
m = zeros(1)
h = full((numPart, 1), eta * boxSize / L)
ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
ds[()] = h
h = zeros(1)
u = full((numPart, 1), internalEnergy)
ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
ds[()] = u
u = zeros(1)
ids = linspace(0, numPart, numPart, endpoint=False).reshape((numPart,1))
ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
ds[()] = ids + 1
x = ids % L;
y = ((ids - x) / L) % L;
z = (ids - x - L * y) / L**2;
coords = zeros((numPart, 3))
coords[:,0] = z[:,0] * boxSize / L + boxSize / (2*L)
coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L)
coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L)
ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
ds[()] = coords
file.close()
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2016 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/>.
#
##############################################################################
import h5py
import sys
from numpy import *
# Generates a swift IC file containing a cartesian distribution of particles
# at a constant density and pressure in a cubic box
# Parameters
periodic= 1 # 1 For periodic box
boxSize = 1.
L = int(sys.argv[1]) # Number of particles along one axis
N = int(sys.argv[2]) # Write N particles at a time to avoid requiring a lot of RAM
rho = 2. # Density
P = 1. # Pressure
gamma = 5./3. # Gas adiabatic index
eta = 1.2349 # 48 ngbs with cubic spline kernel
fileName = "uniformBox_%d.hdf5"%L
#---------------------------------------------------
numPart = L**3
mass = boxSize**3 * rho / numPart
internalEnergy = P / ((gamma - 1.)*rho)
#---------------------------------------------------
n_iterations = numPart / N
remainder = numPart % N
print "Writing", numPart, "in", n_iterations, "iterations of size", N, "and a remainder of", remainder
#---------------------------------------------------
#File
file = h5py.File(fileName, 'w')
# Header
grp = file.create_group("/Header")
grp.attrs["BoxSize"] = boxSize
grp.attrs["NumPart_Total"] = [numPart % (long(1)<<32), 0, 0, 0, 0, 0]
grp.attrs["NumPart_Total_HighWord"] = [numPart / (long(1)<<32), 0, 0, 0, 0, 0]
grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
grp.attrs["Time"] = 0.0
grp.attrs["NumFilesPerSnapshot"] = 1
grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
grp.attrs["Flag_Entropy_ICs"] = 0
#Runtime parameters
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
#Units
grp = file.create_group("/Units")
grp.attrs["Unit length in cgs (U_L)"] = 1.
grp.attrs["Unit mass in cgs (U_M)"] = 1.
grp.attrs["Unit time in cgs (U_t)"] = 1.
grp.attrs["Unit current in cgs (U_I)"] = 1.
grp.attrs["Unit temperature in cgs (U_T)"] = 1.
#Particle group
grp = file.create_group("/PartType0")
# First create the arrays in the file
ds_v = grp.create_dataset('Velocities', (numPart, 3), 'f', chunks=True, compression="gzip")
ds_m = grp.create_dataset('Masses', (numPart,1), 'f', chunks=True, compression="gzip")
ds_h = grp.create_dataset('SmoothingLength', (numPart,1), 'f', chunks=True, compression="gzip")
ds_u = grp.create_dataset('InternalEnergy', (numPart,1), 'f', chunks=True, compression="gzip")
ds_id = grp.create_dataset('ParticleIDs', (numPart, 1), 'L', chunks=True, compression="gzip")
ds_x = grp.create_dataset('Coordinates', (numPart, 3), 'd', chunks=True, compression="gzip")
# Now loop and create parts of the dataset
offset = 0
for n in range(n_iterations):
v = zeros((N, 3))
ds_v[offset:offset+N,:] = v
v = zeros(1)
m = full((N, 1), mass)
ds_m[offset:offset+N] = m
m = zeros(1)
h = full((N, 1), eta * boxSize / L)
ds_h[offset:offset+N] = h
h = zeros(1)
u = full((N, 1), internalEnergy)
ds_u[offset:offset+N] = u
u = zeros(1)
ids = linspace(offset, offset+N, N, endpoint=False).reshape((N,1))
ds_id[offset:offset+N] = ids + 1
x = ids % L;
y = ((ids - x) / L) % L;
z = (ids - x - L * y) / L**2;
ids = zeros(1)
coords = zeros((N, 3))
coords[:,0] = z[:,0] * boxSize / L + boxSize / (2*L)
coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L)
coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L)
ds_x[offset:offset+N,:] = coords
coords = zeros((1,3))
offset += N
print "Done", offset,"/", numPart, "(%.1f %%)"%(100*(float)(offset)/numPart)
# And now, the remainder
v = zeros((remainder, 3))
ds_v[offset:offset+remainder,:] = v
v = zeros(1)
m = full((remainder, 1), mass)
ds_m[offset:offset+remainder] = m
m = zeros(1)
h = full((remainder, 1), eta * boxSize / L)
ds_h[offset:offset+remainder] = h
h = zeros(1)
u = full((remainder, 1), internalEnergy)
ds_u[offset:offset+remainder] = u
u = zeros(1)
ids = linspace(offset, offset+remainder, remainder, endpoint=False).reshape((remainder,1))
ds_id[offset:offset+remainder] = ids + 1
x = ids % L;
y = ((ids - x) / L) % L;
z = (ids - x - L * y) / L**2;
coords = zeros((remainder, 3))
coords[:,0] = z[:,0] * boxSize / L + boxSize / (2*L)
coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L)
coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L)
ods_x[offset:offset+remainder,:] = coords
print "Done", offset+remainder,"/", numPart
file.close()
#!/bin/bash
# Generate the initial conditions if they are not present.
echo "Generating initial conditions for the uniform box example..."
python makeIC.py 10
../swift -s -C -t 16 uniformBox.yml
python energy_plot.py 0
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
......@@ -59,9 +60,11 @@ void cooling_init(const struct swift_params* parameter_file,
double T_min = cooling->creasey_cooling.min_temperature;
double mu = cooling->creasey_cooling.mean_molecular_weight;
double m_p = phys_const->const_proton_mass;
cooling->creasey_cooling.min_internal_energy = k_b * T_min / (hydro_gamma_minus_one * mu * m_p);
double u_floor = k_b * T_min / (hydro_gamma_minus_one * mu * m_p);
double u_floor_cgs = u_floor * units_cgs_conversion_factor(us,UNIT_CONV_ENERGY_PER_UNIT_MASS);
cooling->creasey_cooling.min_internal_energy = u_floor;
cooling->creasey_cooling.min_internal_energy_cgs = u_floor_cgs;
#endif /* CREASEY_COOLING */
}
......@@ -95,11 +98,11 @@ void update_entropy(const struct phys_const* const phys_const, const struct Unit
const struct cooling_data* cooling, struct part* p, double dt){
/*updates the entropy of a particle after integrating the cooling equation*/
float u_old;
float u_new;
float new_entropy;
float old_entropy = p->entropy;
float rho = p->rho;
double u_old;
double u_new;
double new_entropy;
double old_entropy = p->entropy;
double rho = p->rho;
// u_old = old_entropy/(GAMMA_MINUS1) * pow(rho,GAMMA_MINUS1);
u_old = hydro_get_internal_energy(p,0); // dt = 0 because using current entropy
......@@ -111,15 +114,15 @@ void update_entropy(const struct phys_const* const phys_const, const struct Unit
/*This function integrates the cooling equation, given the initial
thermal energy, density and the timestep dt. Returns the final internal energy*/
float calculate_new_thermal_energy(float u_old, float rho, float dt,
double calculate_new_thermal_energy(double u_old, double rho, double dt,
const struct cooling_data* cooling,
const struct phys_const* const phys_const,
const struct UnitSystem* us){
#ifdef CONST_COOLING
/*du/dt = -lambda, independent of density*/
float du_dt = -cooling->const_cooling.lambda;
float u_floor = cooling->const_cooling.min_energy;
float u_new;
double du_dt = -cooling->const_cooling.lambda;
double u_floor = cooling->const_cooling.min_energy;
double u_new;
if (u_old - du_dt*dt > u_floor){
u_new = u_old + du_dt*dt;
}
......@@ -130,23 +133,22 @@ float calculate_new_thermal_energy(float u_old, float rho, float dt,
#ifdef CREASEY_COOLING
/* rho*du/dt = -lambda*n_H^2 */
float u_new;
float m_p = phys_const->const_proton_mass;
float X_H = cooling->creasey_cooling.hydrogen_mass_abundance;
float lambda_cgs = cooling->creasey_cooling.lambda; //this is always in cgs
float u_floor = cooling->creasey_cooling.min_internal_energy;
double u_new;
double m_p = phys_const->const_proton_mass;
double X_H = cooling->creasey_cooling.hydrogen_mass_abundance;
double lambda_cgs = cooling->creasey_cooling.lambda; //this is always in cgs
double u_floor_cgs = cooling->creasey_cooling.min_internal_energy_cgs;
/*convert from internal code units to cgs*/
float dt_cgs = dt * units_cgs_conversion_factor(us,UNIT_CONV_TIME);
float rho_cgs = rho * units_cgs_conversion_factor(us,UNIT_CONV_DENSITY);
float m_p_cgs = m_p * units_cgs_conversion_factor(us,UNIT_CONV_MASS);
float n_H_cgs = rho_cgs / m_p_cgs;
float u_old_cgs = u_old * units_cgs_conversion_factor(us,UNIT_CONV_ENERGY_PER_UNIT_MASS);
float u_floor_cgs = u_floor * units_cgs_conversion_factor(us,UNIT_CONV_ENERGY_PER_UNIT_MASS);
float du_dt_cgs = -lambda_cgs * n_H_cgs * n_H_cgs / rho_cgs;
float u_new_cgs;
if (u_old_cgs + du_dt_cgs*dt_cgs > u_floor_cgs){
double dt_cgs = dt * units_cgs_conversion_factor(us,UNIT_CONV_TIME);
double rho_cgs = rho * units_cgs_conversion_factor(us,UNIT_CONV_DENSITY);
double m_p_cgs = m_p * units_cgs_conversion_factor(us,UNIT_CONV_MASS);
double n_H_cgs = X_H * rho_cgs / m_p_cgs;
double u_old_cgs = u_old * units_cgs_conversion_factor(us,UNIT_CONV_ENERGY_PER_UNIT_MASS);
double du_dt_cgs = -lambda_cgs * n_H_cgs * n_H_cgs / rho_cgs;
double u_new_cgs;
if (u_old_cgs + du_dt_cgs * dt_cgs > u_floor_cgs){
u_new_cgs = u_old_cgs + du_dt_cgs*dt_cgs;
}
else{
......
......@@ -56,6 +56,7 @@ struct cooling_data {
double hydrogen_mass_abundance;
double mean_molecular_weight;
double min_internal_energy;
double min_internal_energy_cgs;
double cooling_tstep_mult;
} creasey_cooling;
#endif
......@@ -71,13 +72,13 @@ struct cooling_data {
* @param phys_const The physical constants in internal units.
* @param Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float
__attribute__((always_inline)) INLINE static double
cooling_timestep(const struct cooling_data* cooling,
const struct phys_const* const phys_const,
const struct part* const p) {
const float cooling_rate = cooling->const_cooling.lambda;
const float internal_energy = hydro_get_internal_energy(p,0);// dt = 0 because using current entropy
const double cooling_rate = cooling->const_cooling.lambda;
const double internal_energy = hydro_get_internal_energy(p,0);// dt = 0 because using current entropy
return cooling->const_cooling.cooling_tstep_mult * internal_energy / cooling_rate;
}
......@@ -93,7 +94,7 @@ void cooling_init(const struct swift_params* parameter_file,
void cooling_print(const struct cooling_data* cooling);
void update_entropy(const struct phys_const* const phys_const, const struct UnitSystem* us,
const struct cooling_data* cooling, struct part* p, double dt);
float calculate_new_thermal_energy(float u_old, float rho, float dt,
double calculate_new_thermal_energy(double u_old, double rho, double dt,
const struct cooling_data* cooling,
const struct phys_const* const phys_const,
const struct UnitSystem* us);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment