Commit ca675a67 authored by Rainer Weinberger's avatar Rainer Weinberger

added example isolated_galaxy_collisionless_3d

parent d6ca5099
#!/bin/bash # this line only there to enable syntax highlighting in this file
#########################################################
# Enable/Disable compile-time options as needed #
# examples/isolated_galaxy_collisionless_3d/Config.sh #
#########################################################
#--------------------------------------- Gravity treatment
SELFGRAVITY # gravitational intraction between simulation particles/cells
HIERARCHICAL_GRAVITY # use hierarchical splitting of the time integration of the gravity
CELL_CENTER_GRAVITY # uses geometric centers to calculate gravity of cells, only possible with HIERARCHICAL_GRAVITY
ALLOW_DIRECT_SUMMATION # Performed direct summation instead of tree-based gravity if number of active particles < DIRECT_SUMMATION_THRESHOLD (= 3000 unless specified differently here)
DIRECT_SUMMATION_THRESHOLD=500 # Overrides maximum number of active particles for which direct summation is performed instead of tree based calculation
GRAVITY_NOT_PERIODIC # gravity is not treated periodically
RANDOMIZE_DOMAINCENTER # random displacement to position of domain center; avoids correlated force-errors, important mainly for isolated systems (which otherwise might start to drift in some direction).
#--------------------------------------- Gravity softening
NSOFTTYPES=2 # Number of different softening values to which particle types can be mapped.
MULTIPLE_NODE_SOFTENING # If a tree node is to be used which is softened, this is done with the softenings of its different mass components
INDIVIDUAL_GRAVITY_SOFTENING=32 # bitmask with particle types where the softenig type should be chosen with that of parttype 1 as a reference type
ADAPTIVE_HYDRO_SOFTENING # Adaptive softening of gas cells depending on their size
#--------------------------------------- Time integration options
TREE_BASED_TIMESTEPS # non-local timestep criterion (take 'signal speed' into account)
#--------------------------------------- Single/Double Precision
DOUBLEPRECISION=1 # Mode of double precision: not defined: single; 1: full double precision 2: mixed, 3: mixed, fewer single precisions; unless short of memory, use 1.
NGB_TREE_DOUBLEPRECISION # if this is enabled, double precision is used for the neighbor node extension
#--------------------------------------- Things for special behaviour
PROCESS_TIMES_OF_OUTPUTLIST # goes through times of output list prior to starting the simulaiton to ensure that outputs are written as close to the desired time as possible (as opposed to at next possible time if this flag is not active)
#--------------------------------------- Output/Input options
HAVE_HDF5 # needed when HDF5 I/O support is desired (recommended)
#--------------------------------------- Testing and Debugging options
DEBUG # enables core-dumps
EXEC = GalIC
CONFIG = Config.sh
BUILD_DIR = build
SRC_DIR = src
#PARAMFILE = Model_D3.param
#N := 16
ifdef SYSTYPE
SYSTYPE := "$(SYSTYPE)"
-include Makefile.systype
else
include Makefile.systype
endif
MAKEFILES = Makefile config-makefile
ifeq ($(wildcard Makefile.systype), Makefile.systype)
MAKEFILES += Makefile.systype
endif
PERL = /usr/bin/perl
RESULT := $(shell CONFIG=$(CONFIG) PERL=$(PERL) BUILD_DIR=$(BUILD_DIR) make -f config-makefile)
CONFIGVARS := $(shell cat $(BUILD_DIR)/galicconfig.h)
#MPICHLIB = -lmpich
GMPLIB = -lgmp
GSLLIB = -lgsl -lgslcblas
MATHLIB = -lm
###############################
# Determine your SYSTEM here #
###############################
ifeq ($(SYSTYPE),"Darwin")
CC = mpicc # sets the C-compiler
OPTIMIZE = -ggdb -O3 -Wall -Wno-format-security -Wno-unknown-pragmas
ifeq (NUM_THREADS,$(findstring NUM_THREADS,$(CONFIGVARS)))
OPTIMIZE += -fopenmp
MPI_COMPILE_FLAGS = $(shell mpicc --showme:compile)
CC = gcc $(MPI_COMPILE_FLAGS) # to replace clang with gcc (mpicc uses clang for some reason)
endif
GSL_INCL = -I/sw/include -I/opt/local/include
GSL_LIBS = -L/sw/lib -L/opt/local/lib
FFTW_INCL= -I/sw/include -I/opt/local/include
FFTW_LIBS= -L/sw/lib -L/opt/local/lib
MPICHLIB = -lmpi
HDF5INCL = -I/sw/lib -I/opt/local/include -DH5_USE_16_API #-DUSE_SSE
HDF5LIB = -L/sw/lib -L/opt/local/lib -lhdf5 -lz
endif
ifeq ($(SYSTYPE),"APHI")
CC = mpicc
CXX = mpicxx
#OPTIMIZE = -g -w -m64 -O3 -msse3
OPTIMIZE = -g -w -m64 -O3 -march=native
ifeq (NUM_THREADS,$(findstring NUM_THREADS,$(CONFIGVARS)))
OPTIMIZE += -fopenmp
else
OPTIMIZE += -Wno-unknown-pragmas
endif
GSL_INCL =
GSL_LIBS =
FFTW_INCL= -I/opt/local/include -I/usr/local/include
FFTW_LIBS= -L/opt/local/lib -I/usr/local/lib
GMP_INCL =
GMP_LIBS =
MPICHLIB =
HDF5INCL =
HDF5LIB =
#OPT += -DNOCALLSOFSYSTEM
#OPT += -DIMPOSE_PINNING
#OPT += -DUSE_SSE
endif
ifndef LINKER
LINKER = $(CC)
endif
##########################################
#determine the needed object/header files#
##########################################
SUBDIRS = .
OBJS = main.o allocate.o allvars.o disk.o grid.o bulge.o set_particles.o parallel_sort.o \
halo.o init.o io.o mymalloc.o orbit_response.o parameters.o structure.o system.o disp_fields.o \
forcetree/gravtree.o forcetree/forcetree.o forcetree/forcetree_walk.o domain/peano.o domain/pqueue.o \
domain/domain.o domain/domain_balance.o domain/domain_counttogo.o domain/domain_exchange.o \
domain/domain_rearrange.o domain/domain_sort_kernels.o domain/domain_toplevel.o domain/domain_vars.o domain/domain_box.o
INCL += allvars.h proto.h
SUBDIRS += forcetree domain
################################
#determine the needed libraries#
################################
ifneq (HAVE_HDF5,$(findstring HAVE_HDF5,$(CONFIGVARS)))
HDF5LIB =
endif
ifeq (NUM_THREADS,$(findstring NUM_THREADS,$(CONFIGVARS)))
THREAD_LIB =
endif
##########################
#combine compiler options#
##########################
CFLAGS = $(OPTIMIZE) $(OPT) $(HDF5INCL) $(GSL_INCL) $(FFTW_INCL) $(ODE_INCL) $(GMP_INCL) $(MKL_INCL) $(CUDA_INCL) -I$(BUILD_DIR)
LIBS = $(MATHLIB) $(HDF5LIB) $(MPICHLIB) $(GSL_LIBS) $(GSLLIB) $(FFTW_LIB) $(GMP_LIBS) $(GMPLIB) $(ODE_LIB) $(MKL_LIBS) $(THREAD_LIB) $(CUDA_LIBS)
SUBDIRS := $(addprefix $(BUILD_DIR)/,$(SUBDIRS))
OBJS := $(addprefix $(BUILD_DIR)/,$(OBJS)) $(BUILD_DIR)/compile_time_info.o
INCL := $(addprefix $(SRC_DIR)/,$(INCL)) $(BUILD_DIR)/galicconfig.h
################
#create subdirs#
################
RESULT := $(shell mkdir -p $(SUBDIRS) )
#############
#build rules#
#############
all: $(EXEC)
$(EXEC): $(OBJS)
$(LINKER) $(OPTIMIZE) $(OBJS) $(LIBS) -o $(EXEC)
# mpirun -n $(N) -f hostfile ./$(EXEC) $(PARAMFILE)
#bg: $(OBJS)
# $(LINKER) $(OPTIMIZE) $(OBJS) $(LIBS) -o $(EXEC)
# mpirun -n $(N) -f hostfile ./$(EXEC) $(PARAMFILE) 1> log.out.txt 2> log.err.txt &
clean:
rm -f $(OBJS) $(EXEC) lib$(LIBRARY).a
rm -f $(BUILD_DIR)/compile_time_info.c $(BUILD_DIR)/galicconfig.h
$(BUILD_DIR)/%.o: $(SRC_DIR)/%.c $(INCL) $(MAKEFILES)
$(CC) $(CFLAGS) -c $< -o $@
$(BUILD_DIR)/compile_time_info.o: $(BUILD_DIR)/compile_time_info.c $(MAKEFILES)
$(CC) $(CFLAGS) -c $< -o $@
""" @package ./examples/isolated_galaxy_collisionless_3d/check.py
Code that checks results of 3d isolated collisionless galaxy problem
created by Rainer Weinberger, last modified 10.03.2019
"""
""" load libraries """
import sys ## system calls
import numpy as np ## load numpy
import h5py ## load h5py; needed to read snapshots
import matplotlib.pyplot as plt ## plot stuff
from matplotlib.colors import LogNorm
createFigures = False
simulation_directory = str(sys.argv[1])
print("examples/isolated_galaxy_collisionless_3d/check.py: checking simulation output in directory " + simulation_directory)
FloatType = np.float64 # double precision: np.float64, for single use np.float32
Gcgs = 6.6738e-8
status=0
rGrid = np.linspace(1,34,100)
mEnc1Grid_ref = np.zeros(rGrid.shape)
mEnc2Grid_ref = np.zeros(rGrid.shape)
mEnc3Grid_ref = np.zeros(rGrid.shape)
""" loop over all output files """
i_file = 0
while True:
""" try to read in snapshot """
directory = simulation_directory+"/output/"
filename = "snap_%03d.hdf5" % (i_file)
try:
data = h5py.File(directory+filename, "r")
except:
break
""" get simulation data """
time = FloatType(data["Header"].attrs["Time"])
Pos1 = np.array(data["PartType1"]["Coordinates"], dtype = FloatType)
Mass1 = np.array(data["PartType1"]["Masses"], dtype=FloatType)
Pos2 = np.array(data["PartType2"]["Coordinates"], dtype = FloatType)
Mass2 = np.array(data["PartType2"]["Masses"], dtype=FloatType)
Pos3 = np.array(data["PartType3"]["Coordinates"], dtype = FloatType)
Mass3 = np.array(data["PartType3"]["Masses"], dtype=FloatType)
Center = np.zeros(3) ## center defined as center of disk
Center[0] = np.average(Pos2[:,0])
Center[1] = np.average(Pos2[:,1])
Center[2] = np.average(Pos2[:,2])
print('center: %g %g %g'%(Center[0], Center[1], Center[2]))
xPosFromCenter1 = (Pos1[:,0] - Center[0])
yPosFromCenter1 = (Pos1[:,1] - Center[1])
zPosFromCenter1 = (Pos1[:,2] - Center[2])
xPosFromCenter2 = (Pos2[:,0] - Center[0])
yPosFromCenter2 = (Pos2[:,1] - Center[1])
zPosFromCenter2 = (Pos2[:,2] - Center[2])
xPosFromCenter3 = (Pos3[:,0] - Center[0])
yPosFromCenter3 = (Pos3[:,1] - Center[1])
zPosFromCenter3 = (Pos3[:,2] - Center[2])
Radius1 = np.sqrt( xPosFromCenter1**2 + yPosFromCenter1**2 + zPosFromCenter1**2 )
Radius2 = np.sqrt( xPosFromCenter2**2 + yPosFromCenter2**2 + zPosFromCenter2**2 )
Radius3 = np.sqrt( xPosFromCenter3**2 + yPosFromCenter3**2 + zPosFromCenter3**2 )
""" disk height """
h_disp = np.std(zPosFromCenter2)
if h_disp > 1.2:
status = 1
print("ERROR: %d disk height: %g"% (i_file, h_disp) )
""" rotation curve / enclosed mass """
## sort by radius
i_sort1 = np.argsort(Radius1)
i_sort2 = np.argsort(Radius2)
i_sort3 = np.argsort(Radius3)
## calculate enclosed mass
mEnc1 = np.cumsum(Mass1[i_sort1])
mEnc2 = np.cumsum(Mass2[i_sort2])
mEnc3 = np.cumsum(Mass3[i_sort3])
#interpolate on uniform grid
mEnc1Grid = np.interp( rGrid, Radius1[i_sort1], mEnc1 )
mEnc2Grid = np.interp( rGrid, Radius2[i_sort2], mEnc2 )
mEnc3Grid = np.interp( rGrid, Radius3[i_sort3], mEnc3 )
mEncTotGrid = mEnc1Grid + mEnc2Grid + mEnc3Grid
## calculate rotation curve
vCirc1Grid = np.sqrt(mEnc1Grid * Gcgs / rGrid * 6.4459e11)
vCirc2Grid = np.sqrt(mEnc2Grid * Gcgs / rGrid * 6.4459e11)
vCirc3Grid = np.sqrt(mEnc3Grid * Gcgs / rGrid * 6.4459e11)
vTotGrid = np.sqrt(mEncTotGrid * Gcgs / rGrid * 6.4459e11)
if createFigures:
## figure of positions
fig = plt.figure(figsize = np.array([6.9,9.2]))
ax = []
ax.append(plt.axes([0.1,0.41,0.75,0.56]) )
ax.append(plt.axes([0.1,0.07,0.75,0.282]) )
cax = plt.axes([0.86,0.07,0.04,0.90])
xext = 15.0
binsx = np.linspace(-xext,xext,128)
binsy = np.linspace(-xext,xext,128)
binsz = np.linspace(-0.5*xext,0.5*xext,64)
H1, xgrid, ygrid = np.histogram2d(xPosFromCenter2, yPosFromCenter2, bins=[binsx, binsy], density=True)
H2, xgrid, ygrid = np.histogram2d(xPosFromCenter3, yPosFromCenter3, bins=[binsx, binsy], density=True)
He1, xgrid, ygrid = np.histogram2d(zPosFromCenter2, yPosFromCenter2, bins=[binsz, binsy], density=True)
He2, xgrid, ygrid = np.histogram2d(zPosFromCenter3, yPosFromCenter3, bins=[binsz, binsy], density=True)
ext=(-xext,xext,-xext,xext)
ax[0].imshow(H1+H2, \
norm=LogNorm(vmin=3e-4,vmax=3e-1), \
cmap='Greys', \
origin='lower', \
extent=ext \
)
ext=(-xext,xext,-0.5*xext,0.5*xext)
img = ax[1].imshow(He1+He2, \
norm=LogNorm(vmin=3e-4,vmax=3e-1), \
cmap='Greys', \
origin='lower', \
extent=ext \
)
ax[0].axis('equal')
ax[1].axis('equal')
ax[0].set_xticks([-10,-5,0,5,10])
ax[0].set_yticks([-10,-5,0,5,10])
ax[1].set_xticks([-10,-5,0,5,10])
ax[1].set_yticks([-5,0,5])
ax[0].set_xlabel('[kpc]')
ax[0].set_ylabel('[kpc]')
ax[1].set_xlabel('[kpc]')
ax[1].set_ylabel('[kpc]')
fig.colorbar(img, cax=cax)
fig.savefig(directory+"/positions_%03d.pdf"%i_file)
## figure of rotation curve
fig = plt.figure()
ax = plt.axes([0.1, 0.1, 0.85, 0.85])
ax.plot( rGrid, vCirc1Grid )
ax.plot( rGrid, vCirc2Grid )
ax.plot( rGrid, vCirc3Grid )
ax.plot( rGrid, vTotGrid )
ax.set_xlim(0,35)
ax.set_ylim(0,300)
ax.set_xlabel(r"radius")
ax.set_ylabel(r"v$_c$")
fig.savefig(directory+"/rotation_curve_%03d.pdf"%i_file)
## comparision to first snapshot (ICs)
if i_file == 0:
mEnc1Grid_ref[:] = mEnc1Grid[:]
mEnc2Grid_ref[:] = mEnc2Grid[:]
mEnc3Grid_ref[:] = mEnc3Grid[:]
else:
tolerance = 1.0
delta_mEnc = mEnc1Grid_ref - mEnc1Grid
if np.max( abs(delta_mEnc) ) > tolerance:
status = 1
print("ERROR: i_file: %d: enclused mass difference in halo (part1) exceeds tolerance %g" % (i_file, tolerance) )
print(delta_mEnc)
delta_mEnc = mEnc2Grid_ref - mEnc2Grid
if np.max( abs(delta_mEnc) ) > tolerance:
status = 1
print("ERROR: i_file: %d: enclused mass difference in disc (part2) exceeds tolerance %g" % (i_file, tolerance) )
print(delta_mEnc)
delta_mEnc = mEnc3Grid_ref - mEnc3Grid
if np.max( abs(delta_mEnc) ) > tolerance:
status = 1
print("ERROR: i_file: %d: enclused mass difference in bulge (part3) exceeds tolerance %g" % (i_file, tolerance) )
print(delta_mEnc)
i_file += 1
## if everything is ok: 0; else: 1
sys.exit(status)
""" @package examples/collisionless_galaxy_3d/create.py
Code creates the output list; ICs need to be present already
collisionless_galaxy_3d uses GalIC initial condition code, Model M4
created by Rainer Weinberger, last modified 29.04.2018
"""
""" load libraries """
import sys # system calls
import numpy as np # scientific computing package
import h5py # hdf5 format
import os # operating system interface
from subprocess import call # execute bash commands
createICs=False
""" input """
simulation_directory = str(sys.argv[1])
print("examples/collisionless_galaxy_3d/create.py " + simulation_directory)
""" creating ics """
if createICs:
status = call(['git', 'clone', 'https://github.com/djurin/GALIC.git', simulation_directory+'/galic'])
if status != 0:
print('CREATE: ERROR: git clone failed!')
sys.exit(status)
cwd = os.getcwd()
os.chdir(simulation_directory+'/galic/')
## copy makefile and config-file
str = 'SYSTYPE="Darwin"'
file=open('Makefile.systype','w')
file.write(str+'\n')
file.close()
call(['cp', cwd+'/examples/isolated_galaxy_collisionless_3d/Makefile.galic', 'Makefile'])
str = 'DOUBLEPRECISION=1 \nHAVE_HDF5 \nVER_1_1 \n'
file=open('Config.sh','w')
file.write(str+'\n')
file.close()
status = call(['make'])
if status != 0:
print('CREATE: ERROR: make failed!')
sys.exit(status)
status = call(['mpiexec','-np','1','./galic',cwd+'/examples/isolated_galaxy_collisionless_3d/param_galic.txt'])
if status != 0:
print('CREATE: ERROR: execution failed!')
sys.exit(status)
os.chdir(cwd)
else:
call(['cp', './examples/isolated_galaxy_collisionless_3d/ICs/ICs', simulation_directory+'/snap_010'])
""" set output times """
outputTimes = np.linspace(0.0,1.0,10, dtype=np.float64)
ones = np.ones(outputTimes.shape, dtype=np.int)
""" write output list file """
data = np.array([outputTimes, ones]).T
np.savetxt(simulation_directory+"/output_list.txt",data, fmt="%g %1.f" )
""" normal exit """
sys.exit(0)
%---- Relevant files
InitCondFile run/examples/isolated_galaxy_collisionless_3d/snap_010
OutputDir run/examples/isolated_galaxy_collisionless_3d/output
SnapshotFileBase snap
OutputListFilename run/examples/isolated_galaxy_collisionless_3d/output_list.txt
%---- File formats
ICFormat 1
SnapFormat 3
%---- CPU-time limits
TimeLimitCPU 93000
CpuTimeBetRestartFile 12000
ResubmitOn 0
ResubmitCommand my-scriptfile
%----- Memory alloction
MaxMemSize 2500
%---- Caracteristics of run
TimeBegin 0.0
TimeMax 1.0 % End of the simulation
%---- Basic code options that set the type of simulation
ComovingIntegrationOn 0
PeriodicBoundariesOn 0
CoolingOn 0
StarformationOn 0
%---- Cosmological parameters (Planck cosmology)
Omega0 0.0
OmegaLambda 0.0
OmegaBaryon 0.0
HubbleParam 1.0
BoxSize 100000.0
%---- Output frequency and output parameters
OutputListOn 1
TimeBetSnapshot 0.0
TimeOfFirstSnapshot 0.0
TimeBetStatistics 0.01
NumFilesPerSnapshot 1
NumFilesWrittenInParallel 1
%---- Accuracy of time integration
TypeOfTimestepCriterion 0
ErrTolIntAccuracy 0.012
CourantFac 0.3
MaxSizeTimestep 0.05
MinSizeTimestep 2.0e-9
%---- Treatment of empty space and temperature limits
InitGasTemp 244.8095
MinGasTemp 5.0
MinimumDensityOnStartUp 1.0e-20
LimitUBelowThisDensity 0.0
LimitUBelowCertainDensityToThisValue 0.0
MinEgySpec 0.0
%---- Tree algorithm, force accuracy, domain update frequency
TypeOfOpeningCriterion 1
ErrTolTheta 0.7
ErrTolForceAcc 0.0025
MultipleDomains 8
TopNodeFactor 2.5
ActivePartFracForNewDomainDecomp 0.01
%---- Initial density estimate
DesNumNgb 64
MaxNumNgbDeviation 4
%---- System of units
UnitLength_in_cm 3.085678e21 % 1.0 kpc
UnitMass_in_g 1.989e43 % 1.0e10 solar masses
UnitVelocity_in_cm_per_s 1e5 % 1 km/sec
GravityConstantInternal 0
%---- Gravitational softening lengths
SofteningComovingType0 1.0
SofteningComovingType1 1.0
SofteningMaxPhysType0 1.0
SofteningMaxPhysType1 1.0
GasSoftFactor 2.5
SofteningTypeOfPartType0 0
SofteningTypeOfPartType1 1
SofteningTypeOfPartType2 1
SofteningTypeOfPartType3 1
SofteningTypeOfPartType4 1
SofteningTypeOfPartType5 1
MinimumComovingHydroSoftening 1.0
AdaptiveHydroSofteningSpacing 1.2
%----- Mesh regularization options
CellShapingSpeed 0.5
CellShapingFactor 1.0
%------ File and path names, as well as output file format
OutputDir ../
OutputFile snap % Base filename of generated sequence of files
SnapFormat 1 % File format selection
%------ Basic structural parameters of model
CC 10.0 % halo concentration
V200 200.0 % circular velocity v_200 (in km/sec)
LAMBDA 0.035 % spin parameter
MD 0.035 % disk mass fraction
MB 0.05 % bulge mass fraction
MBH 0.0 % black hole mass fraction. If zero, no black
% hole is generated, otherwise one at the centre
% is added.
JD 0.035 % disk spin fraction, typically chosen equal to MD
DiskHeight 0.2 % thickness of stellar disk in units of radial scale length
BulgeSize 0.1 % bulge scale length in units of halo scale length
HaloStretch 1.0 % should be one for a spherical halo, smaller than one corresponds to prolate distortion, otherwise oblate
BulgeStretch 1.0 % should be one for a spherical bulge, smaller than one corresponds to prolate distortion, otherwise oblate
%------ Particle numbers in target model
N_HALO 100000 % desired number of particles in dark halo