Cosmology System Test & Gravity Support in SystemTestRunner Class

builds/make.type.cosmology

@@ -20,6 +20,8 @@ DFLAGS += -DCOSMOLOGY
 # Average Slow cell when the cell delta_t is very small
 DFLAGS += -DAVERAGE_SLOW_CELLS
 
+# Solve the Gas Internal Energy using a Dual Energy Formalism
+DFLAGS    += -DDE
 
 # Print Initial Statistics
 DFLAGS += -DPRINT_INITIAL_STATS

builds/make.type.hydro

@@ -16,7 +16,7 @@ endif
 DFLAGS    += -DDENSITY_FLOOR
 DFLAGS    += -DTEMPERATURE_FLOOR
 
-# Solve the Gas Internal Energy usisng a Dual Energy Formalism
+# Solve the Gas Internal Energy using a Dual Energy Formalism
 #DFLAGS    += -DDE
 
 # Apply cooling on the GPU from precomputed tables

python_scripts/combine_test_files.py

@@ -0,0 +1,156 @@
+#!/usr/bin/env python3
+"""Python script for merging hydro, particles, and gravity HDF5 files when they're used as fiducial data files for testing
+
+Raises
+------
+ValueError
+    Duplicate datasets in destination and particle data files
+ValueError
+    Duplicate datasets in destination and gravity data files
+"""
+
+import h5py
+import pathlib
+import argparse
+import shutil
+
+# =====================================================================================================================
+def main():
+    # Initialize the CLI
+    cli = argparse.ArgumentParser()
+
+    # Required Arguments
+    cli.add_argument('-s', '--source-directory', type=pathlib.Path,  required=True, help='The path to the directory for the source HDF5 files.')
+    cli.add_argument('-o', '--output-directory', type=pathlib.Path,  required=True, help='The path to the directory to write out the concatenated HDF5 file.')
+
+    # Optional Arguments
+    cli.add_argument('-p', '--particles', default=False, action='store_true', help='')
+    cli.add_argument('-g', '--gravity',   default=False, action='store_true', help='')
+
+    # Get the CLI arguments
+    args = cli.parse_args()
+
+    # Check that at least one file is being merged
+    if not (args.particles or args.gravity):
+        cli.error('At least one of the -p/--particles or -g/--gravity arguments are required.')
+
+    # Set the file names
+    hydro_path       = args.source_directory / '1.h5.0'
+    particle_path    = args.source_directory / '1_particles.h5.0'
+    gravity_path     = args.source_directory / '1_gravity.h5.0'
+    destination_path = args.output_directory / 'combined.h5'
+
+    # Setup the destination file
+    shutil.copy(hydro_path, destination_path)
+
+    # Merge the particle data into the hydro data
+    if args.particles:
+        merge_particles(particle_path, destination_path)
+
+    # Merge the gravity data into the hydro data
+    if args.gravity:
+        merge_gravity(gravity_path, destination_path)
+# =====================================================================================================================
+
+# =====================================================================================================================
+def merge_particles(particle_path: h5py.File, destination_path: pathlib.Path):
+    """Merge the particles data file into the destination (i.e. hydro) data file
+
+    Parameters
+    ----------
+    particle_path : h5py.File
+        The path to the source particles file
+    destination_path : pathlib.Path
+        The path to the destination file
+
+    Raises
+    ------
+    ValueError
+        If a dataset with an identical name exists in both the particles and destination file an error will be raised.
+        The only exception to this is the `density` dataset which exists in both files with different meanings, there
+        is special handling for that and the particle version is renamed to `particle_density`.
+    """
+    # Open the files
+    particles_file   = h5py.File(particle_path,'r')
+    destination_file = h5py.File(destination_path,'r+')
+
+    # Now lets get a list of everything in both source files
+    destination_attr_keys = destination_file.attrs.keys()
+    destination_data_keys = destination_file.keys()
+    particles_attr_keys   = particles_file.attrs.keys()
+    particles_data_keys   = particles_file.keys()
+
+    # Copy all the attributes in the particles file that weren't in the destination file
+    for key in particles_attr_keys:
+        if not key in destination_attr_keys:
+            destination_file.attrs[key] = particles_file.attrs[key]
+
+    # Now we're going to copy all the datasets from the particles file. Note that the "density" dataset requires
+    # special care to avoid duplicating names
+    destination_file.copy(particles_file['density'], 'particle_density')
+    for key in particles_data_keys:
+        if key != 'density':
+            if key not in destination_data_keys:
+                destination_file.copy(particles_file[key], key)
+            else:
+                raise ValueError('Duplicate datasets in destination and particle data files')
+
+    # Close the files
+    particles_file.close()
+    destination_file.close()
+# =====================================================================================================================
+
+# =====================================================================================================================
+def merge_gravity(gravity_path: h5py.File, destination_path: pathlib.Path):
+    """Merge the gravity data file into the destination (i.e. hydro) data file
+
+    Parameters
+    ----------
+    gravity_path : h5py.File
+        The path to the source gravity file
+    destination_path : pathlib.Path
+        The path to the destination file
+
+    Raises
+    ------
+    ValueError
+        If a dataset with an identical name exists in both the gravity and destination file an error will be raised
+    """
+    # Open the files
+    gravity_file   = h5py.File(gravity_path,'r')
+    destination_file = h5py.File(destination_path,'r+')
+
+    # Now lets get a list of everything in both source files
+    destination_attr_keys = destination_file.attrs.keys()
+    destination_data_keys = destination_file.keys()
+    gravity_attr_keys     = gravity_file.attrs.keys()
+    gravity_data_keys     = gravity_file.keys()
+
+    # Copy all the attributes in the particles file that weren't in the destination file
+    for key in gravity_attr_keys:
+        if not key in destination_attr_keys:
+            destination_file.attrs[key] = gravity_file.attrs[key]
+
+    # Now we're going to copy all the datasets from the gravity file.
+    for key in gravity_data_keys:
+        if key not in destination_data_keys:
+            destination_file.copy(gravity_file[key], key)
+        else:
+            raise ValueError('Duplicate datasets in destination and gravity data files')
+
+    # Close the files
+    gravity_file.close()
+    destination_file.close()
+# =====================================================================================================================
+
+# =====================================================================================================================
+if __name__ == '__main__':
+    """This just times the execution of the `main()` function
+    """
+    from timeit import default_timer
+    start = default_timer()
+
+    main()
+
+    print(f'\nTime to execute: {round(default_timer()-start,2)} seconds')
+# =====================================================================================================================

src/analysis/io_analysis.cpp

@@ -656,7 +656,7 @@ void AnalysisModule::Load_Scale_Outputs(struct Parameters *P)
 {
   char filename_1[100];
   strcpy(filename_1, P->analysis_scale_outputs_file);
-  chprintf(" Loading Analysis Scale_Factor Outpus: %s\n", filename_1);
+  chprintf(" Loading Analysis Scale_Factor Outputs: %s\n", filename_1);
 
   ifstream file_out(filename_1);
   string line;

src/cosmology/io_cosmology.cpp

@@ -11,7 +11,7 @@ void Cosmology::Load_Scale_Outputs(struct Parameters *P)
   char filename_1[100];
   // create the filename to read from
   strcpy(filename_1, P->scale_outputs_file);
-  chprintf(" Loading Scale_Factor Outpus: %s\n", filename_1);
+  chprintf(" Loading Scale_Factor Outputs: %s\n", filename_1);
 
   std::ifstream file_out(filename_1);
   std::string line;

src/system_tests/cosmology_system_tests.cpp

@@ -0,0 +1,17 @@
+/*!
+ * \file cosmology_system_tests.cpp
+ * \brief Contains all the system tests for the GRAVITY build type
+ *
+ */
+
+// External Libraries and Headers
+#include <gtest/gtest.h>
+
+// Local includes
+#include "../system_tests/system_tester.h"
+
+TEST(tCOSMOLOGYSYSTEM50Mpc, CorrectInputExpectCorrectOutput)
+{
+  system_test::SystemTestRunner cosmo(true, true, true, true, true);
+  cosmo.runTest(true, 1.0e-07, 0.0006);
+}

src/system_tests/input_files/tCOSMOLOGYSYSTEM50Mpc_CorrectInputExpectCorrectOutput.txt

@@ -0,0 +1,47 @@
+#
+# Parameter File for a 3D Cosmological Simulation.
+#
+
+######################################
+# number of grid cells in the x dimension
+nx=64
+# number of grid cells in the y dimension
+ny=64
+# number of grid cells in the z dimension
+nz=64s
+# output time
+tout=1000
+# how often to output
+outstep=1000
+# value of gamma
+gamma=1.66666667
+# name of initial conditions
+init=Read_Grid
+nfile=0
+#Cosmological Parameters
+H0=67.66
+Omega_M=0.3111
+Omega_L=0.6889
+Omega_b=0.0497
+#Temperature and density floor
+temperature_floor=1e-3
+density_floor=1e-5
+#Outputs file
+scale_outputs_file=./src/system_tests/input_files/tCOSMOLOGYSYSTEM50Mpc_CorrectInputExpectCorrectOutput_a_values.txt
+# domain properties
+xmin=0.0
+ymin=0.0
+zmin=0.0
+xlen=50000.0
+ylen=50000.0
+zlen=50000.0
+# type of boundary conditions
+xl_bcnd=1
+xu_bcnd=1
+yl_bcnd=1
+yu_bcnd=1
+zl_bcnd=1
+zu_bcnd=1
+# path to output directory
+indir=./cholla-tests-data/initial_conditions/tCOSMOLOGYSYSTEM50Mpc_CorrectInputExpectCorrectOutput/
+outdir=./

src/system_tests/input_files/tCOSMOLOGYSYSTEM50Mpc_CorrectInputExpectCorrectOutput_a_values.txt

@@ -0,0 +1,2 @@
+#Single output at z=0 (a=1.0)
+1.0