From 5ee3465a9373e17f5698a6fa33b04d885c0a712b Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 19 Nov 2024 20:02:08 -0700 Subject: [PATCH] add eos_grid example --- example/CMakeLists.txt | 6 ++ example/eos_grid.cpp | 165 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 171 insertions(+) create mode 100644 example/eos_grid.cpp diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt index c367b8b51e..c070182ab1 100644 --- a/example/CMakeLists.txt +++ b/example/CMakeLists.txt @@ -17,6 +17,12 @@ add_executable(get_sound_speed_press target_link_libraries(get_sound_speed_press PRIVATE singularity-eos::singularity-eos) +if (SINGULARITY_USE_SPINER_WITH_HDF5) + add_executable(eos_grid eos_grid.cpp) + target_link_libraries(eos_grid PRIVATE + singularity-eos::singularity-eos) +endif() + if(SINGULARITY_USE_EOSPAC AND SINGULARITY_USE_SPINER_WITH_HDF5) add_executable(get_sesame_state get_sesame_state.cpp) diff --git a/example/eos_grid.cpp b/example/eos_grid.cpp new file mode 100644 index 0000000000..e89ab6a117 --- /dev/null +++ b/example/eos_grid.cpp @@ -0,0 +1,165 @@ +//------------------------------------------------------------------------------ +// © 2024. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National +// Nuclear Security Administration. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//------------------------------------------------------------------------------ + +// C headers +#include + +// C++ headers +#include +#include +#include + +// HDF5 we'll need this for I/O +#include +#include + +// This library contains portable utilities +#include + +// This contains useful tools for preventing things like divide by zero +#include +// Needed to import the eos models +#include + +// This library contains the spiner table object, which we will use to +// store our output +#include +#include + +// for timing +using duration = std::chrono::nanoseconds; + +// These are the specializations of spiner we will use +using DataBox = Spiner::DataBox; +using RegularGrid1D = Spiner::RegularGrid1D; + +// Set the EOS you want to use here. +using EOS = singularity::StellarCollapse; + +int main(int argc, char *argv[]) { + // This is needed for Kokkos +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::initialize(); +#endif + // note the scoping here. This means Kokkos objects are cleaned up + // before finalization + { + if (argc < 8) { + std::cerr << "Usage: " << argv[0] + << " savefilename log10rhomin log10rhomax numrho log10Tmin log10Tmax " + "numT eosargs..." + << std::endl; + std::exit(1); + } + const std::string savefilename = argv[1]; + const double lrho_min = atof(argv[2]); + const double lrho_max = atof(argv[3]); + const int nrho = atoi(argv[4]); + RegularGrid1D lrhoGrid(lrho_min, lrho_max, nrho); + + const double lT_min = atof(argv[5]); + const double lT_max = atof(argv[6]); + const int nT = atoi(argv[7]); + RegularGrid1D lTGrid(lT_min, lT_max, nT); + + // This is the databox we will evaluate onto + DataBox press_d(Spiner::AllocationTarget::Device, nrho, nT); + press_d.setRange(0, lTGrid); // note 0 is rightmost index + press_d.setRange(1, lrhoGrid); + // This is the host mirror, we will copy into it + DataBox press_h(Spiner::AllocationTarget::Host, nrho, nT); + press_h.setRange(0, lTGrid); + press_h.setRange(1, lrhoGrid); + + // These are the pre-evaluated density and temperature points + std::cout << "Filling grid points..." << std::endl; + DataBox rhos(Spiner::AllocationTarget::Device, nrho); + portableFor( + "Set rho", 0, nrho, + PORTABLE_LAMBDA(const int i) { rhos(i) = std::pow(10., lrhoGrid.x(i)); }); + DataBox Ts(Spiner::AllocationTarget::Device, nT); + portableFor( + "Set T", 0, nT, + PORTABLE_LAMBDA(const int i) { Ts(i) = std::pow(10., lTGrid.x(i)); }); + + // if you have arguments you want to pass into your EOS load them here + const std::string loadname = argv[8]; + const double Ye = atof(argv[9]); + + std::cout << "Initializing EOS..." << std::endl; + EOS eos_h(loadname); + EOS eos_d = eos_h.GetOnDevice(); + + std::cout << "Evaluating..." << std::endl; + // start timers +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); +#endif + auto start = std::chrono::high_resolution_clock::now(); + portableFor( + "P(rho, T)", 0, nrho, 0, nT, PORTABLE_LAMBDA(const int j, const int i) { + // Some EOS objects take lambdas. If you want to set + // them locally, do it like this. Otherwise, you may + // want to create a device-side array for them to + // avoid race conditions. + Real lambda[2]; + lambda[0] = Ye; + press_d(j, i) = eos_d.PressureFromDensityTemperature(rhos(j), Ts(i), lambda); + }); + // stop timers +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); +#endif + auto stop = std::chrono::high_resolution_clock::now(); + auto duration_device = std::chrono::duration_cast(stop - start); + + // copy data from device to host + std::cout << "Copying to host..." << std::endl; + portableCopyToHost(press_h.data(), press_d.data(), press_h.sizeBytes()); + + // For fun lets also evalaute all on host + std::cout << "Generate arrays on host..." << std::endl; + DataBox rhos_h(Spiner::AllocationTarget::Host, nrho); + for (int i = 0; i < nrho; ++i) { + rhos_h(i) = std::pow(10., lrhoGrid.x(i)); + } + DataBox Ts_h(Spiner::AllocationTarget::Host, nT); + for (int i = 0; i < nT; ++i) { + Ts_h(i) = std::pow(10., lTGrid.x(i)); + } + std::cout << "Looping on host" << std::endl; + start = std::chrono::high_resolution_clock::now(); + for (int j = 0; j < nT; ++j) { + for (int i = 0; i < nrho; ++i) { + Real lambda[2]; + lambda[0] = Ye; + press_h(j, i) = eos_h.PressureFromDensityTemperature(rhos_h(j), Ts_h(i), lambda); + } + } + stop = std::chrono::high_resolution_clock::now(); + auto duration_host = std::chrono::duration_cast(stop - start); + + std::cout << "Saving file..." << std::endl; + press_h.saveHDF(savefilename); // if you have HDF5 enabled, thats all it takes + std::cout << "Time per point, (ns):\n" + << "host = " << duration_host.count() / (static_cast(nrho) * nT) + << "\n" + << "device = " << duration_device.count() / (static_cast(nrho) * nT) + << std::endl; + } +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::finalize(); +#endif + return 0; +}