diff --git a/.gitmodules b/.gitmodules index c5213734..f2322e41 100644 --- a/.gitmodules +++ b/.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 diff --git a/CHANGES.md b/CHANGES.md index 6017bdfb..bf3777e0 100644 --- a/CHANGES.md +++ b/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 diff --git a/Dockerfile b/Dockerfile index bc578e3c..1c2d49fd 100644 --- a/Dockerfile +++ b/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 diff --git a/Exec/.gitignore b/Exec/.gitignore index e02b640b..7797b25a 100644 --- a/Exec/.gitignore +++ b/Exec/.gitignore @@ -5,3 +5,6 @@ plt* tmp_build_dir Backtrace.* *.png +*.dat +*.h5 +*.pdf diff --git a/Exec/GNUmakefile b/Exec/GNUmakefile deleted file mode 100644 index 3577d5e8..00000000 --- a/Exec/GNUmakefile +++ /dev/null @@ -1,28 +0,0 @@ -EMU_HOME ?= ../ -AMREX_HOME ?= ../amrex - -SHAPE_FACTOR_ORDER ?= 2 -NUM_FLAVORS ?= 3 - -COMP = gnu - -DEBUG = FALSE - -USE_MPI = TRUE -USE_OMP = FALSE -USE_ACC = FALSE -USE_CUDA = FALSE -AMREX_CUDA_ARCH=60 - -USE_HDF5=FALSE -HDF5_HOME=/usr/local/hdf5-1.12.2_gnu9.4.0 - -TINY_PROFILE = TRUE -USE_PARTICLES = TRUE - -PRECISION = DOUBLE - -Bpack := -Blocs := . - -include $(EMU_HOME)/Make.Emu diff --git a/Jenkinsfile b/Jenkinsfile index 4d401eb9..84c48283 100644 --- a/Jenkinsfile +++ b/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 , 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 } } diff --git a/Make.Emu b/Make.Emu index 1b756472..d1081c91 100644 --- a/Make.Emu +++ b/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) diff --git a/README.md b/README.md index d32d487c..bcabfc9c 100644 --- a/README.md +++ b/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 diff --git a/Scripts/babysitting/avgfee.py b/Scripts/babysitting/avgfee.py index 4b17d038..58f6c09e 100755 --- a/Scripts/babysitting/avgfee.py +++ b/Scripts/babysitting/avgfee.py @@ -1,4 +1,5 @@ -# Run from /ocean/projects/phy200048p/shared to generate plot showing time evolution of at different dimensionalities +# plots and as a function of time +# without reference to the reduced data file outputs import os os.environ['HDF5_USE_FILE_LOCKING'] = 'FALSE' diff --git a/Scripts/babysitting/avgfee_HDF5.py b/Scripts/babysitting/avgfee_HDF5.py new file mode 100755 index 00000000..57c2f8eb --- /dev/null +++ b/Scripts/babysitting/avgfee_HDF5.py @@ -0,0 +1,76 @@ +# plots and 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") diff --git a/Scripts/babysitting/avgfee_gnuplot.plt b/Scripts/babysitting/avgfee_gnuplot.plt new file mode 100644 index 00000000..234139eb --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Scripts/data_reduction/reduce_data_fft.py b/Scripts/data_reduction/reduce_data_fft.py index 9a5527b8..7abbe3fe 100755 --- a/Scripts/data_reduction/reduce_data_fft.py +++ b/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 = [] diff --git a/Scripts/initial_conditions/initial_condition_tools.py b/Scripts/initial_conditions/initial_condition_tools.py index 16c0f625..efa8ef6c 100644 --- a/Scripts/initial_conditions/initial_condition_tools.py +++ b/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 diff --git a/Scripts/initial_conditions/st4_linear_moment_ffi_3F.py b/Scripts/initial_conditions/st4_linear_moment_ffi_3F.py new file mode 100644 index 00000000..93bdc96c --- /dev/null +++ b/Scripts/initial_conditions/st4_linear_moment_ffi_3F.py @@ -0,0 +1,45 @@ +import h5py +import numpy as np +import sys +import os +importpath = os.path.dirname(os.path.realpath(__file__)) +sys.path.append(importpath) +sys.path.append(importpath+"/../data_analysis") +from initial_condition_tools import uniform_sphere, write_particles, moment_interpolate_particles, linear_interpolate +import amrex_plot_tools as amrex + +# generation parameters +# MUST MATCH THE INPUTS IN THE EMU INPUT FILE! +NF = 3 +nphi_equator = 16 +theta = 0 +thetabar = 3.14159265359 +phi = 0 +phibar=0 +ndens = 4.891290819e+32 +ndensbar = 4.891290819e+32 +fluxfac = .333333333333333 +fluxfacbar = .333333333333333 +energy_erg = 50 * 1e6*amrex.eV + +# flux factor vectors +fhat = np.array([np.cos(phi) *np.sin(theta ), + np.sin(phi) *np.sin(theta ), + np.cos(theta )]) +fhatbar = np.array([np.cos(phibar)*np.sin(thetabar), + np.sin(phibar)*np.sin(thetabar), + np.cos(thetabar)]) + +nnu = np.zeros((2,NF)) +nnu[0,0] = ndens +nnu[1,0] = ndensbar +nnu[:,1:] = 0 + +fnu = np.zeros((2,NF,3)) +fnu[0,0,:] = ndens * fluxfac * fhat +fnu[1,0,:] = ndensbar * fluxfacbar * fhatbar +fnu[:,1:,:] = 0 + +particles = moment_interpolate_particles(nphi_equator, nnu, fnu, energy_erg, uniform_sphere, linear_interpolate) # [particle, variable] + +write_particles(np.array(particles), NF, "particle_input.dat") diff --git a/Scripts/symbolic_hermitians/generate_code.py b/Scripts/symbolic_hermitians/generate_code.py index b0f00a57..98d536a6 100755 --- a/Scripts/symbolic_hermitians/generate_code.py +++ b/Scripts/symbolic_hermitians/generate_code.py @@ -99,6 +99,7 @@ def delete_generated_files(): for v in vars: A = HermitianMatrix(args.N, v+"{}{}_{}"+t) code += A.header() + code += ["TrHf"] code = [code[i]+"," for i in range(len(code))] write_code(code, os.path.join(args.emu_home, "Source/generated_files", "FlavoredNeutrinoContainer.H_fill")) @@ -115,6 +116,7 @@ def delete_generated_files(): for v in vars: A = HermitianMatrix(args.N, v+"{}{}_{}"+t) code += A.header() + code += ["TrHf"] code_string = 'attribute_names = {"time", "x", "y", "z", "pupx", "pupy", "pupz", "pupt", ' code = ['"{}"'.format(c) for c in code] code_string = code_string + ", ".join(code) + "};" @@ -168,6 +170,51 @@ def delete_generated_files(): code.append(string1+deplist[icomp]+string2+flist[icomp]+string3+string4[ivar]) write_code(code, os.path.join(args.emu_home, "Source/generated_files", "Evolve.cpp_deposit_to_mesh_fill")) + #================================# + # DataReducer.cpp_fill_particles # + #================================# + tails = ["","bar"] + code = [] + for t in tails: + # diagonal averages + N = HermitianMatrix(args.N, "p.rdata(PIdx::f{}{}_{}"+t+")") + Nlist = N.header_diagonals(); + for i in range(len(Nlist)): + code.append("Trf += "+Nlist[i]+";") + + write_code(code, os.path.join(args.emu_home, "Source/generated_files", "DataReducer.cpp_fill_particles")) + + #======================# + # DataReducer.cpp_fill # + #======================# + tails = ["","bar"] + code = [] + for t in tails: + # diagonal averages + N = HermitianMatrix(args.N, "a(i\,j\,k\,GIdx::N{}{}_{}"+t+")") + Nlist = N.header_diagonals(); + Fx = HermitianMatrix(args.N, "a(i\,j\,k\,GIdx::Fx{}{}_{}"+t+")") + Fxlist = Fx.header_diagonals(); + Fy = HermitianMatrix(args.N, "a(i\,j\,k\,GIdx::Fy{}{}_{}"+t+")") + Fylist = Fy.header_diagonals(); + Fz = HermitianMatrix(args.N, "a(i\,j\,k\,GIdx::Fz{}{}_{}"+t+")") + Fzlist = Fz.header_diagonals(); + for i in range(len(Nlist)): + code.append("Ndiag"+t+"["+str(i)+"] = "+Nlist[i]+";") + code.append("Fxdiag"+t+"["+str(i)+"] = "+Fxlist[i]+";") + code.append("Fydiag"+t+"["+str(i)+"] = "+Fylist[i]+";") + code.append("Fzdiag"+t+"["+str(i)+"] = "+Fzlist[i]+";") + + # off-diagonal magnitude + mag2 = 0 + for i in range(N.size): + for j in range(i+1,N.size): + re,im = N.H[i,j].as_real_imag() + mag2 += re**2 + im**2 + code.append("N_offdiag_mag2 += "+sympy.cxxcode(sympy.simplify(mag2))+";") + + write_code(code, os.path.join(args.emu_home, "Source/generated_files", "DataReducer.cpp_fill")) + #==================# # Evolve.H_M2_fill # #==================# @@ -384,6 +431,11 @@ def sgn(t,var): # Write out dFdt->F code.append(dFdt.code()) + + # store Tr(H*F) for estimating numerical errors + TrHf = (H*F).trace(); + code.append(["p.rdata(PIdx::TrHf) += p.rdata(PIdx::N"+t+") * ("+sympy.cxxcode(sympy.simplify(TrHf))+");"]) + code = [line for sublist in code for line in sublist] write_code(code, os.path.join(args.emu_home, "Source/generated_files", "Evolve.cpp_dfdt_fill")) diff --git a/Source/ArithmeticArray.H b/Source/ArithmeticArray.H new file mode 100644 index 00000000..36d32403 --- /dev/null +++ b/Source/ArithmeticArray.H @@ -0,0 +1,93 @@ +template +class ArithmeticArray{ +public: + std::array data; + + AMREX_GPU_HOST_DEVICE + ArithmeticArray(){} + + AMREX_GPU_HOST_DEVICE + ArithmeticArray(T in){ + for(size_t i=0; i + AMREX_GPU_HOST_DEVICE + ArithmeticArray(ArithmeticArray in){ + for(size_t i=0; i + AMREX_GPU_HOST_DEVICE + ArithmeticArray& operator=(const ArithmeticArray& input){ + for(size_t i=0; i& operator=(const double input){ + for(size_t i=0; i operator*(const double scale) const{ + ArithmeticArray result; + for(size_t i=0; i operator/(const double scale) const{ + double inv_scale = 1./scale; + return operator*(inv_scale); + } + template + AMREX_GPU_HOST_DEVICE + const ArithmeticArray operator+(const ArithmeticArray& input) const{ + ArithmeticArray result; + for(size_t i=0; i + AMREX_GPU_HOST_DEVICE + const ArithmeticArray operator-(const ArithmeticArray& input) const{ + ArithmeticArray result; + for(size_t i=0; i operator-() const{ + ArithmeticArray result; + for(size_t i=0; i + AMREX_GPU_HOST_DEVICE + void operator+=(const ArithmeticArray& input){ + for(size_t i=0; i + AMREX_GPU_HOST_DEVICE + bool operator==(const ArithmeticArray& input){ + bool isequal = true; + for(size_t i=0; i +#include +#include +#include +#include +#include +#ifdef AMREX_USE_HDF5 +#include <../submodules/HighFive/include/highfive/H5File.hpp> +#include <../submodules/HighFive/include/highfive/H5DataSpace.hpp> +#include <../submodules/HighFive/include/highfive/H5DataSet.hpp> +#endif + +class DataReducer{ +private: + +#ifdef AMREX_USE_HDF5 + const char* filename0D = "reduced0D.h5"; + const HighFive::DataSpace dataspace = HighFive::DataSpace({0}, {HighFive::DataSpace::UNLIMITED}); +#else + const char* filename0D = "reduced0D.dat"; +#endif + +public: + void + InitializeFiles(); + + void + WriteReducedData0D(const amrex::Geometry& geom, + const MultiFab& state, + const FlavoredNeutrinoContainer& neutrinos, + amrex::Real time, int step); +}; + +#endif diff --git a/Source/DataReducer.cpp b/Source/DataReducer.cpp new file mode 100644 index 00000000..cb283d11 --- /dev/null +++ b/Source/DataReducer.cpp @@ -0,0 +1,219 @@ +#include "Evolve.H" +#include "Constants.H" +#include "DataReducer.H" +#include "ArithmeticArray.H" +#include +#include +#ifdef AMREX_USE_HDF5 +#include <../submodules/HighFive/include/highfive/H5File.hpp> +#include <../submodules/HighFive/include/highfive/H5DataSpace.hpp> +#include <../submodules/HighFive/include/highfive/H5DataSet.hpp> +#endif + +#ifdef AMREX_USE_HDF5 +// append a single scalar to the specified file and dataset +template +void append_0D(HighFive::File& file0D, const std::string& datasetname, const T value){ + HighFive::DataSet dataset_step = file0D.getDataSet(datasetname); + std::vector dims = dataset_step.getDimensions(); + dims[0] ++; + dataset_step.resize(dims); + dataset_step.select({dims[0]-1},{1}).write((T[1]){value}); +} +#endif + +void DataReducer::InitializeFiles() +{ + if (ParallelDescriptor::IOProcessor()) + { + +#ifdef AMREX_USE_HDF5 + + using namespace HighFive; + File file0D(filename0D, File::Truncate | File::Create); + + DataSetCreateProps props; + props.add(Chunking(std::vector{1})); + file0D.createDataSet("step", dataspace, create_datatype(), props); + file0D.createDataSet("time(s)", dataspace, create_datatype(), props); + file0D.createDataSet("Ntot(1|ccm)", dataspace, create_datatype(), props); + file0D.createDataSet("Ndiff(1|ccm)", dataspace, create_datatype(), props); + for (int i = 0; i < NUM_FLAVORS; i++) + { + file0D.createDataSet(std::string("N") + std::to_string(i) + std::to_string(i) + std::string("(1|ccm)"), dataspace, create_datatype(), props); + file0D.createDataSet(std::string("N") + std::to_string(i) + std::to_string(i) + std::string("bar(1|ccm)"), dataspace, create_datatype(), props); + file0D.createDataSet(std::string("Fx") + std::to_string(i) + std::to_string(i) + std::string("(1|ccm)"), dataspace, create_datatype(), props); + file0D.createDataSet(std::string("Fy") + std::to_string(i) + std::to_string(i) + std::string("(1|ccm)"), dataspace, create_datatype(), props); + file0D.createDataSet(std::string("Fz") + std::to_string(i) + std::to_string(i) + std::string("(1|ccm)"), dataspace, create_datatype(), props); + file0D.createDataSet(std::string("Fx") + std::to_string(i) + std::to_string(i) + std::string("bar(1|ccm)"), dataspace, create_datatype(), props); + file0D.createDataSet(std::string("Fy") + std::to_string(i) + std::to_string(i) + std::string("bar(1|ccm)"), dataspace, create_datatype(), props); + file0D.createDataSet(std::string("Fz") + std::to_string(i) + std::to_string(i) + std::string("bar(1|ccm)"), dataspace, create_datatype(), props); + } + file0D.createDataSet("N_offdiag_mag(1|ccm)", dataspace, create_datatype(), props); + file0D.createDataSet("sumTrf", dataspace, create_datatype(), props); + file0D.createDataSet("sumTrHf", dataspace, create_datatype(), props); + +#else + + std::ofstream outfile; + outfile.open(filename0D, std::ofstream::out); + int j = 0; + j++; + outfile << j << ":step\t"; + j++; + outfile << j << ":time(s)\t"; + j++; + outfile << j << ":Ntot(1|ccm)\t"; + j++; + outfile << j << ":Ndiff(1|ccm)\t"; + for (int i = 0; i < NUM_FLAVORS; i++) + { + j++; + outfile << j << ":N" << i << i << "(1|ccm)\t"; + j++; + outfile << j << ":N" << i << i << "bar(1|ccm)\t"; + j++; + outfile << j << ":Fx" << i << i << "(1|ccm)\t"; + j++; + outfile << j << ":Fy" << i << i << "(1|ccm)\t"; + j++; + outfile << j << ":Fz" << i << i << "(1|ccm)\t"; + j++; + outfile << j << ":Fx" << i << i << "bar(1|ccm)\t"; + j++; + outfile << j << ":Fy" << i << i << "bar(1|ccm)\t"; + j++; + outfile << j << ":Fz" << i << i << "bar(1|ccm)\t"; + } + j++; + outfile << j << ":N_offdiag_mag(1|ccm)\t"; + j++; + outfile << j << ":sumTrf\t"; + j++; + outfile << j << ":sumTrHf\t"; + outfile << std::endl; + outfile.close(); + +#endif + } +} + +void +DataReducer::WriteReducedData0D(const amrex::Geometry& geom, + const MultiFab& state, + const FlavoredNeutrinoContainer& neutrinos, + const amrex::Real time, const int step) +{ + // get index volume of the domain + int ncells = geom.Domain().volume(); + + //==================================// + // Do reductions over the particles // + //==================================// + using PType = typename FlavoredNeutrinoContainer::ParticleType; + amrex::ReduceOps reduce_ops; + auto particleResult = amrex::ParticleReduce< ReduceData >(neutrinos, + [=] AMREX_GPU_DEVICE(const PType& p) noexcept -> amrex::GpuTuple { + Real TrHf = p.rdata(PIdx::TrHf); + Real Trf = 0; +#include "generated_files/DataReducer.cpp_fill_particles" + return GpuTuple{Trf,TrHf}; + }, reduce_ops); + Real Trf = amrex::get<0>(particleResult); + Real TrHf = amrex::get<1>(particleResult); + ParallelDescriptor::ReduceRealSum(Trf); + ParallelDescriptor::ReduceRealSum(TrHf); + + //=============================// + // Do reductions over the grid // + //=============================// + // first, get a reference to the data arrays + auto const& ma = state.const_arrays(); + IntVect nghost(AMREX_D_DECL(0, 0, 0)); + + // use the ParReduce function to define reduction operator + GpuTuple< ArithmeticArray, ArithmeticArray, Real , ArithmeticArray, ArithmeticArray, ArithmeticArray, ArithmeticArray, ArithmeticArray, ArithmeticArray > result = + ParReduce(TypeList{}, + TypeList, ArithmeticArray, Real , ArithmeticArray, ArithmeticArray, ArithmeticArray, ArithmeticArray, ArithmeticArray, ArithmeticArray >{}, + state, nghost, + [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept -> + GpuTuple, ArithmeticArray, Real , ArithmeticArray, ArithmeticArray, ArithmeticArray, ArithmeticArray, ArithmeticArray, ArithmeticArray > { + Array4 const& a = ma[box_no]; + + // Doing the actual work + ArithmeticArray Ndiag, Ndiagbar; + ArithmeticArray Fxdiag, Fxdiagbar; + ArithmeticArray Fydiag, Fydiagbar; + ArithmeticArray Fzdiag, Fzdiagbar; + Real N_offdiag_mag2 = 0; + #include "generated_files/DataReducer.cpp_fill" + return {Ndiag, Ndiagbar, N_offdiag_mag2, Fxdiag, Fydiag, Fzdiag, Fxdiagbar,Fydiagbar,Fzdiagbar}; + + }); + + // retrieve the reduced data values + ArithmeticArray N = amrex::get<0>(result) / ncells; + ArithmeticArray Nbar = amrex::get<1>(result) / ncells; + Real N_offdiag_mag = sqrt(amrex::get<2>(result)) / ncells; + ArithmeticArray Fx = amrex::get<3>(result) / ncells; + ArithmeticArray Fy = amrex::get<4>(result) / ncells; + ArithmeticArray Fz = amrex::get<5>(result) / ncells; + ArithmeticArray Fxbar = amrex::get<6>(result) / ncells; + ArithmeticArray Fybar = amrex::get<7>(result) / ncells; + ArithmeticArray Fzbar = amrex::get<8>(result) / ncells; + + // calculate net number of neutrinos and antineutrinos + Real Ntot=0, Ndiff=0; + for(int i=0; i > read_particle_data(std::strin ss = std::stringstream(line); int NF_in; ss >> NF_in; + if(NF_in != NUM_FLAVORS) amrex::Print() << "Error: number of flavors in particle data file does not match the number of flavors Emu was compiled for." << std::endl; AMREX_ASSERT(NF_in == NUM_FLAVORS); while(std::getline(file, line)){ diff --git a/Source/IO.H b/Source/IO.H index e851f0a2..c159b065 100644 --- a/Source/IO.H +++ b/Source/IO.H @@ -6,6 +6,15 @@ #include #include +void +InitializeReducedData0D(); + +void +WriteReducedData0D(const amrex::MultiFab& state, + const FlavoredNeutrinoContainer& neutrinos, + const amrex::Geometry& geom, amrex::Real time, + int step); + void WritePlotFile (const amrex::MultiFab& state, const FlavoredNeutrinoContainer& neutrinos, diff --git a/Source/Make.package b/Source/Make.package index 191e7a7f..b038507d 100644 --- a/Source/Make.package +++ b/Source/Make.package @@ -1,4 +1,5 @@ CEXE_sources += main.cpp +CEXE_sources += DataReducer.cpp CEXE_sources += IO.cpp CEXE_sources += FlavoredNeutrinoContainerInit.cpp CEXE_sources += Evolve.cpp @@ -11,3 +12,5 @@ CEXE_headers += Particles.H CEXE_headers += IO.H CEXE_headers += Parameters.H CEXE_headers += ParticleInterpolator.H +CEXE_headers += DataReducer.H +CEXE_headers += ArithmeticArray.H \ No newline at end of file diff --git a/Source/main.cpp b/Source/main.cpp index b17266da..0b6eac4b 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -31,6 +31,7 @@ #include "Evolve.H" #include "Constants.H" #include "IO.H" +#include "DataReducer.H" using namespace amrex; @@ -118,10 +119,12 @@ void evolve_flavor(const TestParams* parms) deposit_to_mesh(neutrinos_old, state, geom); // Write plotfile after initialization + DataReducer rd; if (not parms->do_restart) { // If we have just initialized, then always save the particle data for reference - const int write_particles_after_init = 1; + const int write_particles_after_init = (parms->write_plot_particles_every>0); WritePlotFile(state, neutrinos_old, geom, initial_time, initial_step, write_particles_after_init); + rd.InitializeFiles(); } amrex::Print() << "Done. " << std::endl; @@ -176,7 +179,7 @@ void evolve_flavor(const TestParams* parms) const int step = integrator.get_step_number(); const Real time = integrator.get_time(); - amrex::Print() << "Completed time step: " << step << " t = " << time << " s. ct = " << PhysConst::c * time << " cm" << std::endl; + rd.WriteReducedData0D(geom, state, neutrinos, time, step+1); run_fom += neutrinos.TotalNumberOfParticles(); diff --git a/amrex b/amrex deleted file mode 160000 index ff9e8324..00000000 --- a/amrex +++ /dev/null @@ -1 +0,0 @@ -Subproject commit ff9e832483f142a34da72f9bea5fed72e49fea33 diff --git a/makefiles/GNUmakefile_cassowary b/makefiles/GNUmakefile_cassowary new file mode 100644 index 00000000..4ed76b6f --- /dev/null +++ b/makefiles/GNUmakefile_cassowary @@ -0,0 +1,19 @@ +NUM_FLAVORS = 2 +SHAPE_FACTOR_ORDER = 2 + +COMP = gnu + +DEBUG = FALSE + +USE_MPI = TRUE +USE_OMP = FALSE +USE_ACC = FALSE +USE_CUDA = FALSE +AMREX_CUDA_ARCH=60 + +USE_HDF5=FALSE +HDF5_HOME=/usr/lib/x86_64-linux-gnu/hdf5/serial + +EMU_HOME = .. +AMREX_HOME = ../submodules/amrex +include ../Make.Emu diff --git a/makefiles/GNUmakefile_default b/makefiles/GNUmakefile_default new file mode 100644 index 00000000..74d83b04 --- /dev/null +++ b/makefiles/GNUmakefile_default @@ -0,0 +1,16 @@ +NUM_FLAVORS = 2 +SHAPE_FACTOR_ORDER = 2 + +COMP = gnu + +DEBUG = FALSE + +USE_MPI = FALSE +USE_OMP = FALSE +USE_ACC = FALSE +USE_CUDA = FALSE +USE_HDF5 = FALSE + +EMU_HOME = .. +AMREX_HOME = ../submodules/amrex +include ../Make.Emu diff --git a/makefiles/GNUmakefile_ganon b/makefiles/GNUmakefile_ganon deleted file mode 100644 index e6c6fa45..00000000 --- a/makefiles/GNUmakefile_ganon +++ /dev/null @@ -1,33 +0,0 @@ -EMU_HOME ?= ../ -AMREX_HOME ?= ../amrex - -DIM = 3 - -NUM_FLAVORS = 2 - -COMP = gnu - -DEBUG = TRUE - -USE_MPI = TRUE -USE_OMP = FALSE -USE_ACC = FALSE -USE_CUDA = FALSE -USE_HDF5 = FALSE - -TINY_PROFILE = TRUE -USE_PARTICLES = TRUE - -PRECISION = DOUBLE - -Bpack := -Blocs := . - -ifeq ($(USE_HDF5), TRUE) -HDF5_HOME = /usr/local/hdf5-1.12.0_gnu7.5.0 -DEFINES += -DAMREX_USE_HDF5 -INCLUDE_LOCATIONS += $(HDF5_HOME)/include -LIBRARIES += -L$(HDF5_HOME)/lib -lhdf5 -lz -ldl -endif - -include ../Make.Emu diff --git a/makefiles/GNUmakefile_jenkins b/makefiles/GNUmakefile_jenkins index da1d700d..431901c8 100644 --- a/makefiles/GNUmakefile_jenkins +++ b/makefiles/GNUmakefile_jenkins @@ -1,15 +1,9 @@ -EMU_HOME ?= ../ -AMREX_HOME ?= ../amrex - -DIM = 3 - NUM_FLAVORS = 2 - SHAPE_FACTOR_ORDER = 2 COMP = gnu -DEBUG = FALSE +DEBUG = FALSE USE_MPI = TRUE USE_OMP = FALSE @@ -17,12 +11,6 @@ USE_ACC = FALSE USE_CUDA = TRUE AMREX_CUDA_ARCH=75 -TINY_PROFILE = TRUE -USE_PARTICLES = TRUE - -PRECISION = DOUBLE - -Bpack := -Blocs := . - +EMU_HOME = .. +AMREX_HOME = ../submodules/amrex include ../Make.Emu diff --git a/makefiles/GNUmakefile_jenkins_HDF5 b/makefiles/GNUmakefile_jenkins_HDF5 new file mode 100644 index 00000000..1a5b3993 --- /dev/null +++ b/makefiles/GNUmakefile_jenkins_HDF5 @@ -0,0 +1,19 @@ +NUM_FLAVORS = 3 +SHAPE_FACTOR_ORDER = 2 + +COMP = gnu + +DEBUG = FALSE + +USE_MPI = TRUE +USE_OMP = FALSE +USE_ACC = FALSE +USE_CUDA = FALSE +AMREX_CUDA_ARCH=75 + +USE_HDF5=TRUE +HDF5_HOME=/usr/lib/x86_64-linux-gnu/hdf5/openmpi + +EMU_HOME = .. +AMREX_HOME = ../submodules/amrex +include ../Make.Emu diff --git a/sample_inputs/inputs_1d_fiducial b/sample_inputs/inputs_1d_fiducial index 833421e8..c3e89d37 100644 --- a/sample_inputs/inputs_1d_fiducial +++ b/sample_inputs/inputs_1d_fiducial @@ -23,7 +23,7 @@ particle_data_filename = "particle_input.dat" max_grid_size = 16 # Number of steps to run -nsteps = 5000 +nsteps = 1000 # Simulation end time end_time = 5.0e-9 diff --git a/sample_inputs/inputs_fast_flavor_nonzerok b/sample_inputs/inputs_fast_flavor_nonzerok index 1996e71a..dc858a8e 100644 --- a/sample_inputs/inputs_fast_flavor_nonzerok +++ b/sample_inputs/inputs_fast_flavor_nonzerok @@ -21,7 +21,7 @@ nppc = (1, 1, 1) particle_data_filename = "particle_input.dat" # Maximum size of each grid in the domain -max_grid_size = 1000 +max_grid_size = 10 # Number of steps to run nsteps = 5000 diff --git a/submodules/HighFive b/submodules/HighFive new file mode 160000 index 00000000..4d29deeb --- /dev/null +++ b/submodules/HighFive @@ -0,0 +1 @@ +Subproject commit 4d29deebf939ba47f8f28f0594b813fd9289c0fd diff --git a/submodules/amrex b/submodules/amrex new file mode 160000 index 00000000..eefe2463 --- /dev/null +++ b/submodules/amrex @@ -0,0 +1 @@ +Subproject commit eefe2463884081a27025973d4a99dc909657b4f0