added example cosmo_box_gravity_only_3d

examples/cosmo_box_gravity_only_3d/Config.sh

+#!/bin/bash        	# this line only there to enable syntax highlighting in this file
+
+##################################################
+#  Enable/Disable compile-time options as needed #
+#  examples/cosmo_box_gravity_only_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 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
+
+#--------------------------------------- TreePM Options
+PMGRID=256                               # Enables particle mesh; number of cells used for grid in each dimension
+RCUT=5.0                                 # This can be used to override the maximum radius in which the short-range tree-force is evaluated (in case the TreePM algorithm is used). The default value is 4.5, given in mesh-cells.
+
+#--------------------------------------- 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.
+DOUBLEPRECISION_FFTW                     # FFTW calculation in double precision
+OUTPUT_COORDINATES_IN_DOUBLEPRECISION    # will always output coordinates in double precision
+NGB_TREE_DOUBLEPRECISION                 # if this is enabled, double precision is used for the neighbor node extension
+
+#--------------------------------------- On the fly FOF groupfinder
+FOF                                      # enable FoF output
+FOF_PRIMARY_LINK_TYPES=2                 # 2^type for the primary dark matter type
+FOF_SECONDARY_LINK_TYPES=1+16+32         # 2^type for the types linked to nearest primaries
+
+#--------------------------------------- Subfind
+SUBFIND                                  # enables substructure finder
+SAVE_HSML_IN_SNAPSHOT                    # stores hsml, density, and velocity dispersion values in the snapshot files
+SUBFIND_CALC_MORE                        # calculates also the velocity dispersion in the local density estimate (this is automatically enabled by several other options, e.g. SAVE_HSML_IN_SNAPSHOT)
+SUBFIND_EXTENDED_PROPERTIES              # adds calculation of further quantities related to angular momentum in different components
+
+#-------------------------------------------- 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
+REDUCE_FLUSH                             # only flush output to log-files in predefined intervals
+HAVE_HDF5                                # needed when HDF5 I/O support is desired (recommended)
+
+#--------------------------------------- Testing and Debugging options
+DEBUG                                    # enables core-dumps

examples/cosmo_box_gravity_only_3d/Masses_L50n32_z0.txt

+4.222264453125000000e+14
+1.860935937500000000e+14
+9.680791992187500000e+13
+9.419149414062500000e+13
+8.993979492187500000e+13
+8.568809570312500000e+13
+8.372577148437500000e+13
+7.293299316406250000e+13
+6.671896972656250000e+13
+6.606486328125000000e+13
+6.148611328125000000e+13
+6.115905761718750000e+13
+4.578752929687500000e+13
+4.088172363281250000e+13
+3.957350830078125000e+13
+3.761118408203125000e+13
+3.630297119140625000e+13
+3.139716308593750000e+13
+2.878073242187500000e+13
+2.420197998046875000e+13
+2.387492675781250000e+13
+2.354787353515625000e+13
+2.256671142578125000e+13
+2.223965820312500000e+13
+2.191260253906250000e+13
+2.093144287109375000e+13
+2.093144287109375000e+13
+2.060438720703125000e+13
+2.027733520507812500e+13
+1.962322753906250000e+13
+1.962322753906250000e+13
+1.896911987304687500e+13
+1.831501220703125000e+13
+1.831501220703125000e+13
+1.733385009765625000e+13
+1.733385009765625000e+13
+1.733385009765625000e+13
+1.700679687500000000e+13
+1.667974243164062500e+13
+1.602563598632812500e+13
+1.602563598632812500e+13
+1.569858154296875000e+13
+1.537152832031250000e+13
+1.504447387695312500e+13
+1.504447387695312500e+13
+1.373625854492187500e+13
+1.275509765625000000e+13
+1.275509765625000000e+13
+1.242804321289062500e+13
+1.242804321289062500e+13
+1.210098999023437500e+13
+1.210098999023437500e+13
+1.144688232421875000e+13
+1.144688232421875000e+13
+1.144688232421875000e+13
+1.111982910156250000e+13
+1.111982910156250000e+13
+1.111982910156250000e+13
+1.046572143554687500e+13
+1.046572143554687500e+13
+1.046572143554687500e+13
+1.046572143554687500e+13

examples/cosmo_box_gravity_only_3d/Masses_L50n32_z1.txt

+1.040031054687500000e+14
+6.868129394531250000e+13
+6.606486328125000000e+13
+5.429092773437500000e+13
+3.663002441406250000e+13
+3.466770019531250000e+13
+2.976189453125000000e+13
+2.910778564453125000e+13
+2.845367919921875000e+13
+2.191260253906250000e+13
+2.060438720703125000e+13
+1.929617309570312500e+13
+1.896911987304687500e+13
+1.700679687500000000e+13
+1.635268920898437500e+13
+1.504447387695312500e+13
+1.504447387695312500e+13
+1.504447387695312500e+13
+1.504447387695312500e+13
+1.308215087890625000e+13
+1.275509765625000000e+13
+1.242804321289062500e+13
+1.177393676757812500e+13
+1.111982910156250000e+13
+1.079277465820312500e+13

examples/cosmo_box_gravity_only_3d/check.py

+""" @package ./examples/cosmo_box_gravity_only_3d/check.py
+Code that checks results of gravity only structure formation simulation
+
+created by Rainer Weinberger, last modified 25.02.2019
+"""
+
+""" load libraries """
+import sys    ## load sys; needed for exit codes
+import numpy as np    ## load numpy
+import h5py    ## load h5py; needed to read snapshots
+import matplotlib.pyplot as plt    ## plot stuff
+
+simulation_directory = str(sys.argv[1])
+print("examples/cosmo_box_gravity_only_3d/check.py: checking simulation output in directory " + simulation_directory) 
+
+FloatType = np.float64  # double precision: np.float64, for single use np.float32
+
+Boxsize = 50 ##Mpc/h
+HubbleParam = 0.6774
+UnitMass = 1.0e10
+Volume = Boxsize * Boxsize * Boxsize / HubbleParam / HubbleParam / HubbleParam
+
+Redshifts = [1, 0]
+status = 0
+
+CompareAgainstReferenceRun = True ## comparison for small L50m32 box; deactivate this when comparing against self-created ICs
+createFigures = True
+
+for i_file, z in enumerate(Redshifts):
+    """ try to read in snapshot """
+    directory = simulation_directory+"/output/"
+    filename = "fof_subhalo_tab_%03d.hdf5" % (i_file)
+    try:
+        data = h5py.File(directory+filename, "r")
+    except:
+        print("could not open "+directory+filename)
+        sys.exit(1)
+    
+    """ get simulation data """
+    ## simulation data
+    GrpPos = np.array(data["Group"]["GroupPos"], dtype=FloatType) / HubbleParam / 1000.
+    GrpR200c = np.array(data["Group"]["Group_R_Crit200"], dtype=FloatType) / HubbleParam / 1000.
+    M200c = np.array(data["Group"]["GroupMass"], dtype=FloatType) * UnitMass
+    M200c = np.sort(M200c)[::-1]
+    CumMassFunction = np.cumsum(np.ones(M200c.shape) ) / Volume
+    
+    if CompareAgainstReferenceRun:
+        ## comparison to reference run (sorted list of M200)
+        M200c_ref = np.loadtxt("./examples/cosmo_box_gravity_only_3d/Masses_L50n32_z%.1d.txt"% z)
+        
+        minLen = np.min([len(M200c), len(M200c_ref)])
+        i_select = np.arange(minLen)
+        
+        delta = (M200c[i_select]-M200c_ref[i_select]) / M200c_ref[i_select]
+        
+        ## empirically based tolerances
+        tolerance_average = 0.01
+        tolerance_std = 0.05
+        if np.abs(np.average(delta)) > tolerance_average or np.abs(np.std(delta)) > tolerance_std:
+            status = 1
+            print("ERROR: z=%g difference in halo masses exceeding limits!" % z)
+            print("relative mass error (=delta)")
+            print(delta)
+            print("average delta (tolerance: %g)" % tolerance_average)
+            print(np.average(delta))
+            print("stddev delta (tolerance: %g)" % tolerance_std)
+            print(np.std(delta))
+    
+    """ optional figure """
+    if createFigures:
+        filename = "snap_%03d.hdf5" % (i_file)
+        try:
+            data = h5py.File(directory+filename, "r")
+        except:
+            print("could not open "+directory+filename)
+            sys.exit(1)
+        pos = np.array(data["PartType1"]["Coordinates"], dtype=FloatType) / HubbleParam / 1000.
+    
+        fig = plt.figure(figsize=(6.9,6.9))
+        ax = plt.axes([0.1,0.1,0.87,0.87])
+        
+        if(pos.shape[0] > 32**3):
+            i_select = np.random.uniform(low=0.0, high=pos.shape[0], size=32**3).astype(np.int)
+        else:
+            i_select = np.arange(pos.shape[0])
+        ax.scatter(pos[i_select, 0], pos[i_select, 1], marker='.', s=0.05, alpha=0.5, rasterized=True)
+        for i in np.arange(GrpR200c.shape[0]):
+            ax.add_artist(plt.Circle((GrpPos[i,0], GrpPos[i,1]), GrpR200c[i], color='k', fill=False))
+        
+        ax.set_xlim([0,Boxsize/HubbleParam])
+        ax.set_ylim([0,Boxsize/HubbleParam])
+        ax.set_xlabel('[Mpc]')
+        ax.set_ylabel('[Mpc]')
+        
+        bx = plt.axes([0.70,0.74,0.26,0.22])
+        bx.plot(M200c, CumMassFunction)
+        bx.set_xscale('log')
+        bx.set_yscale('log')
+        bx.set_xlim([9e12,5e14])
+        bx.set_ylim([9e-7,2e-4])
+        bx.set_xlabel(r"M$_{200,c}$ [M$_\odot$]")
+        bx.set_ylabel(r"$n$(>M) [Mpc$^{-3}$]")
+        fig.savefig(simulation_directory+'/largeScaleStructure_z%.1d.pdf'%z, dpi=300)
+
+## if everything is ok: 0 else: 1
+sys.exit(status)

examples/cosmo_box_gravity_only_3d/create.py

+""" @package examples/cosmo_box_gravity_only_3d/create.py
+Code that creates 3d cosmologiacl box ICs;
+some ICs provided, but can be created also at startup if 
+ic_creation is set to 'music' or 'ngenic'.
+
+created by Rainer Weinberger, last modified 25.02.2019
+"""
+
+
+""" load libraries """
+import sys    # needed for exit codes
+import numpy as np    # scientific computing package
+import h5py    # hdf5 format
+import os  # operating system interface
+from subprocess import call # execute bash commands
+
+## create new ics with 'music' or 'ngenic' or just 'copy' existing ones to the run directory
+ic_creation='copy'
+
+
+""" input """
+simulation_directory = str(sys.argv[1])
+print("examples/cosmo_box_gravity_only_3d/create.py " + simulation_directory)
+
+
+""" initial conditions: either copy or create with code """
+if ic_creation == 'copy':
+    ## copy and use provided initial conditions
+    call(['cp', './examples/cosmo_box_gravity_only_3d/L50n32/ics', simulation_directory+'/ics'])
+    
+elif ic_creation == 'music':
+    ## create new initial conditions with the MUSIC code
+    status = call(['hg', 'clone', 'https://bitbucket.org/ohahn/music', simulation_directory+'/music'])
+    if status != 0:
+        print('CREATE: ERROR: hg clone failed!')
+        sys.exit(status)
+    cwd = os.getcwd()
+    os.chdir(simulation_directory+'/music/')
+    status = call(['make'])
+    if status != 0:
+        print('CREATE: ERROR: make failed!')
+        sys.exit(status)
+    status = call(['./MUSIC',cwd+'/examples/cosmo_box_gravity_only_3d/param_music.txt'])
+    if status != 0:
+        print('CREATE: ERROR: execution failed!')
+        sys.exit(status)
+    os.chdir(cwd)
+    
+elif ic_creation == 'ngenic':
+    ## create new initial conditions with the N-GenIC code
+    status = call(['git', 'clone', 'https://gitlab.mpcdf.mpg.de/ext-c2c74fbfcdff/ngenic.git', simulation_directory+'/ngenic'])
+    if status != 0:
+        print('CREATE: ERROR: git clone failed!')
+        sys.exit(status)
+    cwd = os.getcwd()
+    os.chdir(simulation_directory+'/ngenic/')
+    status = call(['make'])
+    if status != 0:
+        print('CREATE: ERROR: make failed!')
+        sys.exit(status)
+    status = call(['mpiexec','-np','1','./N-GenIC',cwd+'/examples/cosmo_box_gravity_only_3d/param_ngenic.txt'])
+    if status != 0:
+        print('CREATE: ERROR: execution failed!')
+        sys.exit(status)
+    os.chdir(cwd)
+else:
+    print("CREATE: ERROR: no valid option for ic creation! choose 'copy', 'music' or 'ngenic'")
+    exit(1)
+
+
+""" set output times z=1,0"""
+outputTimes = np.array([0.5, 1], 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) 

examples/cosmo_box_gravity_only_3d/param.txt

+%----  Relevant files
+InitCondFile                          run/examples/cosmo_box_gravity_only_3d/ics
+OutputDir                             run/examples/cosmo_box_gravity_only_3d/output
+SnapshotFileBase                      snap
+OutputListFilename                    run/examples/cosmo_box_gravity_only_3d/output_list.txt
+
+%---- File formats
+ICFormat                              1
+SnapFormat                            3
+
+%---- CPU-time limits
+TimeLimitCPU                          93000
+CpuTimeBetRestartFile                 12000
+FlushCpuTimeDiff                      120
+
+ResubmitOn                            0
+ResubmitCommand                       my-scriptfile
+
+%----- Memory alloction
+MaxMemSize                            2500
+
+%---- Caracteristics of run
+TimeBegin                             0.0078125 % Begin of the simulation z=127
+TimeMax                               1.0 % End of the simulation z=0
+
+%---- Basic code options that set the type of simulation
+ComovingIntegrationOn                 1
+PeriodicBoundariesOn                  1
+CoolingOn                             0
+StarformationOn                       0
+
+%---- Cosmological parameters (Planck 2016 cosmology)
+Omega0                                0.3089
+OmegaLambda                           0.6911
+OmegaBaryon                           0.0     %0.0486
+HubbleParam                           0.6774
+BoxSize                               50000.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.005
+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                       4
+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                40.0
+SofteningComovingType1                40.0
+
+SofteningMaxPhysType0                 20.0
+SofteningMaxPhysType1                 20.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
+
+%----- Subfind
+ErrTolThetaSubfind                    0.7
+DesLinkNgb                            20

examples/cosmo_box_gravity_only_3d/param_music.txt

+[setup]
+boxlength          = 50
+zstart             = 127
+levelmin           = 5
+levelmin_TF        = 5
+levelmax           = 5
+padding            = 8
+overlap            = 4
+ref_center         = 0.5, 0.5, 0.5
+ref_extent         = 0.4, 0.4, 0.4
+align_top          = no
+baryons            = no
+use_2LPT           = no
+use_LLA            = no
+periodic_TF        = yes
+[cosmology]
+Omega_m            = 0.3089
+Omega_L            = 0.6911
+w0                 = -1.0
+wa                 = 0.0
+Omega_b            = 0.0486
+H0                 = 67.74
+sigma_8            = 0.831
+nspec              = 0.9645
+transfer           = eisenstein
+[random]
+seed[7]            = 12345
+seed[8]            = 23456
+seed[9]            = 34567
+seed[10]           = 45678
+seed[11]           = 56789
+seed[12]           = 67890
+[output]
+#Gadget-2 (type=1: high-res particles, type=5: rest)
+format             = gadget2
+filename           = ../ics
+gadget_coarsetype  = 2
+gadget_usekpc      = yes
+[poisson]
+fft_fine           = yes
+accuracy           = 1e-5
+pre_smooth         = 3
+post_smooth        = 3
+smoother           = gs
+laplace_order      = 6
+grad_order         = 6
\ No newline at end of file

examples/cosmo_box_gravity_only_3d/param_ngenic.txt

+
+Nmesh            128                        % This is the size of the FFT grid used to 
+                                            % compute the displacement field. One
+                                            % should have Nmesh >= Nsample.
+
+Nsample          128                        % sets the maximum k that the code uses,
+                                            % i.e. this effectively determines the
+                                            % Nyquist frequency that the code assumes,
+                                            % k_Nyquist = 2*PI/Box * Nsample/2
+                                            % Normally, one chooses Nsample such that
+                                            % Ntot =  Nsample^3, where Ntot is the
+                                            % total number of particles
+ 
+
+Box              50000.0                    % Periodic box size of simulation
+FileBase         ics                        % Base-filename of output files
+OutputDir        ../                        % Directory for output
+GlassFile        dummy_glass.dat            % File with unperturbed glass or
+                                            % Cartesian grid
+TileFac          4                          % Number of times the glass file is
+                                            % tiled in each dimension (must be
+                                            % an integer)
+
+
+Omega            0.3089                     % Total matter density  (at z=0)
+OmegaLambda      0.6911                     % Cosmological constant (at z=0)
+OmegaBaryon      0.0                        % Baryon density        (at z=0)
+HubbleParam      0.6774                     % Hubble paramater (may be used for power spec parameterization)
+Redshift         127                        % Starting redshift
+Sigma8           0.831                      % power spectrum normalization
+
+
+SphereMode       1                          % if "1" only modes with |k| < k_Nyquist are
+                                            % used (i.e. a sphere in k-space), otherwise modes with
+                                            % |k_x|,|k_y|,|k_z| < k_Nyquist are used
+                                            % (i.e. a cube in k-space)
+
+
+WhichSpectrum    1                          % "1" selects Eisenstein & Hu spectrum,
+                                            % "2" selects a tabulated power spectrum in
+                                            % the file 'FileWithInputSpectrum'
+                                            % otherwise, Efstathiou parametrization is used
+
+
+FileWithInputSpectrum   input_spectrum.txt  % filename of tabulated input
+                                            % spectrum (if used)
+InputSpectrum_UnitLength_in_cm  3.085678e24 % defines length unit of tabulated
+                                            % input spectrum in cm/h. 
+                                            % Note: This can be chosen different from UnitLength_in_cm
+
+ReNormalizeInputSpectrum   1                % if set to zero, the
+                                            % tabulated spectrum is
+                                            % assumed to be normalized
+                                            % already in its amplitude to
+                                            % the starting redshift,
+                                            % otherwise this is recomputed
+                                            % based on the specified sigma8
+
+
+ShapeGamma       0.21                       % only needed for Efstathiou power spectrum 
+PrimordialIndex  0.9645                     % may be used to tilt the primordial index, 
+                                            % primordial spectrum is k^PrimordialIndex
+
+
+Seed             123456                     %  seed for IC-generator
+
+
+NumFilesWrittenInParallel 1                 % limits the number of files that are
+                                            % written in parallel when outputting
+
+
+UnitLength_in_cm          3.085678e21       % defines length unit of output (in cm/h) 
+UnitMass_in_g             1.989e43          % defines mass unit of output (in g/cm)
+UnitVelocity_in_cm_per_s  1e5               % defines velocity unit of output (in cm/sec)
+
+
+
+
+
+
+
+

test.sh

@@ -23,6 +23,7 @@ TESTS+="yee_2d "
 
 ## available 3d examples
 TESTS+="noh_3d "
+TESTS+="cosmo_box_gravity_only_3d "
 
 ## loop over all tests
 for TEST in $TESTS