diff --git a/.gitignore b/.gitignore new file mode 100644 index 000000000..c70fdda43 --- /dev/null +++ b/.gitignore @@ -0,0 +1,10 @@ +a.out +*.o +test.bin +convergence.bin +*.hpp.gch +convergence.dat +.markdown-preview.html +*.sp5 +*.lst +*.s diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 000000000..f4cabc62c --- /dev/null +++ b/.gitmodules @@ -0,0 +1,4 @@ +[submodule "Catch2"] + path = Catch2 + url = https://github.com/catchorg/Catch2.git + branch = v2.x diff --git a/Catch2 b/Catch2 new file mode 160000 index 000000000..68975e3ff --- /dev/null +++ b/Catch2 @@ -0,0 +1 @@ +Subproject commit 68975e3ff3a9d026965ad142b04f548e685b5095 diff --git a/Makefile b/Makefile new file mode 100644 index 000000000..f1110d9dc --- /dev/null +++ b/Makefile @@ -0,0 +1,71 @@ +# © (or copyright) 2019-2021. 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. + +USE_HDF5=no + +CC=g++ +INC= spiner_types.hpp \ + databox.hpp \ + interpolation.hpp +GCH=$(INC:.hpp=.hpp.gch) + +H5CC=h5c++ +H5_DIR=/usr/local/hdf5-parallel +H5_INCLUDE=-I${H5_DIR}/include +H5_LIB=-L${H5_DIR}/lib +H5_LINK=-lhdf5_hl -lhdf5 -Wl,-rpath=${H5_DIR}/lib/ + +CATCH_INCLUDE=-ICatch2/single_include/catch2 +PORTS_INCLUDE=-Iports-of-call +BASE_FLAGS=-std=c++14 -Wall -fdiagnostics-color=always -g + +ifeq ($(USE_HDF5),yes) + CC=${H5CC} + INCLUDE_FLAGS=${CATCH_INCLUDE} ${PORTS_INCLUDE} ${H5_INCLUDE} + LIB_FLAGS=${H5_LIB} + LINK_FLAGS=${H5_LINK} + CFLAGS=${BASE_FLAGS} ${INCLUDE_FLAGS} -DSPINER_USE_HDF +else + INCLUDE_FLAGS=${CATCH_INCLUDE} ${PORTS_INCLUDE} + LIB_FLAGS= + LINK_FLAGS= + CFLAGS=${BASE_FLAGS} ${INCLUDE_FLAGS} +endif + +LFLAGS=${LIB_FLAGS} ${LINK_FLAGS} + +default: test + +all: test convergence + +test: test.bin + ./test.bin + +convergence: convergence.bin + ./convergence.bin + python plot_convergence.py + +%.bin: %.o + $(CC) ${CFLAGS} ${LFLAGS} -g -o $@ $^ + +%.o: %.cpp ${INC} + $(CC) ${CFLAGS} -c -o $@ $< + +.PHONY: default all test convergence clean info + +clean: + $(RM) test.bin convergence.bin + $(RM) test.o convergence.o databox.o + $(RM) ${GCH} + $(RM) convergence.dat + $(RM) databox.sp5 diff --git a/Makefile.kokkos b/Makefile.kokkos new file mode 100644 index 000000000..5d41bca23 --- /dev/null +++ b/Makefile.kokkos @@ -0,0 +1,78 @@ +# © (or copyright) 2019-2021. 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. + +KOKKOS_PATH = ${HOME}/kokkos +#SRC = test.cpp +SRC = convergence.cpp +vpath %.cpp $(sort $(dir $(SRC))) + +KOKKOS_CXX_STANDARD=c++14 + +default: build + echo "Start Build" + +#H5_DIR=/usr/local/hdf5-parallel +#H5_INCLUDE=-I${H5_DIR}/include +#H5_LIB=-L${H5_DIR}/lib +#H5_LINK=-lhdf5_hl -lhdf5 -Wl,-rpath=${H5_DIR}/lib/ + +ifneq (,$(findstring Cuda,$(KOKKOS_DEVICES))) +CXX = ${KOKKOS_PATH}/bin/nvcc_wrapper +CXXFLAGS = -O3 --expt-relaxed-constexpr +CXXFLAGS += -DPORTABILITY_STRATEGY_KOKKOS #-DSPINER_USE_HDF +LINK = ${CXX} +#LDFLAGS = +#EXE = simple_test_cuda +KOKKOS_DEVICES = "Cuda" +KOKKOS_ARCH = "Volta70" +KOKKOS_CUDA_OPTIONS += "enable_lambda;rdc" +#CXXFLAGS = -O3 --expt-relaxed-constexpr #-x cuda --cuda-gpu-arch=sm_37 --stdlib=libstdc++ -lineinfo +LINK = ${CXX} +LDFLAGS = +EXE = simple_test_cuda +else +CXX = g++ +CXXFLAGS = -O3 -g -DPORTABILITY_STRATEGY_KOKKOS #-DSPINER_USE_HDF +LINK = ${CXX} +LDFLAGS = #-lhdf5_hl -lhdf5 -Wl,-rpath=${H5_DIR}/lib/ +#EXE = simple_test +EXE = convergence_test +KOKKOS_DEVICES = "OpenMP" +KOKKOS_ARCH = "" +#KOKKOS_DEBUG=yes +endif + +EXTRA_INC += -ICatch2/single_include/catch2 -Iports-of-call +#EXTRA_INC += ${H5_INCLUDE} +DEPFLAGS = -M + +OBJ = $(notdir $(SRC:.cpp=.o)) +#LIB = -L${H5_DIR}/lib + +include $(KOKKOS_PATH)/Makefile.kokkos + +build: $(EXE) + +test: $(EXE) + ./$(EXE) + +$(EXE): $(OBJ) $(KOKKOS_LINK_DEPENDS) + $(LINK) $(KOKKOS_LDFLAGS) $(LDFLAGS) $(EXTRA_PATH) $(OBJ) $(KOKKOS_LIBS) $(LIB) -o $(EXE) + +clean: kokkos-clean + rm -f *.o *.cuda *.host *.rocm + +# Compilation rules + +%.o:%.cpp $(KOKKOS_CPP_DEPENDS) + $(CXX) $(KOKKOS_CPPFLAGS) $(KOKKOS_CXXFLAGS) $(CXXFLAGS) $(EXTRA_INC) -c $< -o $(notdir $@) diff --git a/README.md b/README.md new file mode 100644 index 000000000..192b2c2d0 --- /dev/null +++ b/README.md @@ -0,0 +1,106 @@ +Spiner +=== + +Performance portable utilities for representing and interpolating +tabulated data. Named for [Brent +Spiner](https://en.wikipedia.org/wiki/Brent_Spiner). + +## Installation + +`Spiner` is self-contained. Simply clone it as +```bash +git clone git@gitlab.lanl.gov:jonahm/spiner.git +``` +If you want to use unit tests, clone with submodules to include `Catch2`. +```bash +git clone --recurse-submodules git@gitlab.lanl.gov:jonahm/spiner.git +``` +To build the relevant object file, call +```bash +make databox +``` +To run unit tests, +```bash +make test +``` +and to do convergence testing, +```bash +make convergence +``` +At the moment, `Spiner` cannot be installed into a system directory. + +**Note that you may have to edit the `Makefile` to set paths to, e.g., `hdf5`.** + +## Dependencies + +`Spiner` has no dependencies for the `databox` tool. Simply include it in your project. It is header-only and requires only a few files: + +- `databox.hpp` +- `interpolation.hpp` +- `spiner_types.hpp` +- `sp5.hpp` + +The testing tooling requires a few different pieces: + +- Unit testing requires [Catch2](https://github.com/catchorg/Catch2), + which is header only. This is included via a git submodule. +- Convergence testing requires the scientific python stack, including: + - python3 + - numpy + - matplotlib + +## HDF5 + +`Spiner` supports reading and writing DataBox objects into a custom HDF5 format called `SP5`. +To enable this, compile with the appropriate `HDF5` linking and the flag `-DSPINER_USE_HDF5`. + +## Features + +- Spiner supports interpolation in arbitrary dimensions, and it's fast in 3d and fewer. +- Spiner supports interpolation onto "subtables" + +## Interpolation + +Interpolation is linear. Here's an example of interpolation in 3D (2D +slice shown). Convergence is second-order, as expected. + +![convergence plot](figs/convergence.png) + +## Copyright + +© (or copyright) 2019-2021. 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. + +This program is open source under the BSD-3 License. Redistribution +and use in source and binary forms, with or without modification, are +permitted provided that the following conditions are met: +1. Redistributions of source code must retain the above copyright +notice, this list of conditions and the following disclaimer. +2. Redistributions in binary form must reproduce the above copyright +notice, this list of conditions and the following disclaimer in the +documentation and/or other materials provided with the distribution. +3. Neither the name of the copyright holder nor the names of its +contributors may be used to endorse or promote products derived from +this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE diff --git a/convergence.cpp b/convergence.cpp new file mode 100644 index 000000000..e230ac278 --- /dev/null +++ b/convergence.cpp @@ -0,0 +1,115 @@ +// © (or copyright) 2019-2021. 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. + +#include +#include +#include +#include + +#include "portability.hpp" +#include "spiner_types.hpp" +#include "databox.hpp" +#include "interpolation.hpp" + +using Spiner::DataBox; +using Spiner::RegularGrid1D; + +const std::string outname = "convergence.dat"; + +constexpr Real KX = 2; +constexpr Real KY = 3; +constexpr Real KZ = 4; + +constexpr Real xmin = 0; +constexpr Real xmax = 1; + +constexpr int NGRIDS = 4; +Real errors[NGRIDS]; +constexpr Real NCOARSE[NGRIDS] = {8,32,128,512}; +constexpr int NFINE = 1024; + +inline Real testFunction(Real z, Real y, Real x) { + return sin(2*M_PI*KX*x)*sin(2*M_PI*KY*y)*sin(2*M_PI*KZ*z); +} + +int main() { + RegularGrid1D gfine(xmin,xmax,NFINE); + + std::cout << "Convergence test for Spiner" << std::endl; + std::cout << std::scientific; + std::cout.precision(15); + + #ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::initialize(); + #endif + { + for (int i = 0; i < NGRIDS; i++) { + int n = NCOARSE[i]; + Real* data = (Real*)PORTABLE_MALLOC(sizeof(Real)*n*n*n); + std::cout << "\tCoarse resolution = " << n << std::endl; + DataBox db(data,n,n,n); + for(int d = 0; d < db.rank(); d++) { + db.setRange(d,xmin,xmax,n); + } + std::cout << "\t\tFilling DataBox" << std::endl; + portableFor("filling databox", + 0,n,0,n,0,n, + PORTABLE_LAMBDA(const int iz, + const int iy, + const int ix) { + Real z = db.range(2).x(iz); + Real y = db.range(1).x(iy); + Real x = db.range(0).x(ix); + db(iz,iy,ix) = testFunction(z,y,x); + }); + std::cout << "\t\tInterpolating..." << std::endl; + errors[i] = 0; + portableReduce("interpolating", + 0,NFINE,0,NFINE,0,NFINE, + PORTABLE_LAMBDA(const int iz, + const int iy, + const int ix, + Real& reduced) { + Real z = gfine.x(iz); + Real y = gfine.x(iy); + Real x = gfine.x(ix); + Real f_interp = db.interpToReal(z,y,x); + Real f_true = testFunction(z,y,x); + Real difference = f_interp - f_true; + reduced += difference*difference; + }, + errors[i]); + errors[i] = sqrt(fabs(errors[i])) / (NFINE*NFINE*NFINE); + std::cout << "\t\t...error = " << errors[i] << std::endl; + PORTABLE_FREE(data); + } + } + #ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::finalize(); + #endif + + std::cout << "Saving convergence results to: " << outname << std::endl; + std::ofstream outfile; + outfile.open(outname); + outfile.precision(15); + outfile << std::scientific; + outfile << "# [0]:resolution, [1]:error" << std::endl; + for (int i = 0; i < NGRIDS; i++) { + outfile << NCOARSE[i] << "\t" << errors[i] << std::endl; + } + outfile.close(); + std::cout << "Done!" << std::endl; + + return 0; +} + diff --git a/databox.hpp b/databox.hpp new file mode 100644 index 000000000..0d4a18678 --- /dev/null +++ b/databox.hpp @@ -0,0 +1,912 @@ +#ifndef _SPINER_DATABOX_HPP_ +#define _SPINER_DATABOX_HPP_ +//====================================================================== +// © (or copyright) 2019-2021. 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. +//====================================================================== + +#include +#include +#include +#include +#include +#include +#include + +#ifdef SPINER_USE_HDF +#include "hdf5.h" +#include "hdf5_hl.h" +#endif + +#include "ports-of-call/portability.hpp" +#include "ports-of-call/portable_arrays.hpp" +#include "spiner_types.hpp" +#include "interpolation.hpp" +#include "sp5.hpp" + +// TODO: get named indices working +// TODO: If asserts are too slow, remove them. + +namespace Spiner { + + enum class IndexType { Interpolated = 0, Named = 1, Indexed = 2 }; + enum class DataStatus {empty, shallow_slice, allocated_hostonly}; + constexpr int MAXRANK = PortableMDArray::MAXDIM; + + class DataBox { + public: + + // Base constructor + PORTABLE_INLINE_FUNCTION __attribute__((nothrow)) DataBox() + : rank_(0) + , status_(DataStatus::empty) + , data_(nullptr) + {} + + // Rank constructors w/ pointer + PORTABLE_INLINE_FUNCTION __attribute__((nothrow)) + DataBox(Real* data, int nx1) + : rank_(1) + , status_(DataStatus::shallow_slice) + , data_(data) + { + dataView_.NewPortableMDArray(data, nx1); + setAllIndexed_(); + } + PORTABLE_INLINE_FUNCTION __attribute__((nothrow)) + DataBox(Real* data, int nx2, int nx1) + : rank_(2) + , status_(DataStatus::shallow_slice) + , data_(data) + { + dataView_.NewPortableMDArray(data, nx2,nx1); + setAllIndexed_(); + } + PORTABLE_INLINE_FUNCTION __attribute__((nothrow)) + DataBox(Real* data, + int nx3, int nx2, int nx1) + : rank_(3) + , status_(DataStatus::shallow_slice) + , data_(data) + { + dataView_.NewPortableMDArray(data, nx3,nx2,nx1); + setAllIndexed_(); + } + PORTABLE_INLINE_FUNCTION __attribute__((nothrow)) + DataBox(Real* data, + int nx4, int nx3, + int nx2, int nx1) + : rank_(4) + , status_(DataStatus::shallow_slice) + , data_(data) + { + dataView_.NewPortableMDArray(data, nx4,nx3,nx2,nx1); + setAllIndexed_(); + } + PORTABLE_INLINE_FUNCTION __attribute__((nothrow)) + DataBox(Real* data, + int nx5, int nx4, + int nx3, int nx2, int nx1) + : rank_(5) + , status_(DataStatus::shallow_slice) + , data_(data) + { + dataView_.NewPortableMDArray(data, nx5,nx4,nx3,nx2,nx1); + setAllIndexed_(); + } + PORTABLE_INLINE_FUNCTION __attribute__((nothrow)) + DataBox(Real* data, + int nx6, int nx5, int nx4, + int nx3, int nx2, int nx1) + : rank_(6) + , status_(DataStatus::shallow_slice) + , data_(data) + { + dataView_.NewPortableMDArray(data, nx6,nx5,nx4,nx3,nx2,nx1); + setAllIndexed_(); + } + + // Rank constructors w/o pointer + inline __attribute__((nothrow)) DataBox(int nx1) + : rank_(1) + , status_(DataStatus::allocated_hostonly) + , data_((Real*)malloc(sizeof(Real)*nx1)) + { + dataView_.NewPortableMDArray(data_, nx1); + setAllIndexed_(); + } + inline __attribute__((nothrow)) DataBox(int nx2, int nx1) + : rank_(2) + , status_(DataStatus::allocated_hostonly) + , data_((Real*)malloc(sizeof(Real)*nx2*nx1)) + { + dataView_.NewPortableMDArray(data_, nx2, nx1); + setAllIndexed_(); + } + inline __attribute__((nothrow)) DataBox(int nx3, int nx2, int nx1) + : rank_(3) + , status_(DataStatus::allocated_hostonly) + , data_((Real*)malloc(sizeof(Real)*nx3*nx2*nx1)) + { + dataView_.NewPortableMDArray(data_, nx3, nx2, nx1); + setAllIndexed_(); + } + inline __attribute__((nothrow)) DataBox(int nx4, int nx3, int nx2, int nx1) + : rank_(4) + , status_(DataStatus::allocated_hostonly) + , data_((Real*)malloc(sizeof(Real)*nx4*nx3*nx2*nx1)) + { + dataView_.NewPortableMDArray(data_, nx4, nx3, nx2, nx1); + setAllIndexed_(); + } + inline __attribute__((nothrow)) + DataBox(int nx5, int nx4, + int nx3, int nx2, int nx1) + : rank_(5) + , status_(DataStatus::allocated_hostonly) + , data_((Real*)malloc(sizeof(Real)*nx5*nx4*nx3*nx2*nx1)) + { + dataView_.NewPortableMDArray(data_, nx5, nx4, nx3, nx2, nx1); + setAllIndexed_(); + } + inline __attribute__((nothrow)) + DataBox(int nx6, int nx5, int nx4, + int nx3, int nx2, int nx1) + : rank_(6) + , status_(DataStatus::allocated_hostonly) + , data_((Real*)malloc(sizeof(Real)*nx6*nx5*nx4*nx3*nx2*nx1)) + { + dataView_.NewPortableMDArray(data_, nx6, nx5, nx4, nx3, nx2, nx1); + setAllIndexed_(); + } + + // Copy and move constructors. ALL SHALLOW + inline __attribute__((nothrow)) DataBox(PortableMDArray A) + : rank_(A.GetRank()) + , status_(DataStatus::shallow_slice) + , data_(A.data()) + , dataView_(A) + { + setAllIndexed_(); + } + inline __attribute__((nothrow)) DataBox(PortableMDArray& A) + : rank_(A.GetRank()) + , status_(DataStatus::shallow_slice) + , data_(A.data()) + , dataView_(A) + { + setAllIndexed_(); + } + PORTABLE_INLINE_FUNCTION __attribute__((nothrow)) + DataBox(const DataBox& src) + : rank_(src.rank_) + , status_(DataStatus::shallow_slice) + , data_(src.data_) + { + setAllIndexed_(); + dataView_.InitWithShallowSlice(src.dataView_,6,0,src.dim(6)); + for (int i = 0; i < rank_; i++) { + indices_[i] = src.indices_[i]; + grids_[i] = src.grids_[i]; + } + } + + // Destructor + PORTABLE_INLINE_FUNCTION __attribute__((nothrow)) ~DataBox() { free_(); } + + // Slice constructor + PORTABLE_INLINE_FUNCTION __attribute__((nothrow)) + DataBox(const DataBox& b, + const int dim, + const int indx, + const int nvar); + +#ifdef SPINER_USE_HDF + // HDF5 constructors + inline DataBox(const std::string& filename) { loadHDF(filename); } + inline DataBox(hid_t loc, const std::string& groupname) { + loadHDF(loc, groupname); + } +#endif // SPINER_USE_HDF + + // Read an array, shallow + inline void setArray(PortableMDArray& A); + + // Equivalent to NewPortableMDArray. + // Note that it's destructive! + inline void resize(int nx1); + inline void resize(int nx2, int nx1); + inline void resize(int nx3, int nx2, int nx1); + inline void resize(int nx4, int nx3, int nx2, int nx1); + inline void resize(int nx5, int nx4, int nx3, int nx2, int nx1); + inline void resize(int nx6, int nx5, int nx4, int nx3, int nx2, int nx1); + + // Index operators + // By reference + PORTABLE_INLINE_FUNCTION Real &operator() (const int n) { + return dataView_(n); + } + PORTABLE_INLINE_FUNCTION Real &operator() (const int n2, const int n1) { + return dataView_(n2,n1); + } + PORTABLE_INLINE_FUNCTION Real &operator() (const int n3, const int n2, + const int n1) { + return dataView_(n3,n2,n1); + } + PORTABLE_INLINE_FUNCTION Real &operator() (const int n4, const int n3, + const int n2, const int n1) { + return dataView_(n4,n3,n2,n1); + } + PORTABLE_INLINE_FUNCTION Real &operator() (const int n5, const int n4, + const int n3, const int n2, + const int n1) { + return dataView_(n5,n4,n3,n2,n1); + } + PORTABLE_INLINE_FUNCTION Real &operator() (const int n6, const int n5, + const int n4, const int n3, + const int n2, const int n1) { + return dataView_(n6,n5,n4,n3,n2,n1); + } + + // By value + PORTABLE_INLINE_FUNCTION Real &operator() (const int n) const { + return dataView_(n); + } + PORTABLE_INLINE_FUNCTION Real &operator() (const int n2, + const int n1) const { + return dataView_(n2,n1); + } + PORTABLE_INLINE_FUNCTION Real &operator() (const int n3, + const int n2, + const int n1) const { + return dataView_(n3,n2,n1); + } + PORTABLE_INLINE_FUNCTION Real &operator() (const int n4, + const int n3, + const int n2, + const int n1) const { + return dataView_(n4,n3,n2,n1); + } + PORTABLE_INLINE_FUNCTION Real &operator() (const int n5, + const int n4 , + const int n3, + const int n2, + const int n1) const { + return dataView_(n5,n4,n3,n2,n1); + } + PORTABLE_INLINE_FUNCTION Real &operator() (const int n6, + const int n5, + const int n4, + const int n3, + const int n2, + const int n1) const { + return dataView_(n6,n5,n4,n3,n2,n1); + } + + // Slice operation + PORTABLE_INLINE_FUNCTION + DataBox slice(const int dim, + const int indx, + const int nvar) const { + return DataBox(*this, dim, indx, nvar); + } + PORTABLE_INLINE_FUNCTION DataBox slice(const int indx) const { + return slice(rank_,indx,1); + } + PORTABLE_INLINE_FUNCTION DataBox slice(const int ix2, const int ix1) const { + //DataBox a(*this, rank_, ix2, 1); + //return DataBox(a, a.rank_, ix1, 1); + return slice(ix2).slice(ix1); + } + + // Reshape the view of the array + PORTABLE_INLINE_FUNCTION void reshape(int nx6, int nx5, int nx4, + int nx3, int nx2, int nx1) { + dataView_.Reshape(nx6,nx5,nx4,nx3,nx2,nx1); + } + PORTABLE_INLINE_FUNCTION void reshape(int nx5, int nx4, int nx3, + int nx2, int nx1) { + dataView_.Reshape(nx5,nx4,nx3,nx2,nx1); + } + PORTABLE_INLINE_FUNCTION void reshape(int nx4, int nx3, int nx2, int nx1) { + dataView_.Reshape(nx4,nx3,nx2,nx1); + } + PORTABLE_INLINE_FUNCTION void reshape(int nx3, int nx2, int nx1) { + dataView_.Reshape(nx3,nx2,nx1); + } + PORTABLE_INLINE_FUNCTION void reshape(int nx2, int nx1) { + dataView_.Reshape(nx2,nx1); + } + PORTABLE_INLINE_FUNCTION void reshape(int nx1) { + dataView_.Reshape(nx1); + } + + + // Interpolates whole DataBox to a real number, + // x1 is fastest index. xN is slowest. + PORTABLE_INLINE_FUNCTION Real + __attribute__((nothrow)) __attribute__((always_inline)) + interpToReal(const Real x) const; + PORTABLE_INLINE_FUNCTION Real + __attribute__((nothrow)) __attribute__((always_inline)) + interpToReal(const Real x2, + const Real x1) const; + PORTABLE_INLINE_FUNCTION Real + __attribute__((nothrow)) __attribute__((always_inline)) + interpToReal(const Real x3, + const Real x2, + const Real x1) const; + // Interpolates SLOWEST indices of databox to a new + // DataBox, interpolated at that slowest index. + // WARNING: requires memory to be pre-allocated. + // TODO: add 3d and higher interpFromDB if necessary + PORTABLE_INLINE_FUNCTION void interpFromDB(const DataBox& db, const Real x); + PORTABLE_INLINE_FUNCTION void interpFromDB(const DataBox& db, + const Real x2, const Real x1); + + // Setters + // NOTE: i ranges from 0 to N-1, where 0 is the FASTEST moving + // index, and N-1 is the slowest + // TODO: is this intuitive? + inline void setIndexType(int i, IndexType t) { + assert( 0 <= i && i < rank_ ); + indices_[i] = t; + } + inline void setRange(int i, RegularGrid1D g) { + setIndexType(i, IndexType::Interpolated); + grids_[i] = g; + } + inline void setRange(int i, Real xmin, Real xmax, int N) { + setRange(i,RegularGrid1D(xmin,xmax,N)); + } + + // Reshapes from other databox, but does not allocate memory. + // Does no checks that memory is available. + // Optionally copies shape of source with ndims fewer slowest-moving dimensions + PORTABLE_INLINE_FUNCTION void copyShape(const DataBox& db, const int ndims=0); + // Returns new databox with same memory and metadata + inline void copyMetadata(const DataBox& src); + +#ifdef SPINER_USE_HDF + inline herr_t saveHDF() const { return saveHDF(SP5::DB::FILENAME); } + inline herr_t saveHDF(const std::string& filename) const; + inline herr_t saveHDF(hid_t loc, const std::string& groupname) const; + inline herr_t loadHDF() { return loadHDF(SP5::DB::FILENAME); } + inline herr_t loadHDF(const std::string& filename); + inline herr_t loadHDF(hid_t loc, const std::string& groupname); +#endif + + // Reference accessors + inline IndexType& indexType(const int i) { return indices_[i]; } + inline RegularGrid1D& range(const int i) { return grids_[i]; } + + // Assignment and move, both perform shallow copy + PORTABLE_INLINE_FUNCTION DataBox& operator= (const DataBox& other); + inline void copy(const DataBox& src); + + // utility info + inline DataStatus dataStatus() { return status_; } + inline bool isReference() { return status_ == DataStatus::shallow_slice; } + inline bool ownsAllocatedMemory() { + return (status_ == DataStatus::empty + || status_ == DataStatus::allocated_hostonly); + } + inline bool operator== (const DataBox& other) const; + inline bool operator!= (const DataBox& other) const { return !(*this == other); } + + PORTABLE_INLINE_FUNCTION Real* data() const { return data_; } + PORTABLE_INLINE_FUNCTION Real min() const; + PORTABLE_INLINE_FUNCTION Real max() const; + PORTABLE_INLINE_FUNCTION int rank() const { return rank_; } + PORTABLE_INLINE_FUNCTION int size() const { return dataView_.GetSize(); } + PORTABLE_INLINE_FUNCTION int sizeBytes() const { + return dataView_.GetSizeInBytes(); + } + PORTABLE_INLINE_FUNCTION int dim(int i) const { return dataView_.GetDim(i); } + PORTABLE_INLINE_FUNCTION void range(int i, + Real& min, Real& max, + Real& dx, int& N) const; + PORTABLE_INLINE_FUNCTION RegularGrid1D range(int i) const; + PORTABLE_INLINE_FUNCTION IndexType indexType(const int i) const { + return indices_[i]; + } + + private: + int rank_; // after dataView_ b/c dataView_ should be initialized first + DataStatus status_; + Real* data_; // sometimes we manage this and sometimes we don't + PortableMDArray dataView_; // always used + IndexType indices_[MAXRANK]; + RegularGrid1D grids_[MAXRANK]; + + PORTABLE_INLINE_FUNCTION void setAllIndexed_(); + PORTABLE_INLINE_FUNCTION bool canInterpToReal_(const int interpOrder) const; + inline std::string gridname_(int i) const { + return SP5::DB::GRID_FORMAT[0] + std::to_string(i+1) + SP5::DB::GRID_FORMAT[1]; + } + PORTABLE_INLINE_FUNCTION __attribute__((nothrow)) + void free_() { + if (status_ == DataStatus::allocated_hostonly) free(data_); + } + }; + + // Slice constructor + PORTABLE_INLINE_FUNCTION __attribute__((nothrow)) + DataBox::DataBox(const DataBox& b, + const int dim, + const int indx, + const int nvar) { + data_ = b.data_; + status_ = DataStatus::shallow_slice; + dataView_.InitWithShallowSlice(b.dataView_,dim,indx,nvar); + rank_ = dataView_.GetRank(); + for (int i = 0; i < rank_; i++) { + indices_[i] = b.indices_[i]; + grids_[i] = b.grids_[i]; + } + } + + // Read an array, shallow + inline void DataBox::setArray(PortableMDArray& A) { + dataView_ = A; + rank_ = A.GetRank(); + setAllIndexed_(); + } + + // Allocate memory of constructed DataBox + inline void DataBox::resize(int nx1) { + assert( ownsAllocatedMemory() ) ; + free_(); + rank_ = 1; + data_ = (Real*)malloc(sizeof(Real)*nx1); + setAllIndexed_(); + dataView_.NewPortableMDArray(data_, nx1); + } + inline void DataBox::resize(int nx2, int nx1) { + assert( ownsAllocatedMemory() ); + free_(); + rank_ = 2; + data_ = (Real*)malloc(sizeof(Real)*nx2*nx1); + setAllIndexed_(); + dataView_.NewPortableMDArray(data_, nx2, nx1); + } + inline void DataBox::resize(int nx3, int nx2, int nx1) { + assert( ownsAllocatedMemory() ); + free_(); + rank_ = 3; + data_ = (Real*)malloc(sizeof(Real)*nx3*nx2*nx1); + setAllIndexed_(); + dataView_.NewPortableMDArray(data_, nx3, nx2, nx1); + } + inline void DataBox::resize(int nx4, int nx3, int nx2, int nx1) { + assert( ownsAllocatedMemory() ); + free_(); + rank_ = 4; + data_ = (Real*)malloc(sizeof(Real)*nx4*nx3*nx2*nx1); + setAllIndexed_(); + dataView_.NewPortableMDArray(data_, nx4, nx3, nx2, nx1); + } + inline void DataBox::resize(int nx5, int nx4, int nx3, int nx2, int nx1) { + assert( ownsAllocatedMemory() ); + free_(); + rank_ = 5; + data_ = (Real*)malloc(sizeof(Real)*nx5*nx4*nx3*nx2*nx1); + setAllIndexed_(); + dataView_.NewPortableMDArray(data_, nx5, nx4, nx3, nx2, nx1); + } + inline void DataBox::resize(int nx6, int nx5, int nx4, + int nx3, int nx2, int nx1) { + assert( ownsAllocatedMemory() ); + free_(); + rank_ = 6; + data_ = (Real*)malloc(sizeof(Real)*nx6*nx5*nx4*nx3*nx2*nx1); + setAllIndexed_(); + dataView_.NewPortableMDArray(data_, nx6, nx5, nx4, nx3, nx2, nx1); + } + + PORTABLE_INLINE_FUNCTION Real DataBox::interpToReal(const Real x) const { + assert( canInterpToReal_(1) ); + return grids_[0](x,dataView_); + } + + PORTABLE_INLINE_FUNCTION Real + __attribute__((nothrow)) __attribute__((always_inline)) + DataBox::interpToReal(const Real x2, + const Real x1) const { + assert( canInterpToReal_(2) ); + int ix1, ix2; + weights_t w1, w2; + grids_[0].weights(x1,ix1,w1); + grids_[1].weights(x2,ix2,w2); + // TODO: prefectch corners for speed? + // TODO: re-order access pattern? + return (w2[0]*(w1[0]*dataView_(ix2,ix1) + +w1[1]*dataView_(ix2,ix1+1)) + + w2[1]*(w1[0]*dataView_(ix2+1,ix1) + +w1[1]*dataView_(ix2+1,ix1+1))); + /* + return ( w2[0]*w1[0]*dataView_(ix2, ix1 ) + + w2[0]*w1[1]*dataView_(ix2, ix1+1) + + w2[1]*w1[0]*dataView_(ix2+1, ix1 ) + + w2[1]*w1[1]*dataView_(ix2+1, ix1+1)); + */ + } + + PORTABLE_INLINE_FUNCTION Real + __attribute__((nothrow)) __attribute__((always_inline)) + DataBox::interpToReal(const Real x3, + const Real x2, + const Real x1) const { + assert( canInterpToReal_(3) ); + int ix[3]; + weights_t w[3]; + grids_[0].weights(x1,ix[0],w[0]); + grids_[1].weights(x2,ix[1],w[1]); + grids_[2].weights(x3,ix[2],w[2]); + // TODO: prefect corners for speed? + // TODO: re-order access pattern? + return (w[2][0]*(w[1][0]*(w[0][0]*dataView_(ix[2],ix[1],ix[0]) + +w[0][1]*dataView_(ix[2],ix[1],ix[0]+1)) + +(w[1][1]*(w[0][0]*dataView_(ix[2],ix[1]+1,ix[0]) + +w[0][1]*dataView_(ix[2],ix[1]+1,ix[0]+1)))) + + w[2][1]*(w[1][0]*(w[0][0]*dataView_(ix[2]+1,ix[1],ix[0]) + +w[0][1]*dataView_(ix[2]+1,ix[1],ix[0]+1)) + + w[1][1]*(w[0][0]*dataView_(ix[2]+1,ix[1]+1,ix[0]) + +w[0][1]*dataView_(ix[2]+1,ix[1]+1,ix[0]+1)))); + /* + return ( w[2][0]*w[1][0]*w[0][0]*dataView_(ix[2], ix[1], ix[0] ) + + w[2][0]*w[1][0]*w[0][1]*dataView_(ix[2], ix[1], ix[0]+1) + + w[2][0]*w[1][1]*w[0][0]*dataView_(ix[2], ix[1]+1, ix[0] ) + + w[2][0]*w[1][1]*w[0][1]*dataView_(ix[2], ix[1]+1, ix[0]+1) + + w[2][1]*w[1][0]*w[0][0]*dataView_(ix[2]+1, ix[1], ix[0] ) + + w[2][1]*w[1][0]*w[0][1]*dataView_(ix[2]+1, ix[1], ix[0]+1) + + w[2][1]*w[1][1]*w[0][0]*dataView_(ix[2]+1, ix[1]+1, ix[0] ) + + w[2][1]*w[1][1]*w[0][1]*dataView_(ix[2]+1, ix[1]+1, ix[0]+1)); + */ + } + + PORTABLE_INLINE_FUNCTION void DataBox::interpFromDB(const DataBox& db, + const Real x) { + assert( db.indices_[db.rank_-1] == IndexType::Interpolated ); + assert( db.grids_[db.rank_-1].isWellFormed() ); + assert( size() == (db.size() / db.dim(db.rank_)) ); + + int ix; + weights_t w; + copyShape(db,1); + + db.grids_[db.rank_-1].weights(x,ix,w); + DataBox lower(db.slice(ix)), upper(db.slice(ix+1)); + //lower = db.slice(ix); + //upper = db.slice(ix+1); + for(int i = 0; i < size(); i++) { + dataView_(i) = w[0]*lower(i) + w[1]*upper(i); + } + } + + PORTABLE_INLINE_FUNCTION void DataBox::interpFromDB(const DataBox& db, + const Real x2, + const Real x1) { + assert( db.rank_ >= 2 ); + assert( db.indices_[db.rank_-1] == IndexType::Interpolated ); + assert( db.grids_[db.rank_-1].isWellFormed() ); + assert( db.indices_[db.rank_-2] == IndexType::Interpolated ); + assert( db.grids_[db.rank_-2].isWellFormed() ); + assert( size() == (db.size() / (db.dim(db.rank_)*db.dim(db.rank_-1))) ); + + int ix2,ix1; + weights_t w2,w1; + copyShape(db,2); + + db.grids_[db.rank_-2].weights(x1, ix1, w1); + db.grids_[db.rank_-1].weights(x2, ix2, w2); + DataBox corners[2][2]{{db.slice(ix2, ix1 ), db.slice(ix2+1, ix1 )}, + {db.slice(ix2, ix1+1 ), db.slice(ix2+1, ix1+1 )}}; + // copyShape(db,2); + // + // db.grids_[db.rank_-2].weights(x1, ix1, w1); + // db.grids_[db.rank_-1].weights(x2, ix2, w2); + //corners[0][0] = db.slice(ix2, ix1 ); + //corners[1][0] = db.slice(ix2, ix1+1 ); + //corners[0][1] = db.slice(ix2+1, ix1 ); + //corners[1][1] = db.slice(ix2+1, ix1+1 ); + /* + for (int i = 0; i < size(); i++) { + dataView_(i) = ( w2[0]*w1[0]*corners[0][0](i) + + w2[0]*w1[1]*corners[1][0](i) + + w2[1]*w1[0]*corners[0][1](i) + + w2[1]*w1[1]*corners[1][1](i)); + } + */ + for (int i = 0; i < size(); i++) { + dataView_(i) = ( w2[0]*(w1[0]*corners[0][0](i) + + w1[1]*corners[1][0](i)) + + w2[1]*(w1[0]*corners[0][1](i) + + w1[1]*corners[1][1](i))); + } + } + + // Reshapes from other databox, but does not allocate memory. + // Does no checks that memory is available. + // Optionally copies shape of source with ndims fewer slowest-moving dimensions + PORTABLE_INLINE_FUNCTION void DataBox::copyShape(const DataBox& db, + const int ndims) { + rank_ = db.rank_ - ndims; + int dims[MAXRANK]; + for (int i = 0; i < MAXRANK; i++) dims[i] = 1; + setAllIndexed_(); // TODO: can remove + for (int i = rank_ - 1; i >= 0; i--) dims[i] = db.dim(i+1); + reshape(dims[5],dims[4],dims[3],dims[2],dims[1],dims[0]); + rank_ = db.rank_ - ndims; + for (int i = 0; i < rank_; i++) { + indices_[i] = db.indices_[i]; + grids_[i] = db.grids_[i]; + } + } + // reallocates and then copies shape from other databox + // everything but the actual copy in a deep copy + inline void DataBox::copyMetadata(const DataBox& src) { + resize(src.dim(6),src.dim(5),src.dim(4), + src.dim(3),src.dim(2),src.dim(1)); + rank_ = src.rank_; // resize sets rank + for (int i = 0; i < rank_; i++) { + grids_[i] = src.grids_[i]; + indices_[i] = src.indices_[i]; + } + } + +#ifdef SPINER_USE_HDF + inline herr_t DataBox::saveHDF(const std::string& filename) const { + herr_t status; + hid_t file; + + file = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, + H5P_DEFAULT, H5P_DEFAULT); + status = saveHDF(file, SP5::DB::GRPNAME); + status += H5Fclose(file); + return status; + } + + inline herr_t DataBox::saveHDF(hid_t loc, const std::string& groupname) const { + hid_t group, grids; + herr_t status = 0; + + std::vector dims_int(rank_); + for (int i = 0; i < rank_; i++) dims_int[i] = dim(i+1); + + std::vector dims_hsize_t(rank_); + for (int i = 0; i < rank_; i++) dims_hsize_t[i] = dim(rank_ - i); + + std::vector index_types(rank_); + for (int i = 0; i < rank_; i++) { + index_types[i] = static_cast(indices_[i]); + } + + // Greate group + group = H5Gcreate(loc, groupname.c_str(), + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + + // Record rank as an attribute + status += H5LTset_attribute_int(loc, groupname.c_str(), + SP5::DB::RANKNAME, + &rank_, 1); + + // Index types + status += H5LTset_attribute_int(loc, groupname.c_str(), + SP5::DB::IDXSNAME, + index_types.data(), rank_); + status += H5LTset_attribute_string(loc, groupname.c_str(), + SP5::DB::IDXINFONAME, + SP5::DB::IDXINFO); + + // Dimensions of the PortableMDArray, set as an attribute + status += H5LTset_attribute_int(loc, groupname.c_str(), + SP5::DB::DIMSNAME, + dims_int.data(), + rank_); + + // Save the PortableMDArray, in a human-readable shape + status += H5LTmake_dataset(group, SP5::DB::DSETNAME, + rank_, dims_hsize_t.data(), + H5T_REAL, dataView_.data()); + + // Grids + grids = H5Gcreate(group, SP5::DB::GRIDNAME, + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + for (int i = 0; i < rank_; i++) { + if (indices_[i] == IndexType::Interpolated) { + + status += grids_[i].saveHDF(grids, gridname_(i).c_str()); + } + } + status += H5Gclose(grids); + status += H5Gclose(group); + return status; + } + + inline herr_t DataBox::loadHDF(const std::string& filename) { + herr_t status; + hid_t file = H5Fopen(filename.c_str(), + H5F_ACC_RDONLY, H5P_DEFAULT); + status = loadHDF(file, SP5::DB::GRPNAME); + return status; + } + + inline herr_t DataBox::loadHDF(hid_t loc, const std::string& groupname) { + hid_t group, grids; + herr_t status = 0; + std::vector index_types; + std::vector dims(6,1); + + // Open group + group = H5Gopen(loc, groupname.c_str(), H5P_DEFAULT); + // Get rank + status += H5LTget_attribute_int(loc, groupname.c_str(), + SP5::DB::RANKNAME, + &rank_); + // Resize metadata fields + setAllIndexed_(); + + // Get dimensions + dims.resize(rank_); + status += H5LTget_attribute_int(loc, groupname.c_str(), + SP5::DB::DIMSNAME, + dims.data()); + + // Allocate PortableMDArray and read it in to buffer + free_(); + data_ = (Real*)malloc(sizeof(Real) + *dims[5]*dims[4]*dims[3] + *dims[2]*dims[1]*dims[0]); + dataView_.NewPortableMDArray(data_, + dims[5], dims[4], dims[3], + dims[2], dims[1], dims[0]); + status += H5LTread_dataset(group, SP5::DB::DSETNAME, + H5T_REAL, dataView_.data()); + + // Get index types + index_types.resize(rank_); + status += H5LTget_attribute_int(loc, groupname.c_str(), + SP5::DB::IDXSNAME, + index_types.data()); + for (int i = 0; i < rank_; i++) { + indices_[i] = static_cast(index_types[i]); + } + + // Get grids + grids = H5Gopen(group, SP5::DB::GRIDNAME, H5P_DEFAULT); + for (int i = 0; i < rank_; i++) { + if (indices_[i] == IndexType::Interpolated) { + status += grids_[i].loadHDF(grids, gridname_(i).c_str()); + } + } + status += H5Gclose(grids); + + status += H5Gclose(group); + return status; + } +#endif // SPINER_USE_HDF + + // Performs shallow copy by default + PORTABLE_INLINE_FUNCTION DataBox& DataBox::operator= (const DataBox& src) { + if (this != &src) { + rank_ = src.rank_; + free_(); + status_ = DataStatus::shallow_slice; + data_ = src.data_; + dataView_.InitWithShallowSlice(src.dataView_,6,0,src.dim(6)); + for (int i = 0; i < rank_; i++) { + indices_[i] = src.indices_[i]; + grids_[i] = src.grids_[i]; + } + } + return *this; + } + + // Performs a deep copy + inline void DataBox::copy(const DataBox& src) { + copyMetadata(src); + for (int i = 0; i < src.size(); i++) dataView_(i) = src(i); + } + + inline bool DataBox::operator== (const DataBox& other) const { + if (rank_ != other.rank_) return false; + for (int i = 0; i < rank_; i++) { + if (indices_[i] != other.indices_[i]) return false; + if (indices_[i] == IndexType::Interpolated + && grids_[i] != other.grids_[i]) { + return false; + } + } + return dataView_ == other.dataView_; + } + + // TODO: should this be std::reduce? + PORTABLE_INLINE_FUNCTION Real DataBox::min() const { + Real min = std::numeric_limits::infinity(); + for (int i = 0; i < size(); i++) { + min = std::min(min,dataView_(i)); + } + return min; + } + + PORTABLE_INLINE_FUNCTION Real DataBox::max() const { + Real max = -std::numeric_limits::infinity(); + for (int i = 0; i < size(); i++) { + max = std::max(max,dataView_(i)); + } + return max; + } + + PORTABLE_INLINE_FUNCTION void DataBox::range(int i, + Real& min, Real& max, + Real& dx, int& N) const { + assert( 0 <= i && i < rank_ ); + assert( indices_[i] == IndexType::Interpolated ); + min = grids_[i].min(); + max = grids_[i].max(); + dx = grids_[i].dx(); + N = grids_[i].nPoints(); + } + + PORTABLE_INLINE_FUNCTION RegularGrid1D DataBox::range(int i) const { + assert( 0 <= i && i < rank_ ); + assert( indices_[i] == IndexType::Interpolated ); + return grids_[i]; + } + + PORTABLE_INLINE_FUNCTION void DataBox::setAllIndexed_() { + for (int i = 0; i < rank_; i++) { + indices_[i] = IndexType::Indexed; + } + } + + PORTABLE_INLINE_FUNCTION bool DataBox::canInterpToReal_(const int interpOrder) const { + if (rank_ != interpOrder) return false; + for( int i = 0; i < rank_; i++) { + if ( indices_[i] != IndexType::Interpolated ) return false; + if ( !(grids_[i].isWellFormed()) ) return false; + } + return true; + } + + inline DataBox getOnDeviceDataBox(const DataBox& a_host) + { +#ifdef PORTABILITY_STRATEGY_KOKKOS + using HS = Kokkos::HostSpace; + using DMS = Kokkos::DefaultExecutionSpace::memory_space; + constexpr const bool execution_is_host {std::is_same::value}; + if (execution_is_host) { + return a_host; + } else { + using memUnmanaged = Kokkos::MemoryUnmanaged; + using HostView_t = Kokkos::View; + using DeviceView_t = Kokkos::View; + using Kokkos::deep_copy; + Real* device_data = (Real*)PORTABLE_MALLOC(a_host.sizeBytes()); + DeviceView_t devView(device_data, a_host.size()); + HostView_t hostView(a_host.data(), a_host.size()); + deep_copy(devView, hostView); + DataBox a {devView.data(), + a_host.dim(6), a_host.dim(5), a_host.dim(4), + a_host.dim(3), a_host.dim(2), a_host.dim(1)}; + a.copyShape(a_host); + return a; + } +#else // no kokkos + return a_host; +#endif // kokkos + } + +} +#endif // _SPINER_DATABOX_HPP_ diff --git a/figs/convergence.pdf b/figs/convergence.pdf new file mode 100644 index 000000000..03bf44591 Binary files /dev/null and b/figs/convergence.pdf differ diff --git a/figs/convergence.png b/figs/convergence.png new file mode 100644 index 000000000..1d16840bf Binary files /dev/null and b/figs/convergence.png differ diff --git a/interpolation.hpp b/interpolation.hpp new file mode 100644 index 000000000..204fd37a7 --- /dev/null +++ b/interpolation.hpp @@ -0,0 +1,175 @@ +#ifndef _SPINER_INTERP_ +#define _SPINER_INTERP_ +//====================================================================== +// © (or copyright) 2019-2021. 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. +//====================================================================== + +#include +#include +#include +#include + +#ifdef SPINER_USE_HDF +#include "hdf5.h" +#include "hdf5_hl.h" +#include +#endif + +#include "ports-of-call/portability.hpp" +#include "ports-of-call/portable_arrays.hpp" +#include "spiner_types.hpp" +#include "sp5.hpp" + +namespace Spiner { + // TODO: be more careful about what this number should be + // sqrt machine epsilon or something + constexpr Real EPS = 10.0*std::numeric_limits::epsilon(); + constexpr Real rNaN = std::numeric_limits::signaling_NaN(); + constexpr int iNaN = std::numeric_limits::signaling_NaN(); + + // a poor-man's std::double + struct weights_t { + Real first, second; + PORTABLE_INLINE_FUNCTION Real& operator[] (const int i) { + assert ( 0 <= i && i <= 1 ); + return i == 0 ? first : second; + } + }; + + class RegularGrid1D { + public: + // Constructors + PORTABLE_INLINE_FUNCTION RegularGrid1D() + : min_(rNaN), max_(rNaN), dx_(rNaN), idx_(rNaN), N_(iNaN) + {} + PORTABLE_INLINE_FUNCTION RegularGrid1D(Real min, Real max, size_t N) + : min_(min) + , max_(max) + , dx_((max - min)/((Real)(N-1))) + , idx_(1/dx_) + , N_(N) + {} + + // Assignment operator + PORTABLE_INLINE_FUNCTION RegularGrid1D& operator= (const RegularGrid1D &src) { + if (this != &src) { + min_ = src.min_; + max_ = src.max_; + dx_ = src.dx_; + idx_ = src.idx_; + N_ = src.N_; + } + return *this; + } + + // Forces x in the interval + PORTABLE_INLINE_FUNCTION int bound(int ix) const { + #ifndef SPINER_DISABLE_BOUNDS_CHECKS + if (ix < 0) ix = 0; + if (ix >= (int)N_-1) ix = (int)N_ - 2; // Ensures ix+1 exists + #endif + return ix; + } + + // Gets real value at index + PORTABLE_INLINE_FUNCTION Real x(const int i) const { return i*dx_ + min_; } + PORTABLE_INLINE_FUNCTION int index(const Real x) const { + return bound(idx_*(x - min_)); + } + + // Returns closest index and weights for interpolation + PORTABLE_INLINE_FUNCTION void weights(Real x, int& ix, weights_t& w) const { + ix = index(x); + const Real floor = ix*dx_ + min_; + w[1] = idx_*(x - floor); + w[0] = (1. - w[1]); + } + + // 1D interpolation + PORTABLE_INLINE_FUNCTION Real operator() (const Real x, + const PortableMDArray& A) const { + int ix; + weights_t w; + weights(x,ix,w); + return w[0]*A(ix) + w[1]*A(ix+1); + } + + // utitilies + PORTABLE_INLINE_FUNCTION bool operator== (const RegularGrid1D& other) const { + return (min_ == other.min_ + && max_ == other.max_ + && dx_ == other.dx_ + && idx_ == other.idx_ + && N_ == other.N_); + } + PORTABLE_INLINE_FUNCTION bool operator!= (const RegularGrid1D& other) const { + return !(*this == other); + } + PORTABLE_INLINE_FUNCTION Real min() const { return min_; } + PORTABLE_INLINE_FUNCTION Real max() const { return max_; } + PORTABLE_INLINE_FUNCTION Real dx() const { return dx_; } + PORTABLE_INLINE_FUNCTION Real nPoints() const { return N_; } + PORTABLE_INLINE_FUNCTION bool isnan() const { + return (std::isnan(min_) + || std::isnan(max_) + || std::isnan(dx_) + || std::isnan(idx_) + || std::isnan((Real)N_)); + } + PORTABLE_INLINE_FUNCTION bool isWellFormed() const { return !isnan(); } + +#ifdef SPINER_USE_HDF + inline herr_t saveHDF(hid_t loc, const std::string& name) const { + herr_t status; + Real range[] = {min_,max_,dx_}; + hsize_t range_dims[] = {3}; + int n = static_cast(N_); + status = H5LTmake_dataset(loc, name.c_str(), + SP5::RG1D::RANGE_RANK, + range_dims, + H5T_REAL, range); + status += H5LTset_attribute_int(loc, name.c_str(), + SP5::RG1D::N, + &n, 1); + status += H5LTset_attribute_string(loc, name.c_str(), + SP5::RG1D::RANGE_INFONAME, + SP5::RG1D::RANGE_INFO); + return status; + } + + inline herr_t loadHDF(hid_t loc, const std::string& name) { + herr_t status; + Real range[3]; + int n; + status = H5LTread_dataset(loc, name.c_str(), H5T_REAL, range); + min_ = range[0]; + max_ = range[1]; + dx_ = range[2]; + idx_ = 1./dx_; + status += H5LTget_attribute_int(loc, name.c_str(), + SP5::RG1D::N, + &n); + N_ = n; + return status; + } +#endif + + private: + Real min_, max_; + Real dx_, idx_; + size_t N_; + }; + +} +#endif // _SPINER_INTERP_ diff --git a/plot_convergence.py b/plot_convergence.py new file mode 100755 index 000000000..25f394aaa --- /dev/null +++ b/plot_convergence.py @@ -0,0 +1,48 @@ +#!/usr/bin/env python + +# Generates a convergence plot for Spiner. +# Requires scientific python stack + +import numpy as np +import matplotlib as mpl +from matplotlib import pyplot as plt +from matplotlib import rc +import os +rc('font',size=14) +# mpl.rcParams['xtick.minor.size'] = 0 +# mpl.rcParams['xtick.minor.width'] = 0 + +DFILE = "convergence.dat" +PLOTPATH = os.path.join("figs","convergence") + +KX,KY,KZ = 2,3,4 +xmin,xmax = 0,1 +def test_function(y,x): + return np.sin(2*np.pi*KX*x)*np.sin(2*np.pi*KY*y) + +x = np.linspace(0,1,100) +X,Y = np.meshgrid(x,x) + +data = np.loadtxt(DFILE) + +fig,axarr = plt.subplots(1,2,figsize=(17,8)) + +axarr[0].plot(np.log2(data[:,0]),data[:,1],'ko-',lw=2) +axarr[0].set_yscale('log') +axarr[0].set_xticks(np.log2(data[:,0])) +axarr[0].set_xticklabels([str(d) for d in data[:,0]]) +axarr[0].set_xlabel('Points per axis') +axarr[0].set_ylabel('Error') +axarr[0].grid(True) + +mesh = axarr[1].pcolormesh(X,Y,test_function(Y,X)) +mesh.set_edgecolor('face') +cbar = plt.colorbar(mesh) +cbar.set_label(r'$\sin(2\pi k_x x)\sin(2\pi k_y y)$') +axarr[1].set_xlabel(r'$x$') +axarr[1].set_ylabel(r'$y$') + +for postfix in ['.pdf','.png']: + plt.savefig(PLOTPATH+postfix, + bbox_inches='tight') + diff --git a/ports-of-call/portable_arrays.hpp b/ports-of-call/portable_arrays.hpp index 31aff6e39..a548885b2 100644 --- a/ports-of-call/portable_arrays.hpp +++ b/ports-of-call/portable_arrays.hpp @@ -5,7 +5,7 @@ // Copyright(C) 2014 James M. Stone and other code contributors // Licensed under the 3-clause BSD License, see LICENSE file for details // -// © (or copyright) 2020. Triad National Security, LLC. All rights +// © (or copyright) 2019-2021. 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. diff --git a/sp5.hpp b/sp5.hpp new file mode 100644 index 000000000..b58ae7eda --- /dev/null +++ b/sp5.hpp @@ -0,0 +1,46 @@ +#ifndef _SPINER_SP5_HPP_ +#define _SPINER_SP5_HPP_ +//====================================================================== +// © (or copyright) 2019-2021. 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. +//====================================================================== + +// This file contains strings defined for use accross the SP5 data +// format + +namespace SP5 { + + namespace DB { + constexpr char FILENAME[] = "databox.sp5"; + constexpr char GRPNAME[] = "databox"; + constexpr char DSETNAME[] = "data"; + constexpr char RANKNAME[] = "rank"; + constexpr char DIMSNAME[] = "dims"; + constexpr char IDXSNAME[] = "index_types"; + constexpr char IDXINFONAME[] = "index_types_info"; + constexpr char IDXINFO[] = "Interpolated:0\nNamed:1\nIndexed:2"; + constexpr char GRIDNAME[] = "grids"; + const std::string GRID_FORMAT[] = {"grid_[","]"}; + } + + namespace RG1D { + constexpr char RANGE_NAME[] = "range"; + constexpr char N[] = "npoints"; + constexpr char RANGE_INFONAME[] = "range columns"; + constexpr char RANGE_INFO[] = "[0]:min [1]:max [2]:dx"; + constexpr int RANGE_RANK = 1; + } + +} + +#endif // _SPINER_SP5_HPP_ diff --git a/spiner.py b/spiner.py new file mode 100644 index 000000000..35c18c33a --- /dev/null +++ b/spiner.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python + +#====================================================================== +# © (or copyright) 2019-2021. 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. +#====================================================================== + +from enum import IntEnum + +class IndexType(IntEnum): + Interpolated = 0 + Named = 1 + Indexed = 2 + +class DataBox: + def __init__(self, + indices=[], + grids={}, + data=None): + self._indices=indices + self._grids=grids + self._data=data + + def __getitem__(self,key): + return self._data[key] + + def shape(self): + return self._data.shape + + def grids(self,i): + return self._grids[i] + + def indices(self,i): + return self._indices + + def data(self): + return self._data + + @classmethod + def fromHDF(cls,loc): + import h5py + import numpy as np + indices = loc.attrs['index_types'] + grids = {} + for i in range(indices.shape[0]): + if indices[i] == IndexType.Interpolated: + gridname = 'grids/grid_[{}]'.format(i+1) + xmin,xmax,dx = loc[gridname][()] + nx = loc[gridname].attrs['npoints'][0] + grids[i] = np.linspace(xmin,xmax,nx) + data = loc['data'][()] + return cls(indices,grids,data) diff --git a/spiner_types.hpp b/spiner_types.hpp new file mode 100644 index 000000000..43690e582 --- /dev/null +++ b/spiner_types.hpp @@ -0,0 +1,22 @@ +// © (or copyright) 2019-2021. 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. + +#ifndef _SPINER_TYPES_HPP_ +#define _SPINER_TYPES_HPP_ + +#ifdef SPINER_USE_HDF +#define H5T_REAL H5T_NATIVE_DOUBLE +#define H5_SUCCESS 0 +#endif + +#endif // _SPINER_TYPES_HPP_ diff --git a/test.cpp b/test.cpp new file mode 100644 index 000000000..2ff1b04c4 --- /dev/null +++ b/test.cpp @@ -0,0 +1,435 @@ +// © (or copyright) 2019-2021. 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. + +#include +#include +#include // std::min, std::max +#include // sqrt + +#include "ports-of-call/portability.hpp" +#include "ports-of-call/portable_arrays.hpp" +#include "spiner_types.hpp" +#include "databox.hpp" +#include "interpolation.hpp" + +#define CATCH_CONFIG_RUNNER +#include "catch.hpp" + +using Spiner::DataBox; +using Spiner::RegularGrid1D; +using Spiner::IndexType; +const Real EPSTEST = std::sqrt(Spiner::EPS); + +PORTABLE_INLINE_FUNCTION Real linearFunction(Real z, Real y, Real x) { + return x + y + z; +} + +TEST_CASE( "PortableMDArrays can be allocated from a pointer", + "[PortableMDArray]" ) { + constexpr int N = 2; + constexpr int M = 3; + std::vector data(N*M); + PortableMDArray a; + int tot = 0; + for (int i = 0; i < N*M; i++) { + data[i] = tot; + tot++; + } + a.NewPortableMDArray(data.data(),M,N); + + SECTION( "Shape should be NxM" ) { + REQUIRE( a.GetDim1() == N ); + REQUIRE( a.GetDim2() == M ); + } + + SECTION( "Stride is as set by initialized pointer" ) { + int tot = 0; + for (int j = 0; j < M; j++) { + for (int i = 0; i < N; i++) { + REQUIRE( a(j,i) == tot ); + tot++; + } + } + } + + SECTION( "Identical slices of the same data should compare equal" ) { + PortableMDArray aslc1, aslc2; + aslc1.InitWithShallowSlice(a,1,0,2); + aslc2.InitWithShallowSlice(a,1,0,2); + REQUIRE( aslc1 == aslc2 ); + } +} + +TEST_CASE( "DataBox Basics", "[DataBox]" ) { + + SECTION( "DataBoxes are initialized with correct rank" ) { + DataBox db(2); + DataBox db4(5,4,2,2); + REQUIRE( db.rank() == 1 ); + REQUIRE( db4.rank() == 4 ); + } + + SECTION( "A DataBox can be written to and read from" ) { + + constexpr int M = 3; + constexpr int N = 2; + + DataBox db(M,N); + int tot = 0; + for (int j = 0; j < M; j++) { + for (int i = 0; i < N; i++) { + db(j,i) = tot; + tot++; + } + } + tot = 0; + for (int j = 0; j < M; j++) { + for (int i = 0; i < N; i++) { + REQUIRE( db(j,i) == tot ); + tot++; + } + } + + SECTION( "DataBox min and max can be correctly computed" ) { + REQUIRE( db.max() == tot - 1 ); + REQUIRE( db.min() == 0 ); + } + + SECTION( "DataBox metadata can be copied" ) { + DataBox dbCopy; dbCopy.copyMetadata(db); + REQUIRE( dbCopy.rank() == db.rank() ); + for (int i = 0; i < db.rank(); i++ ) { + REQUIRE( dbCopy.dim(i+1) == db.dim(i+1) ); + REQUIRE( dbCopy.indexType(i) == db.indexType(i) ); + } + REQUIRE( dbCopy != db ); + + SECTION( "DataBoxes can be resized" ) { + dbCopy.resize(5,4,3); + REQUIRE( dbCopy.rank() == 3 ); + REQUIRE( dbCopy.dim(1) == 3 ); + REQUIRE( dbCopy.dim(2) == 4 ); + REQUIRE( dbCopy.dim(3) == 5 ); + } + } + + SECTION("DataBoxes can be shallow copied") { + DataBox db2(db); + REQUIRE( &(db2(0)) == &(db(0)) ); + db2 = db; + REQUIRE( &(db2(0)) == &(db(0)) ); + } + + SECTION("DataBoxes can be deep copied") { + DataBox db2; + db2.copy(db); + REQUIRE( &(db2(0)) != &(db(0)) ); + } + + SECTION( "DataBoxes can be sliced in 2D" ) { + DataBox dbslc = db.slice(0); + DataBox dbslc2(db,1,0,2); + + REQUIRE( dbslc2.rank() == 1 ); + REQUIRE( dbslc2.dim(dbslc2.rank()) == 2 ); + + SECTION( "DataBox slices are correctly indexed" ) { + int tot = 0; + for (int i = 0; i < dbslc.dim(dbslc.rank()); i++) { + REQUIRE( dbslc(i) == tot ); + tot++; + } + } + + SECTION( "DataBox slices are shallow" ) { + REQUIRE( dbslc == dbslc2 ); + REQUIRE( &(dbslc(0)) == &(db(0)) ); + } + } + } +} + +TEST_CASE( "DataBox interpolation", "[DataBox]" ) { + constexpr int NFINE = 100; + constexpr int RANK = 3; + constexpr int NZ = 8; + constexpr int NY = 10; + constexpr int NX = 12; + DataBox db(NZ,NY,NX); + + constexpr Real xmin = 0; + constexpr Real xmax = 1; + constexpr Real ymin = -0.5; + constexpr Real ymax = 0.5; + constexpr Real zmin = -1; + constexpr Real zmax = 0; + + std::array grids = {RegularGrid1D(xmin,xmax,NX), + RegularGrid1D(ymin,ymax,NY), + RegularGrid1D(zmin,zmax,NZ)}; + std::array fine_grids = {RegularGrid1D(xmin,xmax,NFINE), + RegularGrid1D(ymin,ymax,NFINE), + RegularGrid1D(zmin,zmax,NFINE)}; + + for (int i = 0; i < RANK; i++) db.setRange(i, grids[i]); + + for (int iz = 0; iz < NZ; iz++) { + Real z = grids[2].x(iz); + for (int iy = 0; iy < NY; iy++) { + Real y = grids[1].x(iy); + for (int ix = 0; ix < NX; ix++) { + Real x = grids[0].x(ix); + db(iz,iy,ix) = linearFunction(z,y,x); + } + } + } + + SECTION( "interpToReal in 3D is exact for linear functions" ) { + Real error = 0; + for (int iz = 0; iz < NFINE; iz++) { + Real z = fine_grids[2].x(iz); + for (int iy = 0; iy < NFINE; iy++) { + Real y = fine_grids[1].x(iy); + for (int ix = 0; ix < NFINE; ix++) { + Real x = fine_grids[0].x(ix); + Real f_true = linearFunction(z,y,x); + Real difference = db.interpToReal(z,y,x) - f_true; + error += (difference*difference); + } + } + } + error = sqrt(error); + REQUIRE( error <= EPSTEST ); + } + + SECTION( "interpFromDB 3D->2D" ) { + constexpr Real z = (zmax + zmin) / 2.; + + SECTION( "Slicing relevant for interpToDB in slowest index works" ) { + int iz = grids[RANK-1].index(z); + DataBox lower = db.slice(iz); + DataBox upper = db.slice(iz+1); + + Real error = 0; + for (int iy = 0; iy < NY; iy++) { + for (int ix = 0; ix < NX; ix++) { + Real difference = lower(iy, ix) - db(iz, iy, ix); + error += difference*difference; + difference = upper(iy, ix) - db(iz+1, iy, ix); + error += difference*difference; + } + } + error = sqrt(0.5*error); + REQUIRE( error <= EPSTEST ); + } + + DataBox db2d; + db2d.resize(db.size()/db.dim(db.rank())); + db2d.interpFromDB(db,z); + + Real error = 0; + for (int iy = 0; iy < NY; iy++) { + Real y = grids[1].x(iy); + for (int ix = 0; ix < NX; ix++) { + Real x = grids[0].x(ix); + Real f_true = linearFunction(z,y,x); + Real difference = db2d(iy,ix) - f_true; + error += (difference*difference); + } + } + error = sqrt(error); + REQUIRE( error <= EPSTEST ); + + SECTION( "interpToReal 2D" ) { + Real error = 0; + for (int iy = 0; iy < NFINE; iy++) { + Real y = fine_grids[1].x(iy); + for (int ix = 0; ix < NFINE; ix++) { + Real x = fine_grids[0].x(ix); + Real f_true = linearFunction(z,y,x); + Real difference = db2d.interpToReal(y,x) - f_true; + error += (difference*difference); + } + } + error = sqrt(error); + REQUIRE( error <= EPSTEST ); + } + } + + SECTION( "interpFromDB 3D->1D" ) { + constexpr Real z = (zmax + zmin) / 2.; + constexpr Real y = (ymax + ymin) / 2.; + + SECTION( "Slicing in 2D works" ) { + int iz = grids[RANK-1].index(z); + int iy = grids[RANK-2].index(y); + DataBox corner = db.slice(iz,iy); + Real error = 0; + for(int ix = 0; ix < NX; ix++) { + error += (corner(ix) - db(iz,iy,ix))*(corner(ix) - db(iz,iy,ix)); + } + error = sqrt(error); + REQUIRE( error <= EPSTEST ); + } + + DataBox db1d; + db1d.resize(db.size() / (db.dim(db.rank())*db.dim(db.rank()-1))); + db1d.interpFromDB(db,z,y); + REQUIRE( db1d.rank() == 1 ); + REQUIRE( db1d.dim(1) == NX ); + + Real error = 0; + for (int ix = 0; ix < NX; ix++) { + Real x = grids[0].x(ix); + Real f_true = linearFunction(z,y,x); + Real difference = db1d(ix) - f_true; + error += difference*difference; + } + error = sqrt(error); + REQUIRE( error <= EPSTEST ); + } +} + +#if SPINER_USE_HDF +SCENARIO( "DataBox HDF5", "[DataBox],[HDF5]" ) { + constexpr int N = 2; + herr_t status; + + DataBox db(N,N,N); + int tot = 0; + for (int k = 0; k < N; k++) { + for (int j = 0; j < N; j++) { + for (int i = 0; i < N; i++) { + db(k,j,i) = tot++; + } + } + } + db.setRange(0,0,1,10); + + GIVEN( "DataBox can be saved to HDF5" ) { + status = db.saveHDF(); + REQUIRE( status == H5_SUCCESS ); + + WHEN( "DataBox can be loaded from HDF5" ) { + DataBox db2; + status = db2.loadHDF(); + REQUIRE( status == H5_SUCCESS ); + + THEN( "Metadata read in is consistent" ) { + REQUIRE( db2.rank() == db.rank() ); + for (int i = 0; i < db.rank(); i++) { + REQUIRE( db.indexType(i) == db2.indexType(i) ); + REQUIRE( db.dim(i+1) == db2.dim(i+1) ); + if ( db.indexType(i) == IndexType::Interpolated ) { + REQUIRE( db.range(i) == db2.range(i) ); + } + } + AND_THEN( "Data itself is consistent" ) { + for (int k = 0; k < N; k++) { + for (int j = 0; j < N; j++) { + for (int i = 0; i < N; i++) { + REQUIRE( db(k,j,i) == db2(k,j,i) ); + } + } + } + } + } + } + } +} +#endif + + +#ifdef PORTABILITY_STRATEGY_KOKKOS +SCENARIO( "Kokkos functionality: interpolation", + "[DataBox],[Kokkos]" ) { + constexpr int NFINE = 100; + constexpr int RANK = 3; + constexpr int NZ = 8; + constexpr int NY = 10; + constexpr int NX = 12; + DataBox db(NZ,NY,NX); + + constexpr Real xmin = 0; + constexpr Real xmax = 1; + constexpr Real ymin = -0.5; + constexpr Real ymax = 0.5; + constexpr Real zmin = -1; + constexpr Real zmax = 0; + + using Policy3D = Kokkos::MDRangePolicy>; + using DeviceView_t = Kokkos::View; + using HostView_t = Kokkos::View; + + std::array grids = {RegularGrid1D(xmin,xmax,NX), + RegularGrid1D(ymin,ymax,NY), + RegularGrid1D(zmin,zmax,NZ)}; + std::arrayfine_grids = {RegularGrid1D(xmin,xmax,NFINE), + RegularGrid1D(ymin,ymax,NFINE), + RegularGrid1D(zmin,zmax,NFINE)}; + + for (int i = 0; i < RANK; i++) db.setRange(i, grids[i]); + + for (int iz = 0; iz < NZ; iz++) { + Real z = grids[2].x(iz); + for (int iy = 0; iy < NY; iy++) { + Real y = grids[1].x(iy); + for (int ix = 0; ix < NX; ix++) { + Real x = grids[0].x(ix); + db(iz,iy,ix) = linearFunction(z,y,x); + } + } + } + + Real* device_data = (Real*)PORTABLE_MALLOC(db.sizeBytes()); + DeviceView_t deviceView(device_data,db.size()); + HostView_t hostView(db.data(),db.size()); + Kokkos::deep_copy(deviceView, hostView); + DataBox db_dev(device_data,NZ,NY,NX); + db_dev.copyShape(db); + + Real error = 0; + Kokkos::parallel_reduce(Policy3D({0,0,0},{NFINE,NFINE,NFINE}), + PORTABLE_LAMBDA(const int iz, const int iy, const int ix, Real& update) + { + const Real z = fine_grids[2].x(iz); + const Real y = fine_grids[1].x(iy); + const Real x = fine_grids[0].x(ix); + const Real f_true = linearFunction(z,y,x); + const Real difference = db_dev.interpToReal(z,y,x) - f_true; + update += difference*difference; + }, error ); + error = sqrt(error); + REQUIRE( error <= EPSTEST ); + + PORTABLE_FREE(device_data); +} +#endif + +int main(int argc, char* argv[]) { + +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::initialize(); +#endif + int result; + { + result = Catch::Session().run( argc, argv ); + } +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::finalize(); +#endif + return result; + +}