Realtime reduced data (#82)

* output time and step number to a file rather than screen. Will continue with more meaningful outputs

* move all data reduction including file IO to a separate data reduction class

* added missing files for the DataReducer class

* bring amrex up to 23.02

* added parallel reduction operation over grid. ParReduce is a bit more straightforward than doing the array reduction followed by a parallel reduction

* modify reduction operator to deal with an ArithmeticArray, which is an array that includes basic arithmetic operations

* use code generator to make data reduction work for any number of flavors

* add Tr(rho) to the scalar output. Also only output to the reduction file if the processor id is 0

* generate particle reduction code to work with any number of particles

* added HighFive submodule

* first bit of support for hdf5 realtime io

* added hdf5 support for all reduced data

* added support for outputting the net energy of the system. Required adding one more scalar variable to each particle. This is a good metric to track numerical errors - the number should remain constant

* remove commented reduction examples

* total energy should account for the relative number of real neutrinos present in a particular neutrino particle.

* only output particle info at first timestep if it is output at any time

* add fluxes to realtime data output

* comment out print lines so initial conditions functions can be used more cleanly in scripts

* Reduced max grid size on FFI test to force MPI parallelism

* Run fiducial simulation in test suite. Will use it to test data reduction scripts as well.

* Get rid of mpi errors in stdout

* Reduce the number of steps the fiducial test takes to reduce testing time

* install yt to be able to test data reduction scripts

* Add more data reduction scripts to fiducial test

* First attempt at including HFD5 tests in jenkins

* Remove extraneous plot files to avoid accumuation of 'old' files during testing

* move amrex to submodules directory

* Clean up makefiles to make them a little bit more user friendly and transparent

* fix glob search string so reduction script doesnt match to reduced data files

* fix python script filename

* further reduce number of steps for test speed

* fix some paths to get hdf5 working in tests (hopefully)

* Protect reduced data file initialization to only occur for IO processor. Mixed in with some formatting changes from VS Code.

* cache submodules for faster testing

* ignore hdf5 files

* Create a three-flavor version of the initial condition script.

* Add understandable error message that works even when the code is not compiled in debug mode to save future headache

* Add babysitting plot scripts to plot reduced data output. Add to Jenkins to make sure they work and the reduced data look good.

* add gnuplot library

* Rename one of the makefiles to be a clear default

* draft update of the readme

* line spacing

* A few more tweaks going through the readme and testing

* updated change log for next version

---------

Co-authored-by: Sherwood Richers <[email protected]>

.gitmodules

@@ -1,3 +1,6 @@
 [submodule "amrex"]
-	path = amrex
+	path = submodules/amrex
 	url = https://github.com/dwillcox/amrex.git
+[submodule "submodules/HighFive"]
+	path = submodules/HighFive
+	url = https://github.com/BlueBrain/HighFive.git

CHANGES.md

@@ -1,3 +1,16 @@
+# 1.4
+
+   * Restructured initial conditions to have python scripts generate particle distributions instead of doing so inside the C++ code. This improves the code's flexibility.
+
+   * Add option to output data in HDF5 format. Post-processing scripts only work with the original binary format, since yt only reads the binary format.
+
+   * Add realtime output of scalar quantities to make basic analysis many times faster than with the post-processing scripts.
+
+   * Include all of the basic post-processing scripts with Emu itself to avoid keeping multiple incompatible copies of them.
+
+# 1.3
+
+   * Incorporated various feature additions used for _Code Comparison for Fast Flavor Instability Simulations_ (https://doi.org/10.1103/PhysRevD.106.043011)
 
 # 1.2
 

Dockerfile

@@ -1,7 +1,7 @@
 FROM nvidia/cuda:11.4.0-devel-ubuntu20.04
 ARG DEBIAN_FRONTEND=noninteractive
 RUN apt-get update
-RUN apt-get install -y python3 python3-pip gfortran build-essential libhdf5-openmpi-dev openmpi-bin pkg-config libopenmpi-dev openmpi-bin libblas-dev liblapack-dev libpnetcdf-dev git python-is-python3
-RUN pip3 install numpy matplotlib h5py scipy sympy
+RUN apt-get install -y python3 python3-pip gfortran build-essential libhdf5-openmpi-dev openmpi-bin pkg-config libopenmpi-dev openmpi-bin libblas-dev liblapack-dev libpnetcdf-dev git python-is-python3 gnuplot
+RUN pip3 install numpy matplotlib h5py scipy sympy yt
 ENV USER=jenkins
 ENV LOGNAME=jenkins

Exec/.gitignore

@@ -5,3 +5,6 @@ plt*
 tmp_build_dir
 Backtrace.*
 *.png
+*.dat
+*.h5
+*.pdf

Jenkinsfile

@@ -1,6 +1,11 @@
 pipeline {
     triggers { pollSCM('') }  // Run tests whenever a new commit is detected.
     agent { dockerfile {args '--gpus all'}} // Use the Dockerfile defined in the root Flash-X directory
+    environment {
+		// Get rid of Read -1, expected <someNumber>, errno =1 error
+    	// See https://github.com/open-mpi/ompi/issues/4948
+        OMPI_MCA_btl_vader_single_copy_mechanism = 'none'
+    }
     stages {
 
         //=============================//
@@ -27,13 +32,15 @@ pipeline {
 		sh 'python ../Scripts/initial_conditions/st0_msw_test.py'
 	        sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.CUDA.ex ../sample_inputs/inputs_msw_test'
 	        sh 'python ../Scripts/tests/msw_test.py'
+			sh 'rm -rf plt*'
 	    }
 	}}
 
 	stage('Bipolar'){ steps{
 	    dir('Exec'){
-		sh 'python ../Scripts/initial_conditions/st1_bipolar_test.py'
+			sh 'python ../Scripts/initial_conditions/st1_bipolar_test.py'
 	        sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.CUDA.ex ../sample_inputs/inputs_bipolar_test'
+			sh 'rm -rf plt*'
 	    }
 	}}
 
@@ -42,6 +49,7 @@ pipeline {
 		sh 'python ../Scripts/initial_conditions/st2_2beam_fast_flavor.py'
 	        sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.CUDA.ex ../sample_inputs/inputs_fast_flavor'
 	        sh 'python ../Scripts/tests/fast_flavor_test.py'
+			sh 'rm -rf plt*'
 	    }
 	}}
 
@@ -50,9 +58,38 @@ pipeline {
 		sh 'python ../Scripts/initial_conditions/st3_2beam_fast_flavor_nonzerok.py'
 	        sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.CUDA.ex ../sample_inputs/inputs_fast_flavor_nonzerok'
 	        sh 'python ../Scripts/tests/fast_flavor_k_test.py'
+			sh 'rm -rf plt*'
 	    }
 	}}
 
+	stage('Fiducial 2F GPU Binary'){ steps{
+		dir('Exec'){
+			sh 'python ../Scripts/initial_conditions/st4_linear_moment_ffi.py'
+			sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.CUDA.ex ../sample_inputs/inputs_1d_fiducial'
+			sh 'python ../Scripts/data_reduction/reduce_data.py'
+			sh 'python ../Scripts/data_reduction/reduce_data_fft.py'
+			sh 'python ../Scripts/data_reduction/combine_files.py plt _reduced_data.h5'
+			sh 'python ../Scripts/data_reduction/combine_files.py plt _reduced_data_fft_power.h5'
+			sh 'python ../Scripts/babysitting/avgfee.py'
+			sh 'python ../Scripts/babysitting/power_spectrum.py'
+			sh 'python ../Scripts/data_reduction/convertToHDF5.py'
+			sh 'gnuplot ../Scripts/babysitting/avgfee_gnuplot.plt'
+			archiveArtifacts artifacts: '*.pdf'
+			sh 'rm -rf plt*'
+		}
+	}}
+
+	stage('Fiducial 3F CPU HDF5'){ steps{
+		dir('Exec'){
+	    	sh 'cp ../makefiles/GNUmakefile_jenkins_HDF5 GNUmakefile'
+	        sh 'make realclean; make generate; make -j'
+			sh 'python ../Scripts/initial_conditions/st4_linear_moment_ffi_3F.py'
+			sh 'mpirun -np 4 ./main3d.gnu.TPROF.MPI.ex ../sample_inputs/inputs_1d_fiducial'
+			sh 'python3 ../Scripts/babysitting/avgfee_HDF5.py'
+			sh 'rm -rf plt*'
+		}
+	}}
+
     } // stages{
 
     post {
@@ -62,7 +99,7 @@ pipeline {
 		deleteDirs: true,
 		disableDeferredWipeout: false,
 		notFailBuild: true,
-		patterns: [[pattern: 'amrex', type: 'EXCLUDE']] ) // allow amrex to be cached
+		patterns: [[pattern: 'submodules', type: 'EXCLUDE']] ) // allow submodules to be cached
 	}
     }
 

Make.Emu

@@ -1,6 +1,11 @@
-NUM_FLAVORS ?= 2
-SHAPE_FACTOR_ORDER ?= 2
+# things that used to be defined in the makefile in Exec
 DIM = 3
+TINY_PROFILE  = TRUE
+USE_PARTICLES = TRUE
+PRECISION     = DOUBLE
+Bpack   :=
+Blocs   := .
+
 
 TOP := $(EMU_HOME)
 

README.md

@@ -29,30 +29,40 @@ If you would like to run Emu, first clone Emu with the AMReX submodule:
 
 ```
 git clone --recurse-submodules https://github.com/AMReX-Astro/Emu.git
+git submodule update
 ```
 
-Then change directories to `Emu/Exec`.
-
-Before each compilation, you must symbolically generate Emu source code for
-the number of neutrino flavors you wish to use. Do this like:
+Then change directories to `Emu/Exec`. Before each compilation, you must symbolically generate Emu source code for
+the number of neutrino flavors you wish to use and specify a few other compile-time settings in a file called `GNUmakefile`.
 
+Copy in a default makefile. In this file you can specify the number of neutrino flavors, whether to compile for GPUs, etc. We have set the defaults to 2 neutrino flavors, order 2 PIC shape factors, and compiling for a single CPU.
 ```
-make generate NUM_FLAVORS=2
+cp ../makefiles/GNUmakefile_default GNUmakefile
 ```
 
-Then compile Emu with `make`, e.g.:
-
+Compiling occurs in two stages. We first have to generate code according to the number of neutrino flavors.
+```
+make generate
+```
+Then we have to compile Emu.
 ```
-make NUM_FLAVORS=2
+make -j
 ```
 
-Emu parameters are set in an input file, and we provide a series of sample
-input files for various simulation setups in `Emu/sample_inputs`.
+The initial particle distribution is set by an ASCII particle data file. You can generate the data file with our initial condition scripts. For instance, if we want to simulate a two-beam fast flavor instability, generate the initial conditions using
+```
+python3 ../Scripts/initial_conditions/st3_2beam_fast_flavor_nonzerok.py
+```
+You should now see a new file called `particle_input.dat`.
 
-You can run the MSW setup in Emu by doing:
+The parameters for the simulation are set in input files. These include information about things like the size of the domain, the number of grid cells, and fundamental neutrino properties. Run the fast flavor test simulation using the particle distribution generated previously using one of the test input files stored in `sample_inputs`
+```
+./main3d.gnu.TPROF.ex ../sample_inputs/inputs_fast_flavor_nonzerok
+```
 
+We have a number of data reduction, analysis, and visualization scripts in the `Scripts` directory. Generate a PDF file titled `avgfee.pdf` showing the time evolution of the average number density of electron neutrinos using
 ```
-./main3d.gnu.TPROF.MPI.ex inputs_msw_test
+gnuplot ../Scripts/babysitting/avgfee_gnuplot.plt
 ```
 
 # Open Source Development

Scripts/babysitting/avgfee.py

@@ -1,4 +1,5 @@
-# Run from /ocean/projects/phy200048p/shared to generate plot showing time evolution of <fee> at different dimensionalities
+# plots <N00> and <N_offdiag> as a function of time
+# without reference to the reduced data file outputs
 
 import os
 os.environ['HDF5_USE_FILE_LOCKING'] = 'FALSE'

Scripts/babysitting/avgfee_HDF5.py

@@ -0,0 +1,76 @@
+# plots <N00> and <N_offdiag> as a function of time
+# assuming the code was compiled with HDF5 and wrote the file reduced0D.h5
+
+import os
+os.environ['HDF5_USE_FILE_LOCKING'] = 'FALSE'
+import numpy as np
+import matplotlib.pyplot as plt
+import glob
+import h5py
+import matplotlib as mpl
+from matplotlib.ticker import (MultipleLocator, FormatStrFormatter,AutoMinorLocator,LogLocator)
+
+
+base=["N","Fx","Fy","Fz"]
+diag_flavor=["00","11"]#,"22"]
+offdiag_flavor=["01"]#,"02","12"]
+re=["Re","Im"]
+# real/imag
+R=0
+I=1
+
+
+######################
+# read averaged data #
+######################
+avgData = h5py.File("reduced0D.h5","r")
+t=np.array(avgData["time(s)"])*1e9
+N00=np.array(avgData["N00(1|ccm)"])
+Noffdiag = np.array(avgData["N_offdiag_mag(1|ccm)"])
+avgData.close()
+
+################
+# plot options #
+################
+mpl.rcParams['font.size'] = 22
+mpl.rcParams['font.family'] = 'serif'
+#mpl.rc('text', usetex=True)
+mpl.rcParams['xtick.major.size'] = 7
+mpl.rcParams['xtick.major.width'] = 2
+mpl.rcParams['xtick.major.pad'] = 8
+mpl.rcParams['xtick.minor.size'] = 4
+mpl.rcParams['xtick.minor.width'] = 2
+mpl.rcParams['ytick.major.size'] = 7
+mpl.rcParams['ytick.major.width'] = 2
+mpl.rcParams['ytick.minor.size'] = 4
+mpl.rcParams['ytick.minor.width'] = 2
+mpl.rcParams['axes.linewidth'] = 2
+
+
+fig, ax = plt.subplots(1,1, figsize=(6,5))
+
+##############
+# formatting #
+##############
+ax.axhline(1./3., color="green")
+ax.set_ylabel(r"$\langle N\rangle$ (cm$^{-3}$)")
+ax.set_xlabel(r"$t\,(10^{-9}\,\mathrm{s})$")
+ax.tick_params(axis='both', which='both', direction='in', right=True,top=True)
+ax.xaxis.set_minor_locator(AutoMinorLocator())
+ax.yaxis.set_minor_locator(AutoMinorLocator())
+ax.minorticks_on()
+ax.grid(which='both')
+
+#############
+# plot data #
+#############
+ax.plot(t, N00)
+
+############
+# save pdf #
+############
+plt.savefig("avgfee.pdf", bbox_inches="tight")
+
+# same for f_e\mu
+ax.semilogy(t, Noffdiag)
+plt.savefig("avgfemu.pdf", bbox_inches="tight")

Scripts/babysitting/avgfee_gnuplot.plt

@@ -0,0 +1,7 @@
+set term pdf
+set output 'avgfee.pdf'
+set key autotitle columnhead
+set xlabel 'time (s)'
+set ylabel 'N00 (cm^-3)'
+set yrange [0:*]
+plot 'reduced0D.dat' u 2:5 w l

Scripts/data_reduction/reduce_data_fft.py

@@ -9,7 +9,7 @@
 parser.add_argument("-o", "--output", type=str, default="reduced_data_fft.h5", help="Name of the output file (default: reduced_data_fft.h5)")
 args = parser.parse_args()
 
-directories = sorted(glob.glob("plt*"))
+directories = sorted(glob.glob("plt?????"))
 
 t = []
 

Scripts/initial_conditions/initial_condition_tools.py

@@ -116,7 +116,7 @@ def minerbo_Z(fluxfac):
             print("Failed to converge on a solution.")
             assert(False)
 
-    print("fluxfac=",fluxfac," Z=",Z)
+    #print("fluxfac=",fluxfac," Z=",Z)
     return Z
 
 # angular structure as determined by the Levermore closure
@@ -161,7 +161,7 @@ def levermore_v(fluxfac):
             print("Failed to converge on a solution.")
             assert(False)
 
-    print("fluxfac=",fluxfac," v=",v)
+    #print("fluxfac=",fluxfac," v=",v)
     return v
 
 # interpolate the levermore closure
@@ -251,10 +251,10 @@ def moment_interpolate_particles(nphi_equator, nnu, fnu, energy_erg, direction_g
             # double check that the number densities are correct
             particle_n = np.sum(particles[:,rkey[nvarname]] * particles[:,rkey[fvarname]])
             particle_fmag = np.sum(particles[:,rkey[nvarname]] * particles[:,rkey[fvarname]] * mu[:,nu_nubar, flavor])
-            print("nu/nubar,flavor =", nu_nubar, flavor)
-            print("output/input ndens =",particle_n, nnu[nu_nubar,flavor])
-            print("output/input fluxfac =",particle_fmag / particle_n, fluxfac[nu_nubar,flavor])
-            print()
+            #print("nu/nubar,flavor =", nu_nubar, flavor)
+            #print("output/input ndens =",particle_n, nnu[nu_nubar,flavor])
+            #print("output/input fluxfac =",particle_fmag / particle_n, fluxfac[nu_nubar,flavor])
+            #print()
 
     return particles