diff --git a/source/adios2/CMakeLists.txt b/source/adios2/CMakeLists.txt index e131db9c69..e8ed69f360 100644 --- a/source/adios2/CMakeLists.txt +++ b/source/adios2/CMakeLists.txt @@ -228,8 +228,8 @@ if(ADIOS2_HAVE_DataSpaces) add_subdirectory(toolkit/dataspaces) target_sources(adios2_core_mpi PRIVATE toolkit/dataspaces/ds_writer.c - engine/dataspaces/DataSpacesWriter.cpp engine/dataspaces/DataSpacesWriter.tcc - engine/dataspaces/DataSpacesReader.cpp engine/dataspaces/DataSpacesReader.tcc + engine/dataspaces/DataSpacesWriter.cpp engine/dataspaces/DataSpacesWriter.tcc + engine/dataspaces/DataSpacesReader.cpp engine/dataspaces/DataSpacesReader.tcc ) target_link_libraries(adios2_core_mpi PRIVATE DataSpaces::DataSpaces) endif() @@ -260,6 +260,7 @@ if(ADIOS2_HAVE_LIBPRESSIO) endif() if(ADIOS2_HAVE_MGARD) + target_sources(adios2_core PRIVATE operator/compress/CompressMGARDPlus.cpp) target_sources(adios2_core PRIVATE operator/compress/CompressMGARD.cpp) target_link_libraries(adios2_core PRIVATE MGARD::MGARD) endif() diff --git a/source/adios2/core/Operator.h b/source/adios2/core/Operator.h index f82ae0b788..fc99707d07 100644 --- a/source/adios2/core/Operator.h +++ b/source/adios2/core/Operator.h @@ -36,6 +36,7 @@ class Operator COMPRESS_SIRIUS = 5, COMPRESS_SZ = 6, COMPRESS_ZFP = 7, + COMPRESS_MGARDPLUS = 8, CALLBACK_SIGNATURE1 = 51, CALLBACK_SIGNATURE2 = 52, COMPRESS_NULL = 127, diff --git a/source/adios2/operator/OperatorFactory.cpp b/source/adios2/operator/OperatorFactory.cpp index 906a97699a..33a2385092 100644 --- a/source/adios2/operator/OperatorFactory.cpp +++ b/source/adios2/operator/OperatorFactory.cpp @@ -27,6 +27,7 @@ #ifdef ADIOS2_HAVE_MGARD #include "adios2/operator/compress/CompressMGARD.h" +#include "adios2/operator/compress/CompressMGARDPlus.h" #endif #ifdef ADIOS2_HAVE_PNG @@ -62,6 +63,8 @@ std::string OperatorTypeToString(const Operator::OperatorType type) return "libpressio"; case Operator::COMPRESS_MGARD: return "mgard"; + case Operator::COMPRESS_MGARDPLUS: + return "mgardplus"; case Operator::COMPRESS_PNG: return "png"; case Operator::COMPRESS_SIRIUS: @@ -104,6 +107,12 @@ std::shared_ptr MakeOperator(const std::string &type, { #ifdef ADIOS2_HAVE_MGARD ret = std::make_shared(parameters); +#endif + } + else if (typeLowerCase == "mgardplus") + { +#ifdef ADIOS2_HAVE_MGARD + ret = std::make_shared(parameters); #endif } else if (typeLowerCase == "png") diff --git a/source/adios2/operator/compress/CompressMGARDPlus.cpp b/source/adios2/operator/compress/CompressMGARDPlus.cpp new file mode 100644 index 0000000000..13dc8728ac --- /dev/null +++ b/source/adios2/operator/compress/CompressMGARDPlus.cpp @@ -0,0 +1,144 @@ +/* + * Distributed under the OSI-approved Apache License, Version 2.0. See + * accompanying file Copyright.txt for details. + * + * CompressMGARDPlus.cpp : + * + * Created on: Feb 23, 2022 + * Author: Jason Wang jason.ruonan.wang@gmail.com + */ + +#include "CompressMGARDPlus.h" +#include "CompressMGARD.h" +#include "adios2/helper/adiosFunctions.h" + +namespace adios2 +{ +namespace core +{ +namespace compress +{ + +CompressMGARDPlus::CompressMGARDPlus(const Params ¶meters) +: Operator("mgardplus", COMPRESS_MGARDPLUS, "compress", parameters) +{ +} + +size_t CompressMGARDPlus::Operate(const char *dataIn, const Dims &blockStart, + const Dims &blockCount, const DataType type, + char *bufferOut) +{ + size_t bufferOutOffset = 0; + const uint8_t bufferVersion = 1; + + MakeCommonHeader(bufferOut, bufferOutOffset, bufferVersion); + + bufferOutOffset += 32; // TODO: reserve memory space + + CompressMGARD mgard(m_Parameters); + size_t mgardBufferSize = mgard.Operate(dataIn, blockStart, blockCount, type, + bufferOut + bufferOutOffset); + + if (*reinterpret_cast(bufferOut + bufferOutOffset) == + COMPRESS_MGARD) + { + std::vector tmpDecompressBuffer( + helper::GetTotalSize(blockCount, helper::GetDataTypeSize(type))); + + mgard.InverseOperate(bufferOut + bufferOutOffset, mgardBufferSize, + tmpDecompressBuffer.data()); + + // TODO: now the original data is in dataIn, the compressed and then + // decompressed data is in tmpDecompressBuffer.data(). However, these + // are char pointers, you will need to convert them into right types as + // follows: + if (type == DataType::Double || type == DataType::DoubleComplex) + { + // TODO: use reinterpret_cast(dataIn) and + // reinterpret_cast(tmpDecompressBuffer.data()) + // to read original data and decompressed data + } + else if (type == DataType::Float || type == DataType::FloatComplex) + { + // TODO: use reinterpret_cast(dataIn) and + // reinterpret_cast(tmpDecompressBuffer.data()) + // to read original data and decompressed data + } + // TODO: after your algorithm is done, put the result into + // *reinterpret_cast(bufferOut+bufferOutOffset) for your first + // double number *reinterpret_cast(bufferOut+bufferOutOffset+8) + // for your second double number and so on + } + + bufferOutOffset += mgardBufferSize; + return bufferOutOffset; +} + +size_t CompressMGARDPlus::DecompressV1(const char *bufferIn, + const size_t sizeIn, char *dataOut) +{ + // Do NOT remove even if the buffer version is updated. Data might be still + // in lagacy formats. This function must be kept for backward compatibility. + // If a newer buffer format is implemented, create another function, e.g. + // DecompressV2 and keep this function for decompressing lagacy data. + + // TODO: read your results here from + // *reinterpret_cast(bufferIn) for your first double number + // *reinterpret_cast(bufferIn+8) for your second double number and + // so on + + size_t bufferInOffset = 32; // this number needs to be the same as the + // memory space you reserved in Operate() + + CompressMGARD mgard(m_Parameters); + size_t sizeOut = mgard.InverseOperate(bufferIn + bufferInOffset, + sizeIn - bufferInOffset, dataOut); + + // TODO: the regular decompressed buffer is in dataOut, with the size of + // sizeOut. Here you may want to do your magic to change the decompressed + // data somehow to improve its accuracy :) + + return sizeOut; +} + +size_t CompressMGARDPlus::InverseOperate(const char *bufferIn, + const size_t sizeIn, char *dataOut) +{ + size_t bufferInOffset = 1; // skip operator type + const uint8_t bufferVersion = + GetParameter(bufferIn, bufferInOffset); + bufferInOffset += 2; // skip two reserved bytes + + if (bufferVersion == 1) + { + return DecompressV1(bufferIn + bufferInOffset, sizeIn - bufferInOffset, + dataOut); + } + else if (bufferVersion == 2) + { + // TODO: if a Version 2 mgard buffer is being implemented, put it here + // and keep the DecompressV1 routine for backward compatibility + } + else + { + helper::Throw("Operator", "CompressMGARDPlus", + "InverseOperate", + "invalid mgard buffer version"); + } + + return 0; +} + +bool CompressMGARDPlus::IsDataTypeValid(const DataType type) const +{ + if (type == DataType::Double || type == DataType::Float || + type == DataType::DoubleComplex || type == DataType::FloatComplex) + { + return true; + } + return false; +} + +} // end namespace compress +} // end namespace core +} // end namespace adios2 diff --git a/source/adios2/operator/compress/CompressMGARDPlus.h b/source/adios2/operator/compress/CompressMGARDPlus.h new file mode 100644 index 0000000000..83b54c3665 --- /dev/null +++ b/source/adios2/operator/compress/CompressMGARDPlus.h @@ -0,0 +1,74 @@ +/* + * Distributed under the OSI-approved Apache License, Version 2.0. See + * accompanying file Copyright.txt for details. + * + * CompressMGARDPlus.h : + * + * Created on: Feb 23, 2022 + * Author: Jason Wang jason.ruonan.wang@gmail.com + */ + +#ifndef ADIOS2_OPERATOR_COMPRESS_COMPRESSMGARDPLUS_H_ +#define ADIOS2_OPERATOR_COMPRESS_COMPRESSMGARDPLUS_H_ + +#include "adios2/core/Operator.h" + +namespace adios2 +{ +namespace core +{ +namespace compress +{ + +class CompressMGARDPlus : public Operator +{ + +public: + CompressMGARDPlus(const Params ¶meters); + + ~CompressMGARDPlus() = default; + + /** + * @param dataIn + * @param blockStart + * @param blockCount + * @param type + * @param bufferOut + * @return size of compressed buffer + */ + size_t Operate(const char *dataIn, const Dims &blockStart, + const Dims &blockCount, const DataType type, + char *bufferOut) final; + + /** + * @param bufferIn + * @param sizeIn + * @param dataOut + * @return size of decompressed buffer + */ + size_t InverseOperate(const char *bufferIn, const size_t sizeIn, + char *dataOut) final; + + bool IsDataTypeValid(const DataType type) const final; + +private: + /** + * Decompress function for V1 buffer. Do NOT remove even if the buffer + * version is updated. Data might be still in lagacy formats. This function + * must be kept for backward compatibility + * @param bufferIn : compressed data buffer (V1 only) + * @param sizeIn : number of bytes in bufferIn + * @param dataOut : decompressed data buffer + * @return : number of bytes in dataOut + */ + size_t DecompressV1(const char *bufferIn, const size_t sizeIn, + char *dataOut); + + std::string m_VersionInfo; +}; + +} // end namespace compress +} // end namespace core +} // end namespace adios2 + +#endif /* ADIOS2_OPERATOR_COMPRESS_COMPRESSMGARDPLUS_H_ */ diff --git a/source/adios2/toolkit/format/bp/BPBase.cpp b/source/adios2/toolkit/format/bp/BPBase.cpp index 4408a80508..53d0318818 100644 --- a/source/adios2/toolkit/format/bp/BPBase.cpp +++ b/source/adios2/toolkit/format/bp/BPBase.cpp @@ -444,14 +444,20 @@ std::string BPBase::ReadBPString(const std::vector &buffer, // static members const std::set BPBase::m_TransformTypes = { {"unknown", "none", "identity", "bzip2", "sz", "zfp", "mgard", "png", - "blosc", "sirius"}}; + "blosc", "sirius", "mgardplus"}}; const std::map BPBase::m_TransformTypesToNames = { - {transform_unknown, "unknown"}, {transform_none, "none"}, - {transform_identity, "identity"}, {transform_sz, "sz"}, - {transform_zfp, "zfp"}, {transform_mgard, "mgard"}, - {transform_png, "png"}, {transform_bzip2, "bzip2"}, - {transform_blosc, "blosc"}, {transform_sirius, "sirius"}}; + {transform_unknown, "unknown"}, + {transform_none, "none"}, + {transform_identity, "identity"}, + {transform_sz, "sz"}, + {transform_zfp, "zfp"}, + {transform_mgard, "mgard"}, + {transform_png, "png"}, + {transform_bzip2, "bzip2"}, + {transform_blosc, "blosc"}, + {transform_sirius, "sirius"}, + {transform_mgardplus, "mgardplus"}}; BPBase::TransformTypes BPBase::TransformTypeEnum(const std::string transformType) const noexcept diff --git a/source/adios2/toolkit/format/bp/BPBase.h b/source/adios2/toolkit/format/bp/BPBase.h index cddbe68719..2dfb497c83 100644 --- a/source/adios2/toolkit/format/bp/BPBase.h +++ b/source/adios2/toolkit/format/bp/BPBase.h @@ -452,7 +452,8 @@ class BPBase transform_blosc = 11, transform_mgard = 12, transform_png = 13, - transform_sirius = 14 + transform_sirius = 14, + transform_mgardplus = 15 }; /** Supported transform types */ diff --git a/testing/adios2/engine/bp/operations/CMakeLists.txt b/testing/adios2/engine/bp/operations/CMakeLists.txt index ca6442da46..7bf6a90cb6 100644 --- a/testing/adios2/engine/bp/operations/CMakeLists.txt +++ b/testing/adios2/engine/bp/operations/CMakeLists.txt @@ -51,6 +51,7 @@ endif() if(ADIOS2_HAVE_MGARD) bp_gtest_add_tests_helper(WriteReadMGARD MPI_ALLOW) + bp_gtest_add_tests_helper(WriteReadMGARDPlus MPI_ALLOW) endif() if(ADIOS2_HAVE_BZip2) diff --git a/testing/adios2/engine/bp/operations/TestBPWriteReadMGARDPlus.cpp b/testing/adios2/engine/bp/operations/TestBPWriteReadMGARDPlus.cpp new file mode 100644 index 0000000000..bed949d225 --- /dev/null +++ b/testing/adios2/engine/bp/operations/TestBPWriteReadMGARDPlus.cpp @@ -0,0 +1,997 @@ +/* + * Distributed under the OSI-approved Apache License, Version 2.0. See + * accompanying file Copyright.txt for details. + */ +#include +#include + +#include +#include //std::cout +#include //std::iota +#include + +#include + +#include + +std::string engineName; // comes from command line + +void MGARDAccuracy1D(const std::string tolerance) +{ + // Each process would write a 1x8 array and all processes would + // form a mpiSize * Nx 1D array + const std::string fname("BPWRMGARD1D_" + tolerance + ".bp"); + + int mpiRank = 0, mpiSize = 1; + // Number of rows + const size_t Nx = 100; + + // Number of steps + const size_t NSteps = 1; + + std::vector r32s(Nx); + std::vector r64s(Nx); + + // range 0 to 100*50 + std::iota(r32s.begin(), r32s.end(), 0.f); + std::iota(r64s.begin(), r64s.end(), 0.); + +#if ADIOS2_USE_MPI + MPI_Comm_rank(MPI_COMM_WORLD, &mpiRank); + MPI_Comm_size(MPI_COMM_WORLD, &mpiSize); +#endif + +#if ADIOS2_USE_MPI + adios2::ADIOS adios(MPI_COMM_WORLD); +#else + adios2::ADIOS adios; +#endif + { + adios2::IO io = adios.DeclareIO("TestIO"); + + if (!engineName.empty()) + { + io.SetEngine(engineName); + } + + const adios2::Dims shape{static_cast(Nx * mpiSize)}; + const adios2::Dims start{static_cast(Nx * mpiRank)}; + const adios2::Dims count{Nx}; + + auto var_r32 = io.DefineVariable("r32", shape, start, count, + adios2::ConstantDims); + auto var_r64 = io.DefineVariable("r64", shape, start, count, + adios2::ConstantDims); + + var_r32.AddOperation("MGARDPlus", + {{adios2::ops::mgard::key::tolerance, tolerance}}); + var_r64.AddOperation("MGARDPlus", + {{adios2::ops::mgard::key::tolerance, tolerance}}); + + adios2::Engine bpWriter = io.Open(fname, adios2::Mode::Write); + + for (size_t step = 0; step < NSteps; ++step) + { + bpWriter.BeginStep(); + // bpWriter.Put("r32", r32s.data()); + bpWriter.Put("r64", r64s.data()); + bpWriter.EndStep(); + } + + bpWriter.Close(); + } + + { + adios2::IO io = adios.DeclareIO("ReadIO"); + + if (!engineName.empty()) + { + io.SetEngine(engineName); + } + + adios2::Engine bpReader = io.Open(fname, adios2::Mode::Read); + + unsigned int t = 0; + std::vector decompressedR32s; + std::vector decompressedR64s; + + while (bpReader.BeginStep() == adios2::StepStatus::OK) + { + // auto var_r32 = io.InquireVariable("r32"); + // EXPECT_TRUE(var_r32); + // ASSERT_EQ(var_r32.ShapeID(), + // adios2::ShapeID::GlobalArray); ASSERT_EQ(var_r32.Steps(), + // NSteps); ASSERT_EQ(var_r32.Shape()[0], mpiSize * Nx); + // ASSERT_EQ(var_r32.Shape()[1], Ny); + + auto var_r64 = io.InquireVariable("r64"); + EXPECT_TRUE(var_r64); + ASSERT_EQ(var_r64.ShapeID(), adios2::ShapeID::GlobalArray); + ASSERT_EQ(var_r64.Steps(), NSteps); + ASSERT_EQ(var_r64.Shape()[0], mpiSize * Nx); + + const adios2::Dims start{mpiRank * Nx}; + const adios2::Dims count{Nx}; + const adios2::Box sel(start, count); + // var_r32.SetSelection(sel); + var_r64.SetSelection(sel); + + // bpReader.Get(var_r32, decompressedR32s); + bpReader.Get(var_r64, decompressedR64s); + bpReader.EndStep(); + + double maxDiff = 0; + + for (size_t i = 0; i < Nx; ++i) + { + std::stringstream ss; + ss << "t=" << t << " i=" << i << " rank=" << mpiRank; + std::string msg = ss.str(); + + double diff = std::abs(r64s[i] - decompressedR64s[i]); + + if (diff > maxDiff) + { + maxDiff = diff; + } + } + ++t; + + auto itMax = std::max_element(r64s.begin(), r64s.end()); + + const double relativeMaxDiff = maxDiff / *itMax; + ASSERT_LT(relativeMaxDiff, std::stod(tolerance)); + std::cout << "Relative Max Diff " << relativeMaxDiff + << " tolerance " << tolerance << "\n"; + } + + EXPECT_EQ(t, NSteps); + + bpReader.Close(); + } +} + +void MGARDAccuracy2D(const std::string tolerance) +{ + // Each process would write a 1x8 array and all processes would + // form a mpiSize * Nx 1D array + const std::string fname("BPWRMGARD2D_" + tolerance + ".bp"); + + int mpiRank = 0, mpiSize = 1; + // Number of rows + const size_t Nx = 100; + const size_t Ny = 50; + + // Number of steps + const size_t NSteps = 1; + + std::vector r32s(Nx * Ny); + std::vector r64s(Nx * Ny); + + // range 0 to 100*50 + std::iota(r32s.begin(), r32s.end(), 0.f); + std::iota(r64s.begin(), r64s.end(), 0.); + +#if ADIOS2_USE_MPI + MPI_Comm_rank(MPI_COMM_WORLD, &mpiRank); + MPI_Comm_size(MPI_COMM_WORLD, &mpiSize); +#endif + +#if ADIOS2_USE_MPI + adios2::ADIOS adios(MPI_COMM_WORLD); +#else + adios2::ADIOS adios; +#endif + { + adios2::IO io = adios.DeclareIO("TestIO"); + + if (!engineName.empty()) + { + io.SetEngine(engineName); + } + + const adios2::Dims shape{static_cast(Nx * mpiSize), Ny}; + const adios2::Dims start{static_cast(Nx * mpiRank), 0}; + const adios2::Dims count{Nx, Ny}; + + auto var_r32 = io.DefineVariable("r32", shape, start, count, + adios2::ConstantDims); + auto var_r64 = io.DefineVariable("r64", shape, start, count, + adios2::ConstantDims); + + var_r32.AddOperation("MGARDPlus", + {{adios2::ops::mgard::key::tolerance, tolerance}}); + var_r64.AddOperation("MGARDPlus", + {{adios2::ops::mgard::key::tolerance, tolerance}}); + + adios2::Engine bpWriter = io.Open(fname, adios2::Mode::Write); + + for (size_t step = 0; step < NSteps; ++step) + { + bpWriter.BeginStep(); + // bpWriter.Put("r32", r32s.data()); + bpWriter.Put("r64", r64s.data()); + bpWriter.EndStep(); + } + + bpWriter.Close(); + } + + { + adios2::IO io = adios.DeclareIO("ReadIO"); + + if (!engineName.empty()) + { + io.SetEngine(engineName); + } + + adios2::Engine bpReader = io.Open(fname, adios2::Mode::Read); + + unsigned int t = 0; + std::vector decompressedR32s; + std::vector decompressedR64s; + + while (bpReader.BeginStep() == adios2::StepStatus::OK) + { + // auto var_r32 = io.InquireVariable("r32"); + // EXPECT_TRUE(var_r32); + // ASSERT_EQ(var_r32.ShapeID(), + // adios2::ShapeID::GlobalArray); ASSERT_EQ(var_r32.Steps(), + // NSteps); ASSERT_EQ(var_r32.Shape()[0], mpiSize * Nx); + // ASSERT_EQ(var_r32.Shape()[1], Ny); + + auto var_r64 = io.InquireVariable("r64"); + EXPECT_TRUE(var_r64); + ASSERT_EQ(var_r64.ShapeID(), adios2::ShapeID::GlobalArray); + ASSERT_EQ(var_r64.Steps(), NSteps); + ASSERT_EQ(var_r64.Shape()[0], mpiSize * Nx); + ASSERT_EQ(var_r64.Shape()[1], Ny); + + const adios2::Dims start{mpiRank * Nx, 0}; + const adios2::Dims count{Nx, Ny}; + const adios2::Box sel(start, count); + // var_r32.SetSelection(sel); + var_r64.SetSelection(sel); + + // bpReader.Get(var_r32, decompressedR32s); + bpReader.Get(var_r64, decompressedR64s); + bpReader.EndStep(); + + double maxDiff = 0; + + for (size_t i = 0; i < Nx * Ny; ++i) + { + std::stringstream ss; + ss << "t=" << t << " i=" << i << " rank=" << mpiRank; + std::string msg = ss.str(); + + double diff = std::abs(r64s[i] - decompressedR64s[i]); + + if (diff > maxDiff) + { + maxDiff = diff; + } + } + ++t; + + auto itMax = std::max_element(r64s.begin(), r64s.end()); + + const double relativeMaxDiff = maxDiff / *itMax; + ASSERT_LT(relativeMaxDiff, std::stod(tolerance)); + std::cout << "Relative Max Diff " << relativeMaxDiff + << " tolerance " << tolerance << "\n"; + } + + EXPECT_EQ(t, NSteps); + + bpReader.Close(); + } +} + +void MGARDAccuracy3D(const std::string tolerance) +{ + // Each process would write a 1x8 array and all processes would + // form a mpiSize * Nx 1D array + const std::string fname("BPWRMGARD3D_" + tolerance + ".bp"); + + int mpiRank = 0, mpiSize = 1; + // Number of rows + const size_t Nx = 100; + const size_t Ny = 50; + const size_t Nz = 15; + + // Number of steps + const size_t NSteps = 1; + + std::vector r32s(Nx * Ny * Nz); + std::vector r64s(Nx * Ny * Nz); + + // range 0 to 100*50 + std::iota(r32s.begin(), r32s.end(), 0.f); + std::iota(r64s.begin(), r64s.end(), 0.); + +#if ADIOS2_USE_MPI + MPI_Comm_rank(MPI_COMM_WORLD, &mpiRank); + MPI_Comm_size(MPI_COMM_WORLD, &mpiSize); +#endif + +#if ADIOS2_USE_MPI + adios2::ADIOS adios(MPI_COMM_WORLD); +#else + adios2::ADIOS adios; +#endif + { + adios2::IO io = adios.DeclareIO("TestIO"); + + if (!engineName.empty()) + { + io.SetEngine(engineName); + } + + const adios2::Dims shape{static_cast(Nx * mpiSize), Ny, Nz}; + const adios2::Dims start{static_cast(Nx * mpiRank), 0, 0}; + const adios2::Dims count{Nx, Ny, Nz}; + + auto var_r32 = io.DefineVariable("r32", shape, start, count, + adios2::ConstantDims); + auto var_r64 = io.DefineVariable("r64", shape, start, count, + adios2::ConstantDims); + + var_r32.AddOperation("MGARDPlus", + {{adios2::ops::mgard::key::tolerance, tolerance}}); + var_r64.AddOperation("MGARDPlus", + {{adios2::ops::mgard::key::tolerance, tolerance}}); + + adios2::Engine bpWriter = io.Open(fname, adios2::Mode::Write); + + for (size_t step = 0; step < NSteps; ++step) + { + bpWriter.BeginStep(); + // bpWriter.Put("r32", r32s.data()); + bpWriter.Put("r64", r64s.data()); + bpWriter.EndStep(); + } + + bpWriter.Close(); + } + + { + adios2::IO io = adios.DeclareIO("ReadIO"); + + if (!engineName.empty()) + { + io.SetEngine(engineName); + } + + adios2::Engine bpReader = io.Open(fname, adios2::Mode::Read); + + unsigned int t = 0; + std::vector decompressedR32s; + std::vector decompressedR64s; + + while (bpReader.BeginStep() == adios2::StepStatus::OK) + { + // auto var_r32 = io.InquireVariable("r32"); + // EXPECT_TRUE(var_r32); + // ASSERT_EQ(var_r32.ShapeID(), + // adios2::ShapeID::GlobalArray); ASSERT_EQ(var_r32.Steps(), + // NSteps); ASSERT_EQ(var_r32.Shape()[0], mpiSize * Nx); + // ASSERT_EQ(var_r32.Shape()[1], Ny); + // ASSERT_EQ(var_r32.Shape()[2], Nz); + + auto var_r64 = io.InquireVariable("r64"); + EXPECT_TRUE(var_r64); + ASSERT_EQ(var_r64.ShapeID(), adios2::ShapeID::GlobalArray); + ASSERT_EQ(var_r64.Steps(), NSteps); + ASSERT_EQ(var_r64.Shape()[0], mpiSize * Nx); + ASSERT_EQ(var_r64.Shape()[1], Ny); + ASSERT_EQ(var_r64.Shape()[2], Nz); + + const adios2::Dims start{mpiRank * Nx, 0, 0}; + const adios2::Dims count{Nx, Ny, Nz}; + const adios2::Box sel(start, count); + // var_r32.SetSelection(sel); + var_r64.SetSelection(sel); + + // bpReader.Get(var_r32, decompressedR32s); + bpReader.Get(var_r64, decompressedR64s); + bpReader.EndStep(); + + double maxDiff = 0; + + for (size_t i = 0; i < Nx * Ny * Nz; ++i) + { + std::stringstream ss; + ss << "t=" << t << " i=" << i << " rank=" << mpiRank; + std::string msg = ss.str(); + + double diff = std::abs(r64s[i] - decompressedR64s[i]); + + if (diff > maxDiff) + { + maxDiff = diff; + } + } + ++t; + + auto itMax = std::max_element(r64s.begin(), r64s.end()); + + const double relativeMaxDiff = maxDiff / *itMax; + ASSERT_LT(relativeMaxDiff, std::stod(tolerance)); + std::cout << "Relative Max Diff " << relativeMaxDiff + << " tolerance " << tolerance << "\n"; + } + + EXPECT_EQ(t, NSteps); + + bpReader.Close(); + } +} + +void MGARDAccuracy1DSel(const std::string tolerance) +{ + // Each process would write a 1x8 array and all processes would + // form a mpiSize * Nx 1D array + const std::string fname("BPWRMGARD1DSel_" + tolerance + ".bp"); + + int mpiRank = 0, mpiSize = 1; + // Number of rows + const size_t Nx = 1000; + + // Number of steps + const size_t NSteps = 1; + + std::vector r32s(Nx); + std::vector r64s(Nx); + + // range 0 to 999 + std::iota(r32s.begin(), r32s.end(), 0.f); + std::iota(r64s.begin(), r64s.end(), 0.); + +#if ADIOS2_USE_MPI + MPI_Comm_rank(MPI_COMM_WORLD, &mpiRank); + MPI_Comm_size(MPI_COMM_WORLD, &mpiSize); +#endif + +#if ADIOS2_USE_MPI + adios2::ADIOS adios(MPI_COMM_WORLD); +#else + adios2::ADIOS adios; +#endif + { + adios2::IO io = adios.DeclareIO("TestIO"); + + if (!engineName.empty()) + { + io.SetEngine(engineName); + } + + const adios2::Dims shape{static_cast(Nx * mpiSize)}; + const adios2::Dims start{static_cast(Nx * mpiRank)}; + const adios2::Dims count{Nx}; + + auto var_r32 = io.DefineVariable("r32", shape, start, count, + adios2::ConstantDims); + auto var_r64 = io.DefineVariable("r64", shape, start, count, + adios2::ConstantDims); + + var_r32.AddOperation("MGARDPlus", + {{adios2::ops::mgard::key::tolerance, tolerance}}); + var_r64.AddOperation("MGARDPlus", + {{adios2::ops::mgard::key::tolerance, tolerance}}); + + adios2::Engine bpWriter = io.Open(fname, adios2::Mode::Write); + + for (size_t step = 0; step < NSteps; ++step) + { + bpWriter.BeginStep(); + // bpWriter.Put("r32", r32s.data()); + bpWriter.Put("r64", r64s.data()); + bpWriter.EndStep(); + } + + bpWriter.Close(); + } + + { + adios2::IO io = adios.DeclareIO("ReadIO"); + + if (!engineName.empty()) + { + io.SetEngine(engineName); + } + + adios2::Engine bpReader = io.Open(fname, adios2::Mode::Read); + + unsigned int t = 0; + std::vector decompressedR32s; + std::vector decompressedR64s; + + while (bpReader.BeginStep() == adios2::StepStatus::OK) + { + // auto var_r32 = io.InquireVariable("r32"); + // EXPECT_TRUE(var_r32); + // ASSERT_EQ(var_r32.ShapeID(), + // adios2::ShapeID::GlobalArray); ASSERT_EQ(var_r32.Steps(), + // NSteps); ASSERT_EQ(var_r32.Shape()[0], mpiSize * Nx); + + auto var_r64 = io.InquireVariable("r64"); + EXPECT_TRUE(var_r64); + ASSERT_EQ(var_r64.ShapeID(), adios2::ShapeID::GlobalArray); + ASSERT_EQ(var_r64.Steps(), NSteps); + ASSERT_EQ(var_r64.Shape()[0], mpiSize * Nx); + + const adios2::Dims start{mpiRank * Nx + Nx / 2}; + const adios2::Dims count{Nx / 2}; + const adios2::Box sel(start, count); + // var_r32.SetSelection(sel); + var_r64.SetSelection(sel); + + // bpReader.Get(var_r32, decompressedR32s); + bpReader.Get(var_r64, decompressedR64s); + bpReader.EndStep(); + + for (size_t i = 0; i < Nx / 2; ++i) + { + std::stringstream ss; + ss << "t=" << t << " i=" << i << " rank=" << mpiRank; + std::string msg = ss.str(); + + // ASSERT_LT(std::abs(decompressedR32s[i] - + // r32s[Nx / 2 + i]), + // tolerance) + // << msg; + ASSERT_LT(std::abs(decompressedR64s[i] - r64s[Nx / 2 + i]), + std::stod(tolerance)) + << msg; + } + ++t; + } + + EXPECT_EQ(t, NSteps); + + bpReader.Close(); + } +} + +void MGARDAccuracy2DSel(const std::string tolerance) +{ + // Each process would write a 1x8 array and all processes would + // form a mpiSize * Nx 1D array + const std::string fname("BPWRMGARD2DSel_" + tolerance + ".bp"); + + int mpiRank = 0, mpiSize = 1; + // Number of rows + const size_t Nx = 100; + const size_t Ny = 50; + + // Number of steps + const size_t NSteps = 1; + + std::vector r32s(Nx * Ny); + std::vector r64s(Nx * Ny); + + // range 0 to 100*50 + std::iota(r32s.begin(), r32s.end(), 0.f); + std::iota(r64s.begin(), r64s.end(), 0.); + +#if ADIOS2_USE_MPI + MPI_Comm_rank(MPI_COMM_WORLD, &mpiRank); + MPI_Comm_size(MPI_COMM_WORLD, &mpiSize); +#endif + +#if ADIOS2_USE_MPI + adios2::ADIOS adios(MPI_COMM_WORLD); +#else + adios2::ADIOS adios; +#endif + { + adios2::IO io = adios.DeclareIO("TestIO"); + + if (!engineName.empty()) + { + io.SetEngine(engineName); + } + + const adios2::Dims shape{static_cast(Nx * mpiSize), Ny}; + const adios2::Dims start{static_cast(Nx * mpiRank), 0}; + const adios2::Dims count{Nx, Ny}; + + auto var_r32 = io.DefineVariable("r32", shape, start, count, + adios2::ConstantDims); + auto var_r64 = io.DefineVariable("r64", shape, start, count, + adios2::ConstantDims); + + var_r32.AddOperation("MGARDPlus", + {{adios2::ops::mgard::key::tolerance, tolerance}}); + var_r64.AddOperation("MGARDPlus", + {{adios2::ops::mgard::key::tolerance, tolerance}}); + + adios2::Engine bpWriter = io.Open(fname, adios2::Mode::Write); + + for (size_t step = 0; step < NSteps; ++step) + { + bpWriter.BeginStep(); + // bpWriter.Put("r32", r32s.data()); + bpWriter.Put("r64", r64s.data()); + bpWriter.EndStep(); + } + + bpWriter.Close(); + } + + { + adios2::IO io = adios.DeclareIO("ReadIO"); + + if (!engineName.empty()) + { + io.SetEngine(engineName); + } + + adios2::Engine bpReader = io.Open(fname, adios2::Mode::Read); + + unsigned int t = 0; + std::vector decompressedR32s; + std::vector decompressedR64s; + + while (bpReader.BeginStep() == adios2::StepStatus::OK) + { + // auto var_r32 = io.InquireVariable("r32"); + // EXPECT_TRUE(var_r32); + // ASSERT_EQ(var_r32.ShapeID(), + // adios2::ShapeID::GlobalArray); ASSERT_EQ(var_r32.Steps(), + // NSteps); ASSERT_EQ(var_r32.Shape()[0], mpiSize * Nx); + // ASSERT_EQ(var_r32.Shape()[1], Ny); + + auto var_r64 = io.InquireVariable("r64"); + EXPECT_TRUE(var_r64); + ASSERT_EQ(var_r64.ShapeID(), adios2::ShapeID::GlobalArray); + ASSERT_EQ(var_r64.Steps(), NSteps); + ASSERT_EQ(var_r64.Shape()[0], mpiSize * Nx); + ASSERT_EQ(var_r64.Shape()[1], Ny); + + const adios2::Dims start{mpiRank * Nx + Nx / 2, 0}; + const adios2::Dims count{Nx / 2, Ny}; + const adios2::Box sel(start, count); + // var_r32.SetSelection(sel); + var_r64.SetSelection(sel); + + // bpReader.Get(var_r32, decompressedR32s); + bpReader.Get(var_r64, decompressedR64s); + bpReader.EndStep(); + + for (size_t i = 0; i < Nx / 2 * Ny; ++i) + { + std::stringstream ss; + ss << "t=" << t << " i=" << i << " rank=" << mpiRank; + std::string msg = ss.str(); + + // ASSERT_LT(std::abs(decompressedR32s[i] - + // r32s[Nx / 2 * Ny + i]), + // tolerance) + // << msg; + ASSERT_LT(std::abs(decompressedR64s[i] - r64s[Nx / 2 * Ny + i]), + std::stod(tolerance)) + << msg; + } + ++t; + } + + EXPECT_EQ(t, NSteps); + + bpReader.Close(); + } +} + +void MGARDAccuracy3DSel(const std::string tolerance) +{ + // Each process would write a 1x8 array and all processes would + // form a mpiSize * Nx 1D array + const std::string fname("BPWRMGARD3DSel_" + tolerance + ".bp"); + + int mpiRank = 0, mpiSize = 1; + // Number of rows + const size_t Nx = 10; + const size_t Ny = 20; + const size_t Nz = 15; + + // Number of steps + const size_t NSteps = 1; + + std::vector r32s(Nx * Ny * Nz); + std::vector r64s(Nx * Ny * Nz); + + // range 0 to 100*50 + std::iota(r32s.begin(), r32s.end(), 0.f); + std::iota(r64s.begin(), r64s.end(), 0.); + +#if ADIOS2_USE_MPI + MPI_Comm_rank(MPI_COMM_WORLD, &mpiRank); + MPI_Comm_size(MPI_COMM_WORLD, &mpiSize); +#endif + +#if ADIOS2_USE_MPI + adios2::ADIOS adios(MPI_COMM_WORLD); +#else + adios2::ADIOS adios; +#endif + { + adios2::IO io = adios.DeclareIO("TestIO"); + + if (!engineName.empty()) + { + io.SetEngine(engineName); + } + + const adios2::Dims shape{static_cast(Nx * mpiSize), Ny, Nz}; + const adios2::Dims start{static_cast(Nx * mpiRank), 0, 0}; + const adios2::Dims count{Nx, Ny, Nz}; + + auto var_r32 = io.DefineVariable("r32", shape, start, count, + adios2::ConstantDims); + auto var_r64 = io.DefineVariable("r64", shape, start, count, + adios2::ConstantDims); + + var_r32.AddOperation("MGARDPlus", + {{adios2::ops::mgard::key::tolerance, tolerance}}); + var_r64.AddOperation("MGARDPlus", + {{adios2::ops::mgard::key::tolerance, tolerance}}); + + adios2::Engine bpWriter = io.Open(fname, adios2::Mode::Write); + + for (size_t step = 0; step < NSteps; ++step) + { + bpWriter.BeginStep(); + // bpWriter.Put("r32", r32s.data()); + bpWriter.Put("r64", r64s.data()); + bpWriter.EndStep(); + } + + bpWriter.Close(); + } + + { + adios2::IO io = adios.DeclareIO("ReadIO"); + + if (!engineName.empty()) + { + io.SetEngine(engineName); + } + + adios2::Engine bpReader = io.Open(fname, adios2::Mode::Read); + unsigned int t = 0; + std::vector decompressedR32s; + std::vector decompressedR64s; + + while (bpReader.BeginStep() == adios2::StepStatus::OK) + { + + // auto var_r32 = io.InquireVariable("r32"); + // EXPECT_TRUE(var_r32); + // ASSERT_EQ(var_r32.ShapeID(), + // adios2::ShapeID::GlobalArray); ASSERT_EQ(var_r32.Steps(), + // NSteps); ASSERT_EQ(var_r32.Shape()[0], mpiSize * Nx); + // ASSERT_EQ(var_r32.Shape()[1], Ny); + // ASSERT_EQ(var_r32.Shape()[2], Nz); + + auto var_r64 = io.InquireVariable("r64"); + EXPECT_TRUE(var_r64); + ASSERT_EQ(var_r64.ShapeID(), adios2::ShapeID::GlobalArray); + ASSERT_EQ(var_r64.Steps(), NSteps); + ASSERT_EQ(var_r64.Shape()[0], mpiSize * Nx); + ASSERT_EQ(var_r64.Shape()[1], Ny); + ASSERT_EQ(var_r64.Shape()[2], Nz); + + const adios2::Dims start{mpiRank * Nx + Nx / 2, 0, 0}; + const adios2::Dims count{Nx / 2, Ny, Nz}; + const adios2::Box sel(start, count); + // var_r32.SetSelection(sel); + var_r64.SetSelection(sel); + + // bpReader.Get(var_r32, decompressedR32s); + bpReader.Get(var_r64, decompressedR64s); + bpReader.EndStep(); + + double maxDiff = 0; + + for (size_t i = 0; i < Nx / 2 * Ny * Nz; ++i) + { + std::stringstream ss; + ss << "t=" << t << " i=" << i << " rank=" << mpiRank; + std::string msg = ss.str(); + + // ASSERT_LT( + // std::abs(decompressedR32s[i] - r32s[Nx / 2 + // * Ny * Nz + i]), + // tolerance) + // << msg; + double diff = + std::abs(r64s[Nx / 2 * Ny * Nz + i] - decompressedR64s[i]); + + if (diff > maxDiff) + { + maxDiff = diff; + } + } + + auto itMax = std::max_element(r64s.begin(), r64s.end()); + ASSERT_LT(maxDiff / *itMax, std::stod(tolerance)); + + ++t; + } + + EXPECT_EQ(t, NSteps); + + bpReader.Close(); + } +} + +void MGARDAccuracy2DSmallSel(const std::string tolerance) +{ + // Each process would write a 1x8 array and all processes would + // form a mpiSize * Nx 1D array + const std::string fname("BPWRMGARD2DSmallSel_" + tolerance + ".bp"); + + int mpiRank = 0, mpiSize = 1; + // Number of rows + const size_t Nx = 5; + const size_t Ny = 5; + + // Number of steps + const size_t NSteps = 1; + + std::vector r32s = {0.00, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, + 0.07, 0.08, 0.09, 0.10, 0.11, 0.12, 0.13, + 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.20, + 0.21, 0.22, 0.23, 0.24}; + std::vector r64s = {0.00, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, + 0.07, 0.08, 0.09, 0.10, 0.11, 0.12, 0.13, + 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.20, + 0.21, 0.22, 0.23, 0.24}; +#if ADIOS2_USE_MPI + MPI_Comm_rank(MPI_COMM_WORLD, &mpiRank); + MPI_Comm_size(MPI_COMM_WORLD, &mpiSize); +#endif + +#if ADIOS2_USE_MPI + adios2::ADIOS adios(MPI_COMM_WORLD); +#else + adios2::ADIOS adios; +#endif + { + adios2::IO io = adios.DeclareIO("TestIO"); + + if (!engineName.empty()) + { + io.SetEngine(engineName); + } + + const adios2::Dims shape{static_cast(Nx * mpiSize), Ny}; + const adios2::Dims start{static_cast(Nx * mpiRank), 0}; + const adios2::Dims count{Nx, Ny}; + + auto var_r32 = io.DefineVariable("r32", shape, start, count, + adios2::ConstantDims); + auto var_r64 = io.DefineVariable("r64", shape, start, count, + adios2::ConstantDims); + + var_r32.AddOperation("MGARDPlus", + {{adios2::ops::mgard::key::tolerance, tolerance}}); + var_r64.AddOperation("MGARDPlus", + {{adios2::ops::mgard::key::tolerance, tolerance}}); + + adios2::Engine bpWriter = io.Open(fname, adios2::Mode::Write); + + for (size_t step = 0; step < NSteps; ++step) + { + bpWriter.BeginStep(); + // bpWriter.Put("r32", r32s.data()); + bpWriter.Put("r64", r64s.data()); + bpWriter.EndStep(); + } + + bpWriter.Close(); + } + + { + adios2::IO io = adios.DeclareIO("ReadIO"); + + if (!engineName.empty()) + { + io.SetEngine(engineName); + } + + adios2::Engine bpReader = io.Open(fname, adios2::Mode::Read); + + unsigned int t = 0; + std::vector decompressedR32s; + std::vector decompressedR64s; + + while (bpReader.BeginStep() == adios2::StepStatus::OK) + { + // auto var_r32 = io.InquireVariable("r32"); + // EXPECT_TRUE(var_r32); + // ASSERT_EQ(var_r32.ShapeID(), + // adios2::ShapeID::GlobalArray); ASSERT_EQ(var_r32.Steps(), + // NSteps); ASSERT_EQ(var_r32.Shape()[0], mpiSize * Nx); + // ASSERT_EQ(var_r32.Shape()[1], Ny); + + auto var_r64 = io.InquireVariable("r64"); + EXPECT_TRUE(var_r64); + ASSERT_EQ(var_r64.ShapeID(), adios2::ShapeID::GlobalArray); + ASSERT_EQ(var_r64.Steps(), NSteps); + ASSERT_EQ(var_r64.Shape()[0], mpiSize * Nx); + ASSERT_EQ(var_r64.Shape()[1], Ny); + + const adios2::Dims start{static_cast(mpiRank) * Nx + 1, + 1}; + const adios2::Dims count{2, 2}; + const adios2::Box sel(start, count); + // var_r32.SetSelection(sel); + var_r64.SetSelection(sel); + + // bpReader.Get(var_r32, decompressedR32s); + bpReader.Get(var_r64, decompressedR64s); + bpReader.EndStep(); + + // ASSERT_LT(std::abs(decompressedR32s[0] - 0.06), tolerance); + ASSERT_LT(std::abs(decompressedR64s[0] - 0.06), + std::stod(tolerance)); + + // ASSERT_LT(std::abs(decompressedR32s[1] - 0.07), tolerance); + ASSERT_LT(std::abs(decompressedR64s[1] - 0.07), + std::stod(tolerance)); + + // ASSERT_LT(std::abs(decompressedR32s[2] - 0.11), tolerance); + ASSERT_LT(std::abs(decompressedR64s[2] - 0.11), + std::stod(tolerance)); + + // ASSERT_LT(std::abs(decompressedR32s[3] - 0.12), tolerance); + ASSERT_LT(std::abs(decompressedR64s[3] - 0.12), + std::stod(tolerance)); + + ++t; + } + + EXPECT_EQ(t, NSteps); + + bpReader.Close(); + } +} + +class BPWriteReadMGARDPlus : public ::testing::TestWithParam +{ +public: + BPWriteReadMGARDPlus() = default; + virtual void SetUp(){}; + virtual void TearDown(){}; +}; + +TEST_P(BPWriteReadMGARDPlus, BPWRMGARD2D) { MGARDAccuracy2D(GetParam()); } + +TEST_P(BPWriteReadMGARDPlus, BPWRMGARD3D) { MGARDAccuracy3D(GetParam()); } + +INSTANTIATE_TEST_SUITE_P(MGARDAccuracy, BPWriteReadMGARDPlus, + ::testing::Values("0.01", "0.001", "0.0001", + "0.00001")); + +int main(int argc, char **argv) +{ +#if ADIOS2_USE_MPI + MPI_Init(nullptr, nullptr); +#endif + + int result; + ::testing::InitGoogleTest(&argc, argv); + if (argc > 1) + { + engineName = std::string(argv[1]); + } + result = RUN_ALL_TESTS(); + +#if ADIOS2_USE_MPI + MPI_Finalize(); +#endif + + return result; +}