Commit 94ca7df5 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Blackholes swallowing - Gas case

parent 3c4036b8
......@@ -30,6 +30,9 @@ Snapshots:
delta_time: 1.01 # Time difference between consecutive outputs (in internal units)
compression: 1
Scheduler:
links_per_tasks: 500
# Parameters governing the conserved quantities statistics
Statistics:
scale_factor_first: 0.92 # Scale-factor of the first stat dump (cosmological run)
......
......@@ -22,7 +22,10 @@ Cosmology:
Omega_m: 0.307 # Matter density parameter
Omega_lambda: 0.693 # Dark-energy density parameter
Omega_b: 0.0455 # Baryon density parameter
Scheduler:
links_per_tasks: 500
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
......
......@@ -25,6 +25,7 @@ Cosmology:
Scheduler:
max_top_level_cells: 8
links_per_tasks: 500
# Parameters governing the time integration
TimeIntegration:
......
import h5py as h5
import numpy as np
f = h5.File("bh_swallowing_serial_0001.hdf5", "r")
f_mpi = h5.File("bh_swallowing_mpi_0001.hdf5", "r")
f_ics = h5.File("bh_swallowing.hdf5", "r")
m = np.array(f["/PartType5/Masses"])#[:]
m_mpi = np.array(f_mpi["/PartType5/Masses"])#[:]
m_ics = np.array(f_ics["/PartType5/Masses"])#[:]
ids = np.array(f["/PartType5/ParticleIDs"])#[:]
ids_mpi = np.array(f_mpi["/PartType5/ParticleIDs"])#[:]
ids_ics = np.array(f_ics["/PartType5/ParticleIDs"])#[:]
rank = np.argsort(ids)
rank_mpi = np.argsort(ids_mpi)
rank_ics = np.argsort(ids_ics)
m = m[rank]
m_mpi = m_mpi[rank_mpi]
m_ics = m_ics[rank_ics]
ids = ids[rank]
ids_mpi = ids_mpi[rank_mpi]
ids_ics = ids_ics[rank_ics]
# Check
#print ids - ids_mpi
#print ids - ids_ics
count_wrong = 0
print "ID -- ICs mass -- serial mass -- MPI mass"
for i in range(np.size(ids)):
print "%14d %5f %5f %5f"%(ids[i], m_ics[i], m[i], m_mpi[i]),
if m[i] > m_ics[i] * 1.1:
print "#",
else:
print " ",
if m[i] != m_mpi[i]:
count_wrong += 1
print " <---"
else:
print ""
print "Found", count_wrong, "wrong BH masses."
ids_gas_mpi = f_mpi["/PartType0/ParticleIDs"][:]
ids_removed_mpi = np.loadtxt("list_mpi.txt")
for i in range(np.size(ids_removed_mpi)):
result = np.where(ids_gas_mpi == ids_removed_mpi)
print result
ids_gas = f["/PartType0/ParticleIDs"][:]
ids_removed = np.loadtxt("list.txt")
for i in range(np.size(ids_removed)):
result = np.where(ids_gas == ids_removed)
print result
#rho_gas = f["/PartType0/Density"][:]
#print np.mean(rho_gas), np.std(rho_gas)
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5
###############################################################################
# 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
from numpy import *
# Some constants
solar_mass_cgs = 1.988480e33
kpc_in_cm = 3.085678e21
mp_cgs = 1.67e-24
boltzmann_k_cgs = 1.38e-16
# Parameters
gamma = 5./3. # Gas adiabatic index
rho_cgs = mp_cgs # Background density
u0_cgs = 1.2e12 # Desired initial internal energy (1.2e12 ~ 10^4K)
P_cgs = rho_cgs*u0_cgs*(gamma - 1.) # Background pressure
fileName = "bh_swallowing.hdf5"
# Units
unit_l_cgs = 3.085678e24 # kpc
unit_m_cgs = 1.988480e43 # 10^10 Msun
unit_v_cgs = 1e5 # km / s
unit_A_cgs = 1.
unit_T_cgs = 1.
unit_t_cgs = unit_l_cgs / unit_v_cgs
boxsize_cgs = 10. * kpc_in_cm
vol_cgs = boxsize_cgs**3
#---------------------------------------------------
glass = h5py.File("glassCube_32.hdf5", "r")
# Read particle positions and h from the glass
pos = glass["/PartType0/Coordinates"][:,:]
eps = 1e-6
pos = (pos - pos.min()) / (pos.max() - pos.min() + eps) * boxsize_cgs / unit_l_cgs
h = glass["/PartType0/SmoothingLength"][:] * boxsize_cgs / unit_l_cgs
numPart = size(h)
# Generate extra arrays
v = zeros((numPart, 3))
ids = linspace(1, numPart, numPart)
m_cgs = zeros(numPart)
u_cgs = zeros(numPart)
m = zeros(numPart)
u = zeros(numPart)
m_cgs[:] = rho_cgs * vol_cgs / numPart
u_cgs[:] = P_cgs / (rho_cgs * (gamma - 1))
# BHs
bh_pos = zeros((2, 3))
bh_pos[0,:] = 0.5 * boxsize_cgs / unit_l_cgs
#diff = [4.812127e-03, 4.908179e-03, -4.878537e-03] - bh_pos[0,:]
#print diff
bh_pos[1,:] = 0.5 * boxsize_cgs / unit_l_cgs + diff
print bh_pos
# BHs don't move
bh_v = zeros((2, 3))
bh_v[:,:] = 0.
# increase mass to keep it at center
bh_m_cgs = ones(2) * m_cgs[0]
bh_ids = linspace(numPart + 1, numPart + 2, 2, dtype='L')
bh_h = ones(2) * [h.max()]
print bh_ids
#--------------------------------------------------
# Check quantities are correct for debugging
print("boxsize kpc " + str(boxsize_cgs/kpc_in_cm))
print("density cm^-3 " + str(rho_cgs/mp_cgs))
print("initial temperature K " + str(u_cgs[0] / boltzmann_k_cgs*((gamma - 1)*rho_cgs)))
# Convert to internal units
bh_m = bh_m_cgs/unit_m_cgs
m[:] = m_cgs/unit_m_cgs
u[:] = u_cgs*unit_v_cgs**-2
boxsize = boxsize_cgs/unit_l_cgs
#--------------------------------------------------
#File
file = h5py.File(fileName, 'w')
# Header
grp = file.create_group("/Header")
grp.attrs["BoxSize"] = [boxsize]*3
grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 2]
#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, 2]
#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
grp.attrs["Dimension"] = 3
#Runtime parameters
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = 0
#Units
grp = file.create_group("/Units")
grp.attrs["Unit length in cgs (U_L)"] = unit_l_cgs
grp.attrs["Unit mass in cgs (U_M)"] = unit_m_cgs
grp.attrs["Unit time in cgs (U_t)"] = unit_t_cgs
grp.attrs["Unit current in cgs (U_I)"] = unit_A_cgs
grp.attrs["Unit temperature in cgs (U_T)"] = unit_T_cgs
#Particle group
grp = file.create_group("/PartType0")
grp.create_dataset('Coordinates', data=pos, dtype='d')
grp.create_dataset('Velocities', data=v, dtype='f')
grp.create_dataset('Masses', data=m, dtype='f')
grp.create_dataset('SmoothingLength', data=h, dtype='f')
grp.create_dataset('InternalEnergy', data=u, dtype='f')
grp.create_dataset('ParticleIDs', data=ids, dtype='L')
# stellar group
grp = file.create_group("/PartType5")
grp.create_dataset("Coordinates", (2,3),data=bh_pos, dtype="d")
grp.create_dataset('Velocities', (2,3),data=bh_v, dtype='f')
grp.create_dataset('Masses', (2,1),data=bh_m, dtype='f')
grp.create_dataset('SmoothingLength', (2,1),data=bh_h, dtype='f')
grp.create_dataset('ParticleIDs', (2,1),data=bh_ids[::-1], dtype='L')
file.close()
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun in grams
UnitLength_in_cgs: 3.08567758e24 # Mpc in centimeters
UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second
UnitCurrent_in_cgs: 1. # Amperes
UnitTemp_in_cgs: 1. # Kelvin
# Cosmological parameters
Cosmology:
h: 0.6777 # Reduced Hubble constant
a_begin: 0.5 # 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
# Parameters governing the time integration
TimeIntegration:
time_begin: 0 # The starting time of the simulation (in internal units).
time_end: 1.3e-2 # The end time of the simulation (in internal units).
dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-5 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: bh_swallowing_serial # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 3.e-5 # Time difference between consecutive outputs without cosmology (internal units)
compression: 1
# Parameters governing the conserved quantities statistics
Statistics:
time_first: 0.
delta_time: 1.e-5 # non cosmology time between statistics output
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
minimal_temperature: 10. # Kelvin
# Parameters related to the initial conditions
InitialConditions:
file_name: ./bh_swallowing.hdf5 # The file to read
periodic: 1
Scheduler:
max_top_level_cells: 8
tasks_per_cell: 500
# Parameters for the EAGLE "equation of state"
EAGLEEntropyFloor:
Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin.
Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor
Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
# Metallicites read in for the gas and star
EAGLEChemistry:
init_abundance_metal: 0.01
init_abundance_Hydrogen: 0.752
init_abundance_Helium: 0.248
init_abundance_Carbon: 0.0
init_abundance_Nitrogen: 0.0
init_abundance_Oxygen: 0.0
init_abundance_Neon: 0.0
init_abundance_Magnesium: 0.0
init_abundance_Silicon: 0.0
init_abundance_Iron: 0.0
# Standard EAGLE cooling options
EAGLECooling:
dir_name: ./coolingtables/ # Location of the Wiersma+08 cooling tables
H_reion_z: 11.5 # Redshift of Hydrogen re-ionization
H_reion_eV_p_H: 2.0
He_reion_z_centre: 3.5 # Redshift of the centre of the Helium re-ionization Gaussian
He_reion_z_sigma: 0.5 # Spread in redshift of the Helium re-ionization Gaussian
He_reion_eV_p_H: 2.0 # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
# EAGLE AGN model
EAGLEAGN:
subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.98848e43 # 10^10 M_sun in grams
UnitLength_in_cgs: 3.08567758e24 # Mpc in centimeters
UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second
UnitCurrent_in_cgs: 1. # Amperes
UnitTemp_in_cgs: 1. # Kelvin
# Cosmological parameters
Cosmology:
h: 0.6777 # Reduced Hubble constant
a_begin: 0.5 # 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
# Parameters governing the time integration
TimeIntegration:
time_begin: 0 # The starting time of the simulation (in internal units).
time_end: 1.3e-2 # The end time of the simulation (in internal units).
dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-5 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: bh_swallowing_mpi # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 3.e-5 # Time difference between consecutive outputs without cosmology (internal units)
compression: 1
# Parameters governing the conserved quantities statistics
Statistics:
time_first: 0.
delta_time: 1.e-5 # non cosmology time between statistics output
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
minimal_temperature: 10. # Kelvin
# Parameters related to the initial conditions
InitialConditions:
file_name: ./bh_swallowing.hdf5 # The file to read
periodic: 1
Scheduler:
max_top_level_cells: 8
tasks_per_cell: 500
# Parameters for the EAGLE "equation of state"
EAGLEEntropyFloor:
Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin.
Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor
Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
# Metallicites read in for the gas and star
EAGLEChemistry:
init_abundance_metal: 0.01
init_abundance_Hydrogen: 0.752
init_abundance_Helium: 0.248
init_abundance_Carbon: 0.0
init_abundance_Nitrogen: 0.0
init_abundance_Oxygen: 0.0
init_abundance_Neon: 0.0
init_abundance_Magnesium: 0.0
init_abundance_Silicon: 0.0
init_abundance_Iron: 0.0
# Standard EAGLE cooling options
EAGLECooling:
dir_name: ./coolingtables/ # Location of the Wiersma+08 cooling tables
H_reion_z: 11.5 # Redshift of Hydrogen re-ionization
H_reion_eV_p_H: 2.0
He_reion_z_centre: 3.5 # Redshift of the centre of the Helium re-ionization Gaussian
He_reion_z_sigma: 0.5 # Spread in redshift of the Helium re-ionization Gaussian
He_reion_eV_p_H: 2.0 # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
# EAGLE AGN model
EAGLEAGN:
subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
......@@ -53,7 +53,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
star_formation_struct.h star_formation.h star_formation_iact.h \
star_formation_logger.h star_formation_logger_struct.h \
velociraptor_struct.h velociraptor_io.h random.h memuse.h black_holes.h black_holes_io.h \
black_holes_properties.h feedback.h feedback_struct.h feedback_properties.h
black_holes_properties.h black_holes_struct.h feedback.h feedback_struct.h feedback_properties.h
# source files for EAGLE cooling
EAGLE_COOLING_SOURCES =
......@@ -218,6 +218,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
black_holes/Default/black_holes.h black_holes/Default/black_holes_io.h \
black_holes/Default/black_holes_part.h black_holes/Default/black_holes_iact.h \
black_holes/Default/black_holes_properties.h \
black_holes/Default/black_holes_struct.h \
black_holes/EAGLE/black_holes.h black_holes/EAGLE/black_holes_io.h \
black_holes/EAGLE/black_holes_part.h black_holes/EAGLE/black_holes_iact.h \
black_holes/EAGLE/black_holes_properties.h
......
......@@ -23,6 +23,7 @@
/* Local includes */
#include "black_holes_properties.h"
#include "black_holes_struct.h"
#include "dimension.h"
#include "kernel_hydro.h"
#include "minmax.h"
......@@ -139,6 +140,58 @@ black_holes_bpart_has_no_neighbours(struct bpart* restrict bp,
bp->density.wcount_dh = 0.f;
}
/**
* @brief Update the properties of a black hole particles by swallowing
* a gas particle.
*
* @param bp The #bpart to update.
* @param p The #part that is swallowed.
* @param xp The #xpart that is swallowed.
* @param cosmo The current cosmological model.
*/
__attribute__((always_inline)) INLINE static void black_holes_swallow_part(
struct bpart* bp, const struct part* p, const struct xpart* xp,
const struct cosmology* cosmo) {
/* Nothing to do here: No swallowing in the default model */
}
/**
* @brief Update a given #part's BH data field to mark the particle has
* not yet been swallowed.
*
* @param p_data The #part's #black_holes_part_data structure.
*/
__attribute__((always_inline)) INLINE static void
black_holes_mark_as_not_swallowed(struct black_holes_part_data* p_data) {
/* Nothing to do here: No swallowing in the default model */
}
/**
* @brief Update a given #part's BH data field to mark the particle has
* having been been swallowed.
*
* @param p_data The #part's #black_holes_part_data structure.
*/
__attribute__((always_inline)) INLINE static void black_holes_mark_as_swallowed(
struct black_holes_part_data* p_data) {
/* Nothing to do here: No swallowing in the default model */
}
/**
* @brief Return the ID of the BH that should swallow this #part.
*
* @param p_data The #part's #black_holes_part_data structure.
*/
__attribute__((always_inline)) INLINE static long long
black_holes_get_swallow_id(struct black_holes_part_data* p_data) {
/* Return a non-existing ID */
return -1;
}
/**
* @brief Compute the accretion rate of the black hole and all the quantites
* required for the feedback loop.
......
......@@ -63,6 +63,28 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_bh_density(
#endif
}
/**
* @brief Swallowing interaction between two particles (non-symmetric).
*
* Function used to flag the gas particles that will be swallowed
* by the black hole particle.
*
* @param r2 Comoving square distance between the two particles.
* @param dx Comoving vector separating both particles (pi - pj).
* @param hi Comoving smoothing-length of particle i.
* @param hj Comoving smoothing-length of particle j.
* @param bi First particle (black hole).
* @param pj Second particle (gas)
* @param xpj The extended data of the second particle.
* @param cosmo The cosmological model.
* @param ti_current Current integer time value (for random numbers).
*/
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_bh_swallow(
const float r2, const float *dx, const float hi, const float hj,
struct bpart *restrict bi, struct part *restrict pj,
struct xpart *restrict xpj, const struct cosmology *cosmo,
const integertime_t ti_current) {}
/**
* @brief Feedback interaction between two particles (non-symmetric).
*
......
......@@ -19,8 +19,7 @@
#ifndef SWIFT_DEFAULT_BLACK_HOLE_PART_H
#define SWIFT_DEFAULT_BLACK_HOLE_PART_H
/* Some standard headers. */
#include <stdlib.h>
#include "chemistry_struct.h"
/**
* @brief Particle fields for the black hole particles.
......@@ -63,6 +62,10 @@ struct bpart {
} density;
/*! Chemistry information (e.g. metal content at birth, swallowed metal
* content, etc.) */
struct chemistry_bpart_data chemistry_data;
#ifdef SWIFT_DEBUG_CHECKS
/* Time of the last drift */
......
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl)
*
* 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.