diff --git a/.github/workflows/clang.yml b/.github/workflows/clang.yml index b4886ff9623..ce31db83a23 100644 --- a/.github/workflows/clang.yml +++ b/.github/workflows/clang.yml @@ -166,6 +166,82 @@ jobs: ccache -s du -hs ~/.cache/ccache + test_openpmd: + name: openPMD I/O Test + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - name: Dependencies + run: | + .github/workflows/dependencies/dependencies_clang.sh 15 + .github/workflows/dependencies/dependencies_clang-tidy.sh 15 + .github/workflows/dependencies/dependencies_ccache.sh + - name: Set Up Cache + uses: actions/cache@v4 + with: + path: ~/.cache/ccache + key: ccache-${{ github.workflow }}-${{ github.job }}-git-${{ github.sha }} + restore-keys: | + ccache-${{ github.workflow }}-${{ github.job }}-git- + - name: Install openPMD + run: | + export CCACHE_COMPRESS=1 + export CCACHE_COMPRESSLEVEL=10 + export CCACHE_MAXSIZE=40M + + sudo apt-get install -y --no-install-recommends libhdf5-mpi-dev + wget -q https://github.com/openPMD/openPMD-api/archive/refs/tags/0.15.2.tar.gz + tar xfz 0.15.2.tar.gz + cd openPMD-api-0.15.2 + mkdir build + cd build + cmake .. \ + -DCMAKE_CXX_COMPILER_LAUNCHER=ccache \ + -DopenPMD_BUILD_EXAMPLES=OFF \ + -DopenPMD_BUILD_CLI_TOOLS=OFF \ + -DopenPMD_BUILD_TESTING=OFF \ + -DopenPMD_USE_PYTHON=OFF \ + make -j4 + sudo make install + - name: Build & Test + env: {CXXFLAGS: "-fno-operator-names -Werror -Wall -Wextra -Wpedantic -Wnull-dereference -Wfloat-conversion -Wshadow -Woverloaded-virtual -Wextra-semi -Wunreachable-code -Wnon-virtual-dtor"} + run: | + export CCACHE_COMPRESS=1 + export CCACHE_COMPRESSLEVEL=10 + export CCACHE_MAXSIZE=40M + export CCACHE_EXTRAFILES=${{ github.workspace }}/.clang-tidy + export CCACHE_LOGFILE=${{ github.workspace }}/ccache.log.txt + ccache -z + + mkdir build + cd build + cmake .. \ + -DAMReX_OPENPMD=ON \ + -DCMAKE_VERBOSE_MAKEFILE=ON \ + -DAMReX_ENABLE_TESTS=ON \ + -DAMReX_SPACEDIM=3 \ + -DAMReX_MPI=ON \ + -DAMReX_OMP=OFF \ + -DAMReX_FORTRAN=OFF \ + -DAMReX_EB=OFF \ + -DAMReX_LINEAR_SOLVERS=OFF \ + -DAMReX_AMRLEVEL=OFF \ + -DAMReX_PARTICLES=ON \ + -DCMAKE_C_COMPILER=$(which clang-15) \ + -DCMAKE_CXX_COMPILER=$(which clang++-15) \ + -DCMAKE_CXX_COMPILER_LAUNCHER=ccache + make -j 4 + + ${{github.workspace}}/Tools/C_scripts/mmclt.py --input ${{github.workspace}}/ccache.log.txt + make -j4 -k -f clang-tidy-ccache-misses.mak \ + CLANG_TIDY=clang-tidy-15 \ + CLANG_TIDY_ARGS="--config-file=${{github.workspace}}/.clang-tidy --warnings-as-errors=*" + + ctest --output-on-failure + + ccache -s + du -hs ~/.cache/ccache + save_pr_number: if: github.event_name == 'pull_request' runs-on: ubuntu-latest diff --git a/Src/Base/AMReX_PlotFileUtil.H b/Src/Base/AMReX_PlotFileUtil.H index 7e09d4b2e32..724891ac45e 100644 --- a/Src/Base/AMReX_PlotFileUtil.H +++ b/Src/Base/AMReX_PlotFileUtil.H @@ -12,6 +12,10 @@ #include #endif +#ifdef AMREX_USE_OPENPMD +#include +#endif + #include #include #include diff --git a/Src/CMakeLists.txt b/Src/CMakeLists.txt index febdea04b5b..e25a9cd4a64 100644 --- a/Src/CMakeLists.txt +++ b/Src/CMakeLists.txt @@ -200,6 +200,10 @@ if (AMReX_HDF5) add_subdirectory(Extern/HDF5) endif () +if (AMReX_OPENPMD) + add_subdirectory(Extern/openPMD-api) +endif() + # # Print out summary -- do it here so we already linked all # libs at this point diff --git a/Src/Extern/openPMD-api/AMReX_ParticlesOPENPMD.H b/Src/Extern/openPMD-api/AMReX_ParticlesOPENPMD.H new file mode 100644 index 00000000000..1cdba9c76b7 --- /dev/null +++ b/Src/Extern/openPMD-api/AMReX_ParticlesOPENPMD.H @@ -0,0 +1,78 @@ +#ifndef AMREX_PTL_OPENPMD_API_H +#define AMREX_PTL_OPENPMD_API_H + +/** @file This file is included from AMReX_ParticleCounter.H +* and belongs to the class ParticleContainer_impl. +* The namespace amrex was already declared. +* +* Therefore all the definitions in this file belongs to namespace amrex +*/ + +struct AMReX_PtlCounter +{ + int m_MPIRank = 0; + int m_MPISize = 1; + + unsigned long long m_Total = 0; + + std::vector m_ParticleCounterByLevel; + + [[nodiscard]] unsigned long long GetTotalNumParticles () const { return m_Total;} + + std::vector m_ParticleOffsetAtRank; + std::vector m_ParticleSizeAtRank; +}; + + +void CountParticles() +{ + + m_PtlCounter.m_MPISize = amrex::ParallelDescriptor::NProcs(); + m_PtlCounter.m_MPIRank = amrex::ParallelDescriptor::MyProc(); + + m_PtlCounter.m_ParticleCounterByLevel.resize(this->finestLevel()+1); + m_PtlCounter.m_ParticleOffsetAtRank.resize(this->finestLevel()+1); + m_PtlCounter.m_ParticleSizeAtRank.resize(this->finestLevel()+1); + + auto lf_GetParticleOffsetOfProcessor = [&](const long long& numParticles, + unsigned long long& offset, + unsigned long long& sum) -> void + { + std::vector result(m_PtlCounter.m_MPISize, 0); + amrex::ParallelGather::Gather (numParticles, result.data(), -1, amrex::ParallelDescriptor::Communicator()); + + sum = 0; + offset = 0; + for (int i=0; ifinestLevel(); currentLevel++) + { + // numParticles in this processor + auto numParticles = static_cast(this->NumberOfParticlesAtLevel(currentLevel, true, true)); + + unsigned long long offset=0; // offset of this level + unsigned long long sum=0; // numParticles in this level (sum from all processors) + lf_GetParticleOffsetOfProcessor(numParticles, offset, sum); + + m_PtlCounter.m_ParticleCounterByLevel[currentLevel] = sum; + m_PtlCounter.m_ParticleOffsetAtRank[currentLevel] = offset; + m_PtlCounter.m_ParticleSizeAtRank[currentLevel] = numParticles; + + // adjust offset, it should be numbered after particles from previous levels + for (auto lv=0; lv +#include +#include +#include +#include +#include + +#include +#include + +#ifdef AMREX_USE_EB +#include +#endif + +#include +#include +#include + +namespace amrex::openpmd_api { + + + //////////////////////////////////////// + // + // Struct AMReX_VarNameParser + // parser var names to field and comp names + // + //////////////////////////////////////// + + AMReX_VarNameParser::AMReX_VarNameParser(std::string const& varname) + :m_CompName(openPMD::MeshRecordComponent::SCALAR) + { + //auto [varname_no_mode, mode_index] = GetFieldNameModeInt(varname); + GetFieldNameModeInt(varname); + //bool var_in_theta_mode = mode_index != -1; // thetaMode or reconstructed Cartesian 2D slice + //std::string m_FieldName = varname_no_mode; + + // assume fields are scalar unless they match the following match of known vector fields + } + + void AMReX_VarNameParser::GetFieldNameModeInt (const std::string& varname) + { + // mode_index = -1 if varname isn't of form fieldName_mode_realOrImag + // mode_index = 2 * mode - 1 + (realOrImag == 'imag') + // in either case, there is a -1 in mode_index + //int mode_index = -1; + + std::regex e_real_imag("(.*)_([0-9]*)_(real|imag)"); + std::smatch sm; + std::regex_match(varname, sm, e_real_imag, std::regex_constants::match_default); + + if (sm.size() != 4 ) + { + m_ThetaMode = (m_ModeIndex != -1); // thetaMode or reconstructed Cartesian 2D slice + m_FieldName = varname; + } + else + { + // sm = [varname, field_name, mode, real_imag] + int mode = std::stoi(sm[2]); + if (mode == 0) { + m_ModeIndex = 0; + } else { + if (sm[3] == "imag") { + m_ModeIndex += 1; + } + m_ModeIndex += 2 * mode; + } + m_ThetaMode = (m_ModeIndex != -1); // thetaMode or reconstructed Cartesian 2D slice + m_FieldName = std::string(sm[1]); + } + } + + + void AMReX_VarNameParser::GetMeshCompNames (int meshLevel) + { + std::string varname = m_FieldName; + if (varname.size() >= 2U ) + { + std::string const varname_1st = varname.substr(0U, 1U); // 1st character + std::string const varname_2nd = varname.substr(1U, 1U); // 2nd character + + // Check if this field is a vector. If so, then extract the field name + + std::vector< std::string > const vector_fields = {"E", "B", "j"}; + std::vector< std::string > const field_components = getFieldComponentLabels(); + + for( std::string const& vector_field : vector_fields ) + { + for( std::string const& component : field_components ) + { + if( vector_field == varname_1st && component == varname_2nd ) + { + m_FieldName = varname_1st + varname.substr(2); // Strip component + m_CompName = varname_2nd; + } + } + } + } + + + if ( 0 == meshLevel ) { + return; + } + + m_FieldName += std::string("_lvl").append(std::to_string(meshLevel)); + } + + + + // create json options to pass to openpmd + inline std::string + getSeriesOptions (std::string const & operator_type, + std::map< std::string, std::string > const & operator_parameters, + std::string const & engine_type, + std::map< std::string, std::string > const & engine_parameters) + { + if (operator_type.empty() && engine_type.empty()) { + return "{}"; + } + + std::string options; + std::string top_block; + std::string end_block; + std::string op_block; + std::string en_block; + + std::string op_parameters; + for (const auto& kv : operator_parameters) { + if (!op_parameters.empty()) { op_parameters.append(",\n"); } + op_parameters.append(std::string(12, ' ')) /* just pretty alignment */ + .append("\"").append(kv.first).append("\": ") /* key */ + .append("\"").append(kv.second).append("\""); /* value (as string) */ + } + + std::string en_parameters; + for (const auto& kv : engine_parameters) { + if (!en_parameters.empty()) { en_parameters.append(",\n"); } + en_parameters.append(std::string(12, ' ')) /* just pretty alignment */ + .append("\"").append(kv.first).append("\": ") /* key */ + .append("\"").append(kv.second).append("\""); /* value (as string) */ + } + + top_block = R"END( +{ + "adios2": {)END"; + end_block = R"END( + } +})END"; + + if (!operator_type.empty()) { + op_block = R"END( + "dataset": { + "operators": [ + { + "type": ")END"; + op_block += operator_type + "\""; + if (!op_parameters.empty()) { + op_block += R"END(, + "parameters": { +)END"; + op_block += op_parameters + + "\n }"; + } + op_block += R"END( + } + ] + })END"; + + if (!engine_type.empty() || !en_parameters.empty()) { + op_block += ","; + } + } + + if (!engine_type.empty() || !en_parameters.empty()) + { + en_block = R"END( + "engine": {)END"; + if (!engine_type.empty()) { + en_block += R"END( + "type": ")END"; + en_block += engine_type + "\""; + if(!en_parameters.empty()) { + en_block += ","; + } + } + if (!en_parameters.empty()) { + en_block += R"END( + "parameters": { +)END"; + en_block += en_parameters + + "\n }"; + } + en_block += R"END( + })END"; + } + + options = top_block + op_block + en_block + end_block; + return options; + } + + + //////////////////////////////////////// + // + // Class AMReX_openPMDHandler + // files are saved as prefix/openpmd.bp + // + //////////////////////////////////////// + + AMReX_openPMDHandler::AMReX_openPMDHandler(std::string const& prefix) // NOLINT(modernize-pass-by-value) // match to diag_name in warpx + :m_Writer(nullptr) + { + BL_PROFILE("AMReX_openPMDHandler::()"); + CreateWriter(prefix); + } + + void AMReX_openPMDHandler::CreateWriter(const std::string& prefix) + { + ParmParse pp_prefix(prefix); + + // choose backend (e.g. ADIOS, ADIOS2 or HDF5). Default depends on openPMD-api configuration + std::string openpmd_backend {"default"}; + pp_prefix.query("openpmd_backend", openpmd_backend); + + std::string openpmd_encoding {"f"}; + pp_prefix.query("openpmd_encoding", openpmd_encoding); + openPMD::IterationEncoding encoding = openPMD::IterationEncoding::groupBased; + + if ( openpmd_encoding == "v" ) { + encoding = openPMD::IterationEncoding::variableBased; + } else if ( openpmd_encoding == "g" ) { + encoding = openPMD::IterationEncoding::groupBased; + } else if ( openpmd_encoding == "f" ) { + encoding = openPMD::IterationEncoding::fileBased; + } + + auto lf_collect = [&](const char* key, + const std::string& parameter_tag, + std::string& key_type, + std::map< std::string, std::string>& result)->void + { + //std::string key_type; + pp_prefix.query(key, key_type); + std::string const key_prefix = prefix + parameter_tag; + ParmParse pp; + auto entr = ParmParse::getEntries(key_prefix); + + auto const prefix_len = key_prefix.size() + 1; + for (std::string k : entr) { + std::string v; + pp.get(k.c_str(), v); + k.erase(0, prefix_len); + result.insert({k, v}); + } + }; + + std::string operator_type; + std::map< std::string, std::string > operator_parameters; + lf_collect("adios2_operator.type", ".adios2_operator.parameters", operator_type, operator_parameters); + + std::string engine_type; + std::map< std::string, std::string > engine_parameters; + lf_collect("adios2_engine.type", ".adios2_engine.parameters", engine_type, engine_parameters); + + std::string options=getSeriesOptions(operator_type, operator_parameters, + engine_type, engine_parameters); + + m_Writer = std::make_unique(prefix, encoding, openpmd_backend, options); + + pp_prefix.query("file_min_digits", m_Writer->m_openPMDMinDigits); + + } // CreateWriter() + + + void AMReX_openPMDHandler::SetWriter(amrex::openpmd_api::AMReX_openPMDWriter* w) + { + BL_ASSERT ( w != nullptr ); + + // assuer that input key/values are inherited + // so the openpmd filepath assigned from input file is still in use + w->m_openPMDPrefix = m_Writer->m_openPMDPrefix; + w->m_openPMDEncoding = m_Writer->m_openPMDEncoding; + w->m_openPMDFileType = m_Writer->m_openPMDFileType; + w->m_openPMDSeriesOptions = m_Writer->m_openPMDSeriesOptions; + + m_Writer.reset(w); + } + + //////////////////////////////////////// + // + // Class AMReX_openPMDWriter + // + //////////////////////////////////////// + + AMReX_openPMDWriter::AMReX_openPMDWriter (std::string prefix, + openPMD::IterationEncoding ie, + std::string filetype, + std::string options) + :m_openPMDPrefix(std::move(prefix)), + m_openPMDEncoding(ie), + m_openPMDFileType(std::move(filetype)), + m_openPMDSeriesOptions(std::move(options)) + //std::vector fieldPMLdirections // warpx specific + { + if( m_openPMDFileType == "default" ) { +#if openPMD_HAVE_ADIOS2==1 + m_openPMDFileType = "bp"; +#elif openPMD_HAVE_ADIOS1==1 + m_openPMDFileType = "bp"; +#elif openPMD_HAVE_HDF5==1 + m_openPMDFileType = "h5"; +#else + m_openPMDFileType = "json"; +#endif + } + } + + AMReX_openPMDWriter::~AMReX_openPMDWriter () + { + if( m_Series ) + { + m_Series->flush(); + m_Series.reset( nullptr ); + } + } + + void AMReX_openPMDWriter::SetStep(int ts) + { + m_CurrentStep = ts; + + Init(openPMD::Access::CREATE); + } + + void AMReX_openPMDWriter::CloseStep(int /*ts*/) + { + if (m_Series) { + GetIteration(m_CurrentStep).close(); + } + } + + void AMReX_openPMDWriter::Init(openPMD::Access access) + { + std::string filepath = m_openPMDPrefix; + GetFileName(filepath); + + if ( m_openPMDEncoding == openPMD::IterationEncoding::fileBased ) { + m_Series = nullptr; + } else if ( m_Series != nullptr ) { + return; + } + + if (amrex::ParallelDescriptor::NProcs() > 1) + { +#if defined(AMREX_USE_MPI) + m_Series = std::make_unique( + filepath, access, + amrex::ParallelDescriptor::Communicator(), + m_openPMDSeriesOptions + ); +#else + amrex::Abort(Utils::TextMsg::Err("AMReX did not build with MPI support!")); +#endif + } + else + { + m_Series = std::make_unique(filepath, access, m_openPMDSeriesOptions); + } + + m_Series->setIterationEncoding( m_openPMDEncoding ); + + m_Series->setMeshesPath( "fields" ); + // conform to ED-PIC extension of openPMD + + uint32_t const openPMD_ED_PIC = 1U; + m_Series->setOpenPMDextension( openPMD_ED_PIC ); + // meta info + + m_Series->setSoftware( "AMReX", amrex::Version() ); + } + + void AMReX_openPMDWriter::CompSetup(int lev, + openPMD::Container< openPMD::Mesh >& meshes, + amrex::Geometry& full_geom, + const std::vector& varnames, + const amrex::MultiFab* curr_mf) const + { + int const ncomp = curr_mf->nComp(); + for ( int icomp=0; icomp& meshes, + amrex::Geometry& full_geom, + const std::vector& varnames, + const amrex::MultiFab* curr_mf) + { + int const ncomp = curr_mf->nComp(); + amrex::Box const & global_box = full_geom.Domain(); + + + for ( int icomp=0; icompisManaged() || fab.arena()->isDevice()) + { + amrex::BaseFab foo(local_box, 1, amrex::The_Pinned_Arena()); + std::shared_ptr data_pinned(foo.release()); + amrex::Gpu::dtoh_memcpy_async(data_pinned.get(), fab.dataPtr(icomp), local_box.numPts()*sizeof(amrex::Real)); + // intentionally delayed until before we .flush(): amrex::Gpu::streamSynchronize(); + mesh_comp.storeChunk(data_pinned, chunk_offset, chunk_size); + } + else +#endif + { + amrex::Real const *local_data = fab.dataPtr(icomp); + mesh_comp.storeChunkRaw(local_data, chunk_offset, chunk_size); + } + } + } // icomp store loop + + } + + void AMReX_openPMDWriter::WriteMesh(const std::vector& varnames, + const amrex::Vector& mf, + const amrex::Vector& geom, + //const int iteration, + const double time ) const + + { + BL_PROFILE("AMReX_openPMDWriter::WriteMesh()"); + openPMD::Iteration series_iteration = GetIteration(m_CurrentStep); + series_iteration.open(); + + auto meshes = series_iteration.meshes; + series_iteration.setTime( time ); + + if ( varnames.empty() ) { return; } + + auto output_levels = int(geom.size()); + for (int lev=0; lev < output_levels; lev++) + { + amrex::Geometry full_geom = geom[lev]; + + if ( 0 == lev ) { + SetupFields(meshes, full_geom); + } + + CompSetup(lev, meshes, full_geom, varnames, mf[lev]); + CompStorage(lev, meshes, full_geom, varnames, mf[lev]); +#ifdef AMREX_USE_GPU + amrex::Gpu::streamSynchronize(); +#endif + m_Series->flush(); + } // for lev loop + } + + void AMReX_openPMDWriter::GetFileName(std::string& filepath) const + { + if (filepath.empty()) { + filepath.append("."); + } + + filepath.append("/"); + // transform paths for Windows +#ifdef _WIN32 + filepath = detail::replace_all(filepath, "/", "\\"); +#endif + + std::string filename = "openpmd"; + // + // OpenPMD supports timestepped names + // + if (m_openPMDEncoding == openPMD::IterationEncoding::fileBased) + { + std::string fileSuffix = std::string("_%0") + std::to_string(m_openPMDMinDigits) + std::string("T"); + filename = filename.append(fileSuffix); + } + filename.append(".").append(m_openPMDFileType); + filepath.append(filename); + } + + void AMReX_openPMDWriter::SetupFields (openPMD::Container< openPMD::Mesh >& meshes, + amrex::Geometry& full_geom) const + { + //} + // meta data for ED-PIC extension + auto const period = full_geom.periodicity(); + + std::vector fieldBoundary(6, "reflecting"); + std::vector particleBoundary(6, "absorbing"); + fieldBoundary.resize(AMREX_SPACEDIM * 2); + particleBoundary.resize(AMREX_SPACEDIM * 2); + +#if AMREX_SPACEDIM != 3 + fieldBoundary.resize(4); + particleBoundary.resize(4); +#endif + + for (int i = 0; i < int(fieldBoundary.size() / 2); ++i) { + if (period.isPeriodic(i)) { + fieldBoundary.at(2 * i) = "periodic"; + fieldBoundary.at(2 * i + 1) = "periodic"; + particleBoundary.at(2 * i) = "periodic"; + particleBoundary.at(2 * i + 1) = "periodic"; + } + } + + meshes.setAttribute("fieldBoundary", fieldBoundary); + meshes.setAttribute("particleBoundary", particleBoundary); + } + + + void AMReX_openPMDWriter::SetupMeshComp (openPMD::Mesh& mesh, + const amrex::Geometry& full_geom, + amrex::MultiFab const& mf, + const AMReX_VarNameParser& varName + ) const + { + BL_PROFILE("SetupMeshComp(default)"); + + auto mesh_comp = mesh[varName.m_CompName]; + amrex::Box const & global_box = full_geom.Domain(); + auto global_size = helper::getReversedVec(global_box.size() ); + + // - Grid spacing + std::vector const grid_spacing = helper::getReversedVec(full_geom.CellSize()); + mesh.setGridSpacing(grid_spacing); + + // - Global offset + std::vector const global_offset = helper::getReversedVec(full_geom.ProbLo()); + mesh.setGridGlobalOffset(global_offset); + + // - AxisLabels + std::vector axis_labels = varName.getFieldAxisLabels(); + mesh.setAxisLabels(axis_labels); + + // Prepare the type of dataset that will be written + openPMD::Datatype const datatype = openPMD::determineDatatype(); + auto const dataset = openPMD::Dataset(datatype, global_size); + mesh.setDataOrder(openPMD::Mesh::DataOrder::C); + + if (varName.m_ThetaMode) { + mesh.setGeometry("thetaMode"); + } + + mesh.setAttribute("fieldSmoothing", "none"); + mesh_comp.resetDataset(dataset); + + auto relative_cell_pos = helper::getRelativeCellPosition(mf); // AMReX Fortran index order + std::reverse( relative_cell_pos.begin(), relative_cell_pos.end() ); // now in C order + mesh_comp.setPosition( relative_cell_pos ); + } + + +} // namespace amrex::openpmd_api diff --git a/Src/Extern/openPMD-api/AMReX_PlotFileOPENPMD_PTL.cpp b/Src/Extern/openPMD-api/AMReX_PlotFileOPENPMD_PTL.cpp new file mode 100644 index 00000000000..e5e70ef8203 --- /dev/null +++ b/Src/Extern/openPMD-api/AMReX_PlotFileOPENPMD_PTL.cpp @@ -0,0 +1,102 @@ +#include +#include +#include +#include +#include +#include + +#ifdef AMREX_USE_EB +#include +#endif + +#include +#include + +#include +#include +#include + + +namespace amrex::openpmd_api { + + bool AMReX_openPMDWriter::AllocatePtlProperties(openPMD::ParticleSpecies& currSpecies, + const amrex::Vector& write_real_comp, + const amrex::Vector& real_comp_names, + const amrex::Vector& write_int_comp, + const amrex::Vector& int_comp_names, + const unsigned long long np) const + { + SetupPos(currSpecies, np); + + // Allocate _all_ datasets of dtype. + // handle scalar and non-scalar records by name + auto const lf_compRecordInit = [&currSpecies](const amrex::Vector& write_comp, + const amrex::Vector& comp_names, + openPMD::Dataset& dtype) + { + auto const min_counter = std::min(write_comp.size(), comp_names.size()); + for (int i = 0; i < min_counter; ++i) + { + if (write_comp[i]) { + helper::getComponentRecord(currSpecies, comp_names[i]).resetDataset(dtype); + } + } + }; + auto dtype_real = openPMD::Dataset(openPMD::determineDatatype(), {np}, m_openPMDDatasetOptions); + lf_compRecordInit(write_real_comp, real_comp_names, dtype_real); + + auto dtype_int = openPMD::Dataset(openPMD::determineDatatype(), {np}, m_openPMDDatasetOptions); + lf_compRecordInit(write_int_comp, int_comp_names, dtype_int); + + return true; + } + + void AMReX_openPMDWriter::SetupPos(openPMD::ParticleSpecies& currSpecies, + const unsigned long long& np) const + { + auto realType = openPMD::Dataset(openPMD::determineDatatype(), {np}, m_openPMDDatasetOptions); + auto idType = openPMD::Dataset(openPMD::determineDatatype< uint64_t >(), {np}, m_openPMDDatasetOptions); + + auto const positionComponents = /*helper::*/getParticlePositionComponentLabels(); + for( auto const& comp : positionComponents ) + { + currSpecies["position"][comp].resetDataset( realType ); + } + + auto const * const scalar = openPMD::RecordComponent::SCALAR; + currSpecies["id"][scalar].resetDataset( idType ); + } + + + unsigned long long AMReX_openPMDWriter::GetGrandOffset() const + { + return 0; + } + + // from warpx SetConstParticleRecordsEDPIC () + // + // constant values, just call before flushing + // + void AMReX_openPMDWriter::SetupConstant(openPMD::ParticleSpecies& currSpecies, + const unsigned long long& np) const + { + auto realType = openPMD::Dataset(openPMD::determineDatatype(), {np}, m_openPMDDatasetOptions); + + auto const positionComponents = getParticlePositionComponentLabels(); + for( auto const& comp : positionComponents ) { + currSpecies["positionOffset"][comp].resetDataset( realType ); + } + + // make constant + using namespace amrex::literals; + for( auto const& comp : positionComponents ) { + currSpecies["positionOffset"][comp].makeConstant( 0._prt ); + } + + currSpecies["positionOffset"].setUnitDimension( helper::getUnitDimension("positionOffset") ); + currSpecies["position"].setUnitDimension( helper::getUnitDimension("position") ); + + SetParticleSpecieAttributes(currSpecies); + } + +} diff --git a/Src/Extern/openPMD-api/AMReX_PlotFileUtilOPENPMD.H b/Src/Extern/openPMD-api/AMReX_PlotFileUtilOPENPMD.H new file mode 100644 index 00000000000..f14175dda49 --- /dev/null +++ b/Src/Extern/openPMD-api/AMReX_PlotFileUtilOPENPMD.H @@ -0,0 +1,403 @@ +#ifndef AMREX_PLOTFILE_OPENPMD_API_H +#define AMREX_PLOTFILE_OPENPMD_API_H + +#include + +#include +#include +#include + +#ifdef AMREX_PARTICLES +#include +#include +#endif + +#ifdef AMREX_USE_OPENPMD +#include +#endif + +#include +#include + +namespace amrex::openpmd_api +{ + namespace helper + { + // from warpx Utils/RelativeCellPosition.cpp + inline std::vector< double > + getRelativeCellPosition(amrex::MultiFab const& mf) + { + amrex::IndexType const idx_type = mf.ixType(); + + std::vector< double > relative_position(AMREX_SPACEDIM, 0.0); + // amrex::CellIndex::CELL means: 0.5 from lower corner for that index/direction + // amrex::CellIndex::NODE means: at corner for that index/direction + // WarpX::do_nodal means: all indices/directions on CellIndex::NODE + for (int d = 0; d < AMREX_SPACEDIM; d++) + { + if (idx_type.cellCentered(d)) { + relative_position.at(d) = 0.5; + } + } + return relative_position; + } + + inline std::vector + getReversedVec( const IntVect& v ) + { + // Convert the IntVect v to and std::vector u + std::vector u = { + AMREX_D_DECL( + static_cast(v[0]), + static_cast(v[1]), + static_cast(v[2]) + ) + }; + // Reverse the order of elements, if v corresponds to the indices of a + // Fortran-order array (like an AMReX FArrayBox) + // but u is intended to be used with a C-order API (like openPMD) + std::reverse( u.begin(), u.end() ); + + return u; + } + + inline std::vector + getReversedVec( const Real* v ) + { + // Convert Real* v to and std::vector u + std::vector u = { + AMREX_D_DECL( + static_cast(v[0]), + static_cast(v[1]), + static_cast(v[2]) + ) + }; + // Reverse the order of elements, if v corresponds to the indices of a + // Fortran-order array (like an AMReX FArrayBox) + // but u is intended to be used with a C-order API (like openPMD) + + std::reverse( u.begin(), u.end() ); + + return u; + } + } + + struct AMReX_VarNameParser; // forward declaration + class AMReX_openPMDWriter; + class AMReX_openPMDHandler; + + /* + */ + + void UseCustomWriter(AMReX_openPMDWriter* w); + + std::unique_ptr InitUserHandler(const std::string& prefix); + void CloseUserHandler(std::unique_ptr& userHandler); + + void InitHandler(const std::string& prefix); + //void SetCustomizer(AMReX_Optional* f); + + void SetStep(int ts); + void CloseStep(int ts); + void CloseHandler(); + + void WriteSingleLevel (//const std::string &plotfilename, + const MultiFab &mf, + const Vector &varnames, + const Geometry &geom, + Real t, + //int level_step, + const std::string &ignored_versionName = "HyperCLaw-V1.1", + const std::string &ignored_levelPrefix = "Level_", + const std::string &ignored_mfPrefix = "Cell", + const Vector& ignored_extra_dirs = Vector()); + + void WriteMultiLevel ( + //int nlevels, will all levels in mf & geom + const Vector &mf, + const Vector &varnames, + const Vector &geom, + Real time, + //const Vector &level_steps, + const Vector &ref_ratio, + const std::string &ignored_versionName = "HyperCLaw-V1.1", + const std::string &ignored_levelPrefix = "Level_", + const std::string &ignored_mfPrefix = "Cell", + const Vector& ignored_extra_dirs = Vector()); + + + //////////////////////////// + struct AMReX_VarNameParser + { + AMReX_VarNameParser(std::string const& varname); + virtual ~AMReX_VarNameParser () = default; + AMReX_VarNameParser (AMReX_VarNameParser const&) = default; + AMReX_VarNameParser (AMReX_VarNameParser &&) = delete; + AMReX_VarNameParser& operator= (AMReX_VarNameParser const&) = default; + AMReX_VarNameParser& operator= (AMReX_VarNameParser &&) = delete; + + void GetFieldNameModeInt(const std::string& varname); + void GetMeshCompNames (int meshLevel); + + [[nodiscard]] std::vector< std::string > getFieldComponentLabels () const + { + using vs = std::vector< std::string >; + if (m_ThetaMode) + { + // if we write individual modes + vs fieldComponents{"r", "t", "z"}; + return fieldComponents; + } + else + { + // if we just write reconstructed fields at theta=0 or are Cartesian + // note: 1D3V and 2D3V simulations still have 3 components for the fields + vs fieldComponents{"x", "y", "z"}; + return fieldComponents; + } + } + + [[nodiscard]] virtual std::vector< std::string > + getFieldAxisLabels () const + { + using vs = std::vector< std::string >; + + // temporary resolution +#if AMREX_SPACEDIM == 1 + vs const axisLabels{"z"}; +#elif AMREX_SPACEDIM == 2 + vs const axisLabels{"x", "z"}; +#elif AMREX_SPACEDIM == 3 + vs const axisLabels{"x", "y", "z"}; // x varies fastest in memory +#else +# error Unable to label more than 3d + // no labels to be addressed + vs const axisLabels{} +#endif + // revert to C order (fastest varying index last) + return {axisLabels.rbegin(), axisLabels.rend()}; + } // getFieldAxisLabels + + + ///// members + std::string m_FieldName; + std::string m_CompName; + bool m_ThetaMode = false; + int m_ModeIndex=-1; + + }; + + + //////////////////////////// + //////////////////////////// + class AMReX_openPMDWriter { + public: + AMReX_openPMDWriter(std::string prefix, + openPMD::IterationEncoding ie, + std::string filetype, + std::string openpmdOptions); + + AMReX_openPMDWriter() = default; + AMReX_openPMDWriter (AMReX_openPMDWriter const&) = delete; + AMReX_openPMDWriter (AMReX_openPMDWriter &&) = delete; + AMReX_openPMDWriter& operator= (AMReX_openPMDWriter const&) = delete; + AMReX_openPMDWriter& operator= (AMReX_openPMDWriter &&) = delete; + + virtual ~AMReX_openPMDWriter(); + + [[nodiscard]] virtual openPMD::Iteration GetIteration (int const iteration) const + { + return m_Series->writeIterations()[iteration]; + } + + virtual void SetStep(int ts); + virtual void CloseStep(int ts); + virtual void Init(openPMD::Access access); + + virtual void WriteMesh(const std::vector& varnames, + const amrex::Vector& mf, + const amrex::Vector& geom, + //const Vector &iteration, + //const int iteration, /* note: all levels are outputting the same step */ + double time ) const; + + template + void DumpParticles(PC const& pc, + const std::string& name, + //const int iteration, + const amrex::Vector& write_real_comp, // size = NStructReal + NArrayReal + const amrex::Vector& write_int_comp, // size = NStructInt + NArrayInt + const amrex::Vector& real_comp_names, + const amrex::Vector& int_comp_names + ) const; + + template + void DumpParticles(PC const& pc, + const std::string& name, + //const int iteration, + const amrex::Vector& write_real_comp, // size = NStructReal + NArrayReal + const amrex::Vector& write_int_comp, // size = NStructInt + NArrayInt + const amrex::Vector& real_comp_names, + const amrex::Vector& int_comp_names, + FUNC_pc&& pc_func, + FUNC_pit&& posId_func + ) const; + + [[nodiscard]] virtual bool AllocatePtlProperties(openPMD::ParticleSpecies& currSpecies, + const amrex::Vector& write_real_comp, + const amrex::Vector& real_comp_names, + const amrex::Vector& write_int_comp, + const amrex::Vector& int_comp_names, + unsigned long long np) const; + + virtual void SetParticleSpecieAttributes(openPMD::ParticleSpecies& /*currSpecies*/) const + {} + + void SetupPos(openPMD::ParticleSpecies& species, const unsigned long long& np) const; + [[nodiscard]] virtual unsigned long long GetGrandOffset() const; // to adjust offset/total with data already flushed, example use case: BTD + + // called once because it is constant valued record + virtual void SetupConstant(openPMD::ParticleSpecies& currSpecies, + const unsigned long long& np) const; + + [[nodiscard]] virtual std::vector< std::string > + getParticlePositionComponentLabels() const + { +#if (AMREX_SPACEDIM == 3) + std::vector< std::string > positionComponents{"x", "y", "z"}; +#elif (AMREX_SPACEDIM == 2) + std::vector< std::string > positionComponents{"x", "y"}; +#elif (AMREX_SPACEDIM == 1) + std::vector< std::string > positionComponents{"x"}; +#else +# error Unknown dimensionality. +#endif + return positionComponents; + } + + template + void + SetupRealProperties(PC const& pc, + openPMD::ParticleSpecies& currSpecies, + const amrex::Vector& write_real_comp, + const amrex::Vector& real_comp_names, + const amrex::Vector& write_int_comp, + const amrex::Vector& int_comp_names, + unsigned long long np) const; + /* + */ + template + void + SaveRealProperty (PIt& pti, + openPMD::ParticleSpecies& currSpecies, + unsigned long long offset, + const amrex::Vector& write_real_comp, + const amrex::Vector& real_comp_names, + const amrex::Vector& write_int_comp, + const amrex::Vector& int_comp_names) const; + + + std::unique_ptr m_Series = nullptr; + std::string m_openPMDPrefix = std::string(); + + int m_openPMDMinDigits = 6; + std::string m_openPMDDatasetOptions = "{}"; + + openPMD::IterationEncoding m_openPMDEncoding = openPMD::IterationEncoding::fileBased; + + std::string m_openPMDFileType = "bp"; + std::string m_openPMDSeriesOptions = std::string(); + + protected: + void CompSetup(int lev, + openPMD::Container< openPMD::Mesh >& meshes, + amrex::Geometry& full_geom, + const std::vector& varnames, + const amrex::MultiFab* mf) const; + + static void CompStorage(int lev, + openPMD::Container< openPMD::Mesh >& meshes, + amrex::Geometry& full_geom, + const std::vector& varnames, + const amrex::MultiFab* mf); + + + virtual void SetupFields (openPMD::Container< openPMD::Mesh >& meshes, + amrex::Geometry& full_geom) const; + + virtual void SetupMeshComp (openPMD::Mesh& mesh, + const amrex::Geometry& full_geom, + amrex::MultiFab const& mf, + const AMReX_VarNameParser& varName) const; + + void GetFileName(std::string& filepath) const; + + template + void StoreAoS_Real(PIt& pti, + openPMD::ParticleSpecies& currSpecies, + const amrex::Vector& write_real_comp, + const amrex::Vector& real_comp_names, + unsigned long long offset) const; + + template + void + SavePosId(PIt& pti, + openPMD::ParticleSpecies& currSpecies, + unsigned long long offset) const; + + template + void + SaveSpecieId(PIt& pti, + openPMD::ParticleSpecies& currSpecies, + unsigned long long offset) const; + + template + void + StoreAoS_Int(PIt& pti, + openPMD::ParticleSpecies& currSpecies, + const amrex::Vector& write_int_comp, + const amrex::Vector& int_comp_names, + unsigned long long offset) const; + + + template + void StoreSoAReal(PIt& pti, + openPMD::ParticleSpecies& currSpecies, + const amrex::Vector& write_real_comp, + const amrex::Vector& real_comp_names, + unsigned long long offset) const; + + + template + void StoreSoAInt(PIt& pti, + openPMD::ParticleSpecies& currSpecies, + const amrex::Vector& write_int_comp, + const amrex::Vector& int_comp_names, + unsigned long long offset) const; + + int m_CurrentStep = -1; + + }; + + + + //////////////////////////// + class AMReX_openPMDHandler { + public: + AMReX_openPMDHandler(std::string const& prefix = std::string()); + + std::unique_ptr m_Writer; + + void SetWriter(amrex::openpmd_api::AMReX_openPMDWriter* w); + private: + void CreateWriter(const std::string& prefix = std::string()); + + }; // class AMReX_openPMDHandler + + +} // namespace amrex::openpmd_api + + +#include + +#endif // AMREX_PLOTFILE_OPENPMD_API_H diff --git a/Src/Extern/openPMD-api/AMReX_PlotFileUtilOPENPMD.cpp b/Src/Extern/openPMD-api/AMReX_PlotFileUtilOPENPMD.cpp new file mode 100644 index 00000000000..0e020b9ae0f --- /dev/null +++ b/Src/Extern/openPMD-api/AMReX_PlotFileUtilOPENPMD.cpp @@ -0,0 +1,140 @@ +#include +#include +#include +#include +#include +#include + +#ifdef AMREX_USE_EB +#include +#endif + +#include + +#include +#include + +namespace amrex::openpmd_api +{ + /* global handler, activate with InitHandler() & deactivate with CloseHandler() */ + std::unique_ptr< AMReX_openPMDHandler > m_OpenPMDHandler = nullptr; + + std::unique_ptr InitUserHandler(const std::string& prefix) + { + std::string filePath; + if (prefix.empty()) + { + ParmParse pp; + pp.query("openpmd_directory", filePath); + } + else { + filePath = prefix; + } + + return std::make_unique(filePath); + } + + void CloseUserHandler(std::unique_ptr& userHandler) + { + if (userHandler == nullptr) { return; } + + userHandler.reset(nullptr); + } + + void InitHandler(const std::string& prefix) + { + std::string filePath; + if (prefix.empty()) + { + ParmParse pp; + pp.query("openpmd_directory", filePath); + } + else { + filePath = prefix; + } + + if (m_OpenPMDHandler == nullptr) { + m_OpenPMDHandler = std::make_unique(filePath); + } else if (m_OpenPMDHandler->m_Writer != nullptr) { + if (m_OpenPMDHandler->m_Writer->m_openPMDPrefix != filePath) { + m_OpenPMDHandler = std::make_unique(filePath); + } + } + // already using the directory, no action needed + } + + void UseCustomWriter(AMReX_openPMDWriter* w) + { + BL_ASSERT ( m_OpenPMDHandler != nullptr ); + BL_ASSERT ( w != nullptr ); + + m_OpenPMDHandler->SetWriter(w); + } + + void CloseHandler() + { + m_OpenPMDHandler.reset(); + } + + void SetStep(int ts) + { + if ((m_OpenPMDHandler == nullptr) || (m_OpenPMDHandler->m_Writer == nullptr)) { + return; + } + + m_OpenPMDHandler->m_Writer->SetStep(ts); + } + + void CloseStep(int ts) + { + if ((m_OpenPMDHandler == nullptr) || (m_OpenPMDHandler->m_Writer == nullptr)) { + return; + } + + m_OpenPMDHandler->m_Writer->CloseStep(ts); + } + + void WriteSingleLevel (//const std::string &plotfilename, + const MultiFab &mf, + const Vector &varnames, + const Geometry &geom, + Real t, + //int level_step, + const std::string &versionName, + const std::string &levelPrefix, + const std::string &mfPrefix, + const Vector& extra_dirs) + { + Vector v_mf(1,&mf); + Vector v_geom(1,geom); + //Vector v_level_steps(1,level_step); + Vector ref_ratio; + + WriteMultiLevel(v_mf, varnames, v_geom, t, /*v_level_steps,*/ ref_ratio, versionName, levelPrefix, mfPrefix, extra_dirs); + } + + void WriteMultiLevel ( + //int nlevels, // will write all levels in mf & geom + const Vector &mf, + const Vector &varnames, + const Vector &geom, + Real time, + //const Vector &level_steps, + const Vector & /*ref_ratio*/, + const std::string & /*versionName*/, + const std::string & /*levelPrefix*/, + const std::string & /*mfPrefix*/, + const Vector& /*extra_dirs*/) + { + if ((m_OpenPMDHandler == nullptr) || (m_OpenPMDHandler->m_Writer == nullptr)) { return; } + + BL_ASSERT ( geom.size() == mf.size() ); + BL_ASSERT ( mf[0]->nComp() <= varnames.size() ); + + m_OpenPMDHandler->m_Writer->WriteMesh(varnames, + mf, //amrex::GetVecOfConstPtrs(mf), + geom, + //level_steps[0], + time); + } +} // namespace amrex::openpmd_api diff --git a/Src/Extern/openPMD-api/AMReX_PlotFileUtilOPENPMD_PTLImpl.H b/Src/Extern/openPMD-api/AMReX_PlotFileUtilOPENPMD_PTLImpl.H new file mode 100644 index 00000000000..f7cbe9a2774 --- /dev/null +++ b/Src/Extern/openPMD-api/AMReX_PlotFileUtilOPENPMD_PTLImpl.H @@ -0,0 +1,657 @@ +#ifndef AMREX_PLOTFILE_OPENPMD_PTLS_API_H +#define AMREX_PLOTFILE_OPENPMD_PTLS_API_H + +namespace amrex::openpmd_api +{ + extern std::unique_ptr< amrex::openpmd_api::AMReX_openPMDHandler > m_OpenPMDHandler; + //////////////////////////////////////////////////////////////////////// + // + // Global function to write particle container + // each WriteParticles() call has two versions: + // with FUNC(for user supploed saveposid) or not (using default from the writer) + // + //////////////////////////////////////////////////////////////////////// + + template + void WriteParticles(PC const& pc, + const std::string& specieName, + //int iteration, + const Vector& real_comp_names, + const Vector& int_comp_names, + const Vector& do_write_real_comp, + const Vector& do_write_int_comp, + FUNC_pc&& pcFunc, + FUNC_pit&& posIdFunc) + { + BL_PROFILE("amrex::openpmd_api::WriteParticles()"); + + if ((m_OpenPMDHandler == nullptr) || (m_OpenPMDHandler->m_Writer == nullptr)) { + return; + } + + AMREX_ASSERT( do_write_real_comp.size() <= real_comp_names.size() ); + AMREX_ASSERT( do_write_int_comp.size() <= int_comp_names.size() ); + AMREX_ASSERT( real_comp_names.size() <= pc.NStructReal + pc.NumRealComps() ); + + m_OpenPMDHandler->m_Writer->DumpParticles(pc, specieName, //iteration, + do_write_real_comp, do_write_int_comp, real_comp_names, int_comp_names, + pcFunc, posIdFunc); + } + + // use default pos_id func + template + void WriteParticles(PC const& pc, + const std::string& specieName, + //int iteration, + const Vector& real_comp_names, + const Vector& int_comp_names, + const Vector& do_write_real_comp, + const Vector& do_write_int_comp) + { + BL_PROFILE("amrex::openpmd_api::WriteParticles(..)"); + + amrex::Print()<<" Real: "<m_Writer == nullptr)) { + return; + } + + AMREX_ASSERT( do_write_real_comp.size() <= real_comp_names.size() ); + AMREX_ASSERT( do_write_int_comp.size() <= int_comp_names.size() ); + AMREX_ASSERT( real_comp_names.size() <= pc.NStructReal + pc.NumRealComps() ); + + m_OpenPMDHandler->m_Writer->DumpParticles(pc, specieName, //iteration, + do_write_real_comp, do_write_int_comp, real_comp_names, int_comp_names); + } + + template + void WriteParticles(PC const& pc, + const std::string& specieName, + //int iteration, + const Vector& real_comp_names, + const Vector& int_comp_names) + { + Vector write_real_comp; + + for (int i = 0; i < real_comp_names.size(); ++i ) + { + write_real_comp.push_back(1); + } + + Vector write_int_comp; + for (int i = 0; i < int_comp_names.size(); ++i ) + { + write_int_comp.push_back(1); + } + + WriteParticles(pc, specieName, //iteration, + real_comp_names, int_comp_names, + write_real_comp, write_int_comp); + } + + + template + void WriteParticles(PC const& pc, + const std::string& specieName) + { + Vector write_real_comp; + Vector real_comp_names; + + for (int i = 0; i < pc.NStructReal + pc.NumRealComps(); ++i ) + { + write_real_comp.push_back(1); + + std::stringstream ss; + if (i < pc.NStructReal) { + ss << "aos-r_" << i; + } else { + ss << "SoA-r" << i; + } + + real_comp_names.push_back(ss.str()); + } + + Vector write_int_comp; + Vector int_comp_names; + for (int i = 0; i < pc.NStructInt + pc.NumIntComps(); ++i ) + { + write_int_comp.push_back(1); + std::stringstream ss; + + if (i < pc.NStructInt) { + ss << "aos-i_" << i; + } else { + ss << "SoA-i" << i; + } + + int_comp_names.push_back(ss.str()); + } + + WriteParticles(pc, specieName, //iteration, + real_comp_names, int_comp_names, + write_real_comp, write_int_comp); + } + + + template + void WriteParticles(PC const& pc, + const std::string& specieName, + FUNC_pc&& pcFunc, + FUNC_pit&& pitFunc + ) + { + Vector write_real_comp; + Vector real_comp_names; + + for (int i = 0; i < pc.NStructReal + pc.NumRealComps(); ++i ) + { + write_real_comp.push_back(1); + + std::stringstream ss; + if (i < pc.NStructReal) { + ss << "aos-r_" << i; + } else { + ss << "SoA-r" << i; + } + + real_comp_names.push_back(ss.str()); + } + + Vector write_int_comp; + Vector int_comp_names; + for (int i = 0; i < pc.NStructInt + pc.NumIntComps(); ++i ) + { + write_int_comp.push_back(1); + std::stringstream ss; + + if (i < pc.NStructInt) { + ss << "aos-i_" << i; + } else { + ss << "SoA-i" << i; + } + int_comp_names.push_back(ss.str()); + } + + WriteParticles(pc, specieName, //iteration, + real_comp_names, int_comp_names, + write_real_comp, write_int_comp, pcFunc, pitFunc); + } + + namespace helper { + //////////////////////////////////// + // + // all helper functions + // + //////////////////////////////////// + inline std::map< openPMD::UnitDimension, double > + getUnitDimension ( std::string const & record_name ) + { + + if( record_name == "position" ) { + return {{openPMD::UnitDimension::L, 1.}}; + } else if( record_name == "positionOffset" ) { + return {{openPMD::UnitDimension::L, 1.}}; + } else if( record_name == "momentum" ) { + return {{openPMD::UnitDimension::L, 1.}, + {openPMD::UnitDimension::M, 1.}, + {openPMD::UnitDimension::T, -1.}}; + } else if( record_name == "charge" ) { + return {{openPMD::UnitDimension::T, 1.}, + {openPMD::UnitDimension::I, 1.}}; + } else if( record_name == "mass" ) { + return {{openPMD::UnitDimension::M, 1.}}; + } else if( record_name == "E" ) { + return {{openPMD::UnitDimension::L, 1.}, + {openPMD::UnitDimension::M, 1.}, + {openPMD::UnitDimension::T, -3.}, + {openPMD::UnitDimension::I, -1.}}; + } else if( record_name == "B" ) { + return {{openPMD::UnitDimension::M, 1.}, + {openPMD::UnitDimension::I, -1.}, + {openPMD::UnitDimension::T, -2.}}; + } else { + return {}; + } + } + + // + // + // + constexpr uint64_t + localIDtoGlobal(int const id, int const cpu) + { + // from WarpXUtilIO::localIDtoGlobal + static_assert(sizeof(int) * 2U <= sizeof(uint64_t), + "int size might cause collisions in global IDs"); + + return uint64_t(id) | uint64_t(cpu) << 32U; + } + + + inline std::pair< std::string, std::string > + name2openPMD ( std::string const& fullName ) + { + std::string record_name = fullName; + std::string component_name = openPMD::RecordComponent::SCALAR; + + // we use "_" as separator in names to group vector records + std::size_t startComp = fullName.find_last_of('_'); + if( startComp != std::string::npos ) { // non-scalar + record_name = fullName.substr(0, startComp); + component_name = fullName.substr(startComp + 1U); + } + return make_pair(record_name, component_name); + } + + // + // used by both setup/saverealproperty so better be here than a lambda function + // + inline + auto getComponentRecord (openPMD::ParticleSpecies& currSpecies, + std::string const& comp_name) + { + // handle scalar and non-scalar records by name + const auto [record_name, component_name] = name2openPMD(comp_name); + return currSpecies[record_name][component_name]; + } + + } // namespace helper + + + + + + + template + void + AMReX_openPMDWriter::StoreAoS_Real(PIt& pti, + openPMD::ParticleSpecies& currSpecies, + const amrex::Vector& write_real_comp, + const amrex::Vector& real_comp_names, + unsigned long long offset) const + { + auto const& aos = pti.GetArrayOfStructs(); // size = numParticlesOnTile + auto const numParticleOnTile = pti.numParticles(); + auto const numParticleOnTile64 = static_cast( numParticleOnTile ); + + // first we concatenate the AoS into contiguous arrays + { + for( auto idx=0; idx d( + new amrex::ParticleReal[numParticleOnTile], + [](amrex::ParticleReal const *p){ delete[] p; } + ); + + for( auto kk=0; kk + void + AMReX_openPMDWriter::StoreAoS_Int(PIt& pti, + openPMD::ParticleSpecies& currSpecies, + const amrex::Vector& write_int_comp, + const amrex::Vector& int_comp_names, + unsigned long long offset) const + { + auto const& aos = pti.GetArrayOfStructs(); // size = numParticlesOnTile + auto const numParticleOnTile = pti.numParticles(); + auto const numParticleOnTile64 = static_cast( numParticleOnTile ); + + // first we concatenate the AoS into contiguous arrays + { + for( auto idx=0; idx d( + new int[numParticleOnTile], + [](int const *p){ delete[] p; } + ); + + for( auto kk=0; kk + void + AMReX_openPMDWriter::StoreSoAReal(PIt& pti, + openPMD::ParticleSpecies& currSpecies, + const amrex::Vector& write_real_comp, + const amrex::Vector& real_comp_names, + unsigned long long offset) const + { + auto const& soa = pti.GetStructOfArrays(); + auto const numParticleOnTile = pti.numParticles(); + auto const numParticleOnTile64 = static_cast( numParticleOnTile ); + + auto const real_counter = std::min(write_real_comp.size(), real_comp_names.size()) - PIt::ContainerType::NStructReal; + + for (auto idx=0; idx + void AMReX_openPMDWriter::StoreSoAInt(PIt& pti, + openPMD::ParticleSpecies& currSpecies, + const amrex::Vector& write_int_comp, + const amrex::Vector& int_comp_names, + unsigned long long offset) const + { + auto const& soa = pti.GetStructOfArrays(); + auto const numParticleOnTile = pti.numParticles(); + auto const numParticleOnTile64 = static_cast( numParticleOnTile ); + + auto const int_counter = std::min(write_int_comp.size(), int_comp_names.size()) - PIt::ContainerType::NStructInt; + for (auto idx=0; idx + void + AMReX_openPMDWriter::SavePosId(PIt& pti, + openPMD::ParticleSpecies& currSpecies, + unsigned long long offset) const + { + BL_PROFILE("amrex::openpmd_api::SavePosId(default..)"); + auto const numParticleOnTile = pti.numParticles(); + auto const numParticleOnTile64 = static_cast( numParticleOnTile ); + + // get position and particle ID from aos + // note: this implementation iterates the AoS 4x... + // if we flush late as we do now, we can also copy out the data in one go + const auto& aos = pti.GetArrayOfStructs(); // size = numParticlesOnTile + { + // Save positions + auto const positionComponents = /*helper::*/getParticlePositionComponentLabels(); + + for (auto currDim = 0; currDim < AMREX_SPACEDIM; currDim++) + { + //amrex::AllPrint()< curr(new amrex::ParticleReal[numParticleOnTile], + [](amrex::ParticleReal const *p) { delete[] p; } + ); + for (auto i = 0; i < numParticleOnTile; i++) { + //amrex::AllPrint()<<" i = "< + void + AMReX_openPMDWriter::SaveSpecieId(PIt& pti, + openPMD::ParticleSpecies& currSpecies, + unsigned long long offset) const + { + BL_PROFILE("amrex::openpmd_api::SaveSpecieId"); + auto const numParticleOnTile = pti.numParticles(); + auto const numParticleOnTile64 = static_cast( numParticleOnTile ); + + // get position and particle ID from aos + // note: this implementation iterates the AoS 4x... + // if we flush late as we do now, we can also copy out the data in one go + const auto& aos = pti.GetArrayOfStructs(); // size = numParticlesOnTile + { + // save particle ID after converting it to a globally unique ID + std::shared_ptr< uint64_t > ids(new uint64_t[numParticleOnTile], + [](uint64_t const *p){ delete[] p; } + ); + for (auto i=0; i + void + AMReX_openPMDWriter::SaveRealProperty (PIt& pti, + openPMD::ParticleSpecies& currSpecies, + unsigned long long offset, + const amrex::Vector& write_real_comp, + const amrex::Vector& real_comp_names, + const amrex::Vector& write_int_comp, + const amrex::Vector& int_comp_names) const + { + // note: WarpX does not yet use extra AoS Real attributes + StoreAoS_Real(pti, currSpecies, write_real_comp, real_comp_names, offset); + StoreAoS_Int(pti, currSpecies, write_int_comp, int_comp_names, offset); + + //auto const& soa = pti.GetStructOfArrays(); + StoreSoAReal(pti, currSpecies, write_real_comp, real_comp_names, offset); + StoreSoAInt(pti, currSpecies, write_int_comp, int_comp_names, offset); + } + + + + + //////////////////////////////////// + // + // AMReX_openPMDWriter::DumpParticles() + // + //////////////////////////////////// + + template + void AMReX_openPMDWriter::DumpParticles(PC const& pc, + const std::string& name, + //int iteration, + const amrex::Vector& write_real_comp, + const amrex::Vector& write_int_comp, + const amrex::Vector& real_comp_names, + const amrex::Vector& int_comp_names + ) const + { + DumpParticles(pc, name, write_real_comp, write_int_comp, real_comp_names, int_comp_names, + [=] (auto& /*pc*/, openPMD::ParticleSpecies& /*currSpecies*/, unsigned long long /*size*/) + { + }, // no extra pc level info to save + [=] (auto& pti, openPMD::ParticleSpecies& currSpecies, unsigned long long offset) + { + SavePosId(pti, currSpecies, offset); // default way of save pos & id + }); + } + + template + void AMReX_openPMDWriter::DumpParticles(PC const& pc, + const std::string& name, + //int iteration, + const amrex::Vector& write_real_comp, + const amrex::Vector& write_int_comp, + const amrex::Vector& real_comp_names, + const amrex::Vector& int_comp_names, + FUNC_pc&& pc_func, + FUNC_pit&& posId_func + ) const + { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_Series != nullptr, "openPMD series must be initialized"); + + // TODO have to Count before this function due to const restriction + // Count is not a const function + //pc.CountParticles(); + + openPMD::Iteration currIteration = GetIteration( m_CurrentStep ); + openPMD::ParticleSpecies currSpecies = currIteration.particles[name]; + + // + // define positions & offsets + // + /*AllocatePtlProperties(currSpecies, write_real_comp, real_comp_names, write_int_comp, int_comp_names, + pc.m_PtlCounter.GetTotalNumParticles()); */ + + if ( AllocatePtlProperties(currSpecies, write_real_comp, real_comp_names, write_int_comp, int_comp_names, + pc.m_PtlCounter.GetTotalNumParticles()) ) + { + SetupRealProperties(pc, currSpecies, write_real_comp, real_comp_names, write_int_comp, int_comp_names, + pc.m_PtlCounter.GetTotalNumParticles()); + SetupConstant(currSpecies, pc.m_PtlCounter.GetTotalNumParticles()); + } + + pc_func(pc, currSpecies, pc.m_PtlCounter.GetTotalNumParticles()); + // open files from all processors, in case some will not contribute below + m_Series->flush(); + + bool emptyPC = true; + for (auto currentLevel = 0; currentLevel <= pc.finestLevel(); currentLevel++) + { + auto offset = static_cast( pc.m_PtlCounter.m_ParticleOffsetAtRank[currentLevel] ); + using ParIter = typename PC::ParConstIterType; + + //for (WarpXParIter pti(*pc, currentLevel); pti.isValid(); ++pti) { + for (ParIter pti(pc, currentLevel); pti.isValid(); ++pti) + { + auto const numParticleOnTile = pti.numParticles(); + auto const numParticleOnTile64 = static_cast( numParticleOnTile ); + + // Do not call storeChunk() with zero-sized particle tiles: + // https://github.com/openPMD/openPMD-api/issues/1147 + // https://github.com/ECP-WarpX/WarpX/pull/1898#discussion_r745008290 + if (numParticleOnTile == 0) { continue; } + + emptyPC = false; + + //SavePosId(pti, currSpecies, offset); + posId_func(pti, currSpecies, offset); + // save "extra" particle properties in AoS and SoA + SaveRealProperty(pti, + currSpecies, + offset + GetGrandOffset(), + write_real_comp,real_comp_names, + write_int_comp, int_comp_names); + + offset += numParticleOnTile64; + } + } + + if (emptyPC && m_Series->backend() == "ADIOS2") { + // TODO adiosWorkaround with empty rank reading error() for __BTD__ + } + m_Series->flush(); + } + + + template + void AMReX_openPMDWriter::SetupRealProperties (PC const& /*pc*/, + openPMD::ParticleSpecies& currSpecies, + const amrex::Vector& write_real_comp, + const amrex::Vector& real_comp_names, + const amrex::Vector& write_int_comp, + const amrex::Vector& int_comp_names, + const unsigned long long /*np*/) const + { + + //AllocateRealProperties(currSpecies, write_real_comp, real_comp_names, write_int_comp, int_comp_names, np); + + // set attributes for some properties + std::set< std::string > addedRecords; // add meta-data per record only once + + auto const real_counter = std::min(write_real_comp.size(), real_comp_names.size()) - PC::NStructReal; + //for (auto idx=0; idx) + + target_sources( amrex_${D}d + PRIVATE + AMReX_PlotFileUtilOPENPMD.H + AMReX_PlotFileUtilOPENPMD.cpp + AMReX_PlotFileOPENPMD.cpp + AMReX_PlotFileOPENPMD_PTL.cpp + AMReX_ParticlesOPENPMD.H + AMReX_PlotFileUtilOPENPMD_PTLImpl.H + ) +endforeach() diff --git a/Src/Extern/openPMD-api/Make.package b/Src/Extern/openPMD-api/Make.package new file mode 100644 index 00000000000..f6ce7cae466 --- /dev/null +++ b/Src/Extern/openPMD-api/Make.package @@ -0,0 +1,6 @@ +CEXE_sources += AMReX_PlotFileOPENPMD.cpp AMReX_PlotFileOPENPMD_PTL.cpp AMReX_PlotFileUtilOPENPMD.cpp + +CEXE_headers += AMReX_ParticlesOPENPMD.H AMReX_PlotFileUtilOPENPMD.H AMReX_PlotFileUtilOPENPMD_PTLImpl.H + +VPATH_LOCATIONS += $(AMREX_HOME)/Src/Extern/openPMD-api +INCLUDE_LOCATIONS += $(AMREX_HOME)/Src/Extern/openPMD-api diff --git a/Src/Particle/AMReX_ParticleContainer.H b/Src/Particle/AMReX_ParticleContainer.H index 568fe8472e1..98102d63f17 100644 --- a/Src/Particle/AMReX_ParticleContainer.H +++ b/Src/Particle/AMReX_ParticleContainer.H @@ -1355,6 +1355,10 @@ public: #include "AMReX_ParticlesHDF5.H" #endif +#ifdef AMREX_USE_OPENPMD +#include "AMReX_ParticlesOPENPMD.H" +#endif + protected: template diff --git a/Tests/Amr/Advection_AmrLevel/CMakeLists.txt b/Tests/Amr/Advection_AmrLevel/CMakeLists.txt index 95d35436af3..443619455e1 100644 --- a/Tests/Amr/Advection_AmrLevel/CMakeLists.txt +++ b/Tests/Amr/Advection_AmrLevel/CMakeLists.txt @@ -1,3 +1,7 @@ +if (NOT AMReX_AMRLEVEL) + return() +endif () + foreach(D IN LISTS AMReX_SPACEDIM) if (D EQUAL 1) continue() diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt index d15a1e89697..f82f79df81c 100644 --- a/Tests/CMakeLists.txt +++ b/Tests/CMakeLists.txt @@ -135,6 +135,10 @@ else() list(APPEND AMREX_TESTS_SUBDIRS HDF5Benchmark) endif () + if (AMReX_OPENPMD) + list(APPEND AMREX_TESTS_SUBDIRS openPMDTests) + endif () + if (AMReX_FORTRAN_INTERFACES) list(APPEND AMREX_TESTS_SUBDIRS FortranInterface) endif () diff --git a/Tests/EB/CNS/CMakeLists.txt b/Tests/EB/CNS/CMakeLists.txt index 4f4a9d51f93..2a289dfa4ab 100644 --- a/Tests/EB/CNS/CMakeLists.txt +++ b/Tests/EB/CNS/CMakeLists.txt @@ -1,3 +1,7 @@ +if (NOT AMReX_AMRLEVEL) + return() +endif () + if (NOT (3 IN_LIST AMReX_SPACEDIM) OR NOT CMAKE_Fortran_COMPILER_LOADED) return() endif () diff --git a/Tests/EB_CNS/CMakeLists.txt b/Tests/EB_CNS/CMakeLists.txt index 5cd0e0b577d..932f1b93f4c 100644 --- a/Tests/EB_CNS/CMakeLists.txt +++ b/Tests/EB_CNS/CMakeLists.txt @@ -1,3 +1,7 @@ +if (NOT AMReX_AMRLEVEL) + return() +endif () + foreach(D IN LISTS AMReX_SPACEDIM) set(_sources main.cpp CNS_advance.cpp CNS_advance_box.cpp CNS_advance_box_eb.cpp CNS_bcfill.cpp CNSBld.cpp CNS.cpp CNS.H CNS_derive.cpp CNS_derive.H CNS_index_macros.H CNS_init_eb2.cpp CNS_io.cpp diff --git a/Tests/GPU/CNS/CMakeLists.txt b/Tests/GPU/CNS/CMakeLists.txt index 9e037ee86d7..0e1a472bae7 100644 --- a/Tests/GPU/CNS/CMakeLists.txt +++ b/Tests/GPU/CNS/CMakeLists.txt @@ -1,3 +1,7 @@ +if (NOT AMReX_AMRLEVEL) + return() +endif () + if (NOT 3 IN_LIST AMReX_SPACEDIM) return() endif () diff --git a/Tests/openPMDTests/fields/CMakeLists.txt b/Tests/openPMDTests/fields/CMakeLists.txt new file mode 100644 index 00000000000..8e4fd3e7c0d --- /dev/null +++ b/Tests/openPMDTests/fields/CMakeLists.txt @@ -0,0 +1,10 @@ +foreach(D IN LISTS AMReX_SPACEDIM) + set(_sources main.cpp) + set(_input_files inputs ) + + setup_test(${D} _sources _input_files) + + unset(_sources) + unset(_input_files) +endforeach() + diff --git a/Tests/openPMDTests/fields/GNUmakefile b/Tests/openPMDTests/fields/GNUmakefile new file mode 100644 index 00000000000..7e2b86dff4f --- /dev/null +++ b/Tests/openPMDTests/fields/GNUmakefile @@ -0,0 +1,25 @@ +DEBUG = FALSE + +USE_MPI = TRUE +USE_OMP = FALSE + +USE_OPENPMD = TRUE + +COMP = gnu + +DIM = 3 + +AMREX_HOME = ../../.. + +include $(AMREX_HOME)/Tools/GNUMake/Make.defs + +include ./Make.package + +Pdirs := Base + +Ppack += $(foreach dir, $(Pdirs), $(AMREX_HOME)/Src/$(dir)/Make.package) + +include $(Ppack) + +include $(AMREX_HOME)/Tools/GNUMake/Make.rules + diff --git a/Tests/openPMDTests/fields/Make.package b/Tests/openPMDTests/fields/Make.package new file mode 100644 index 00000000000..6b4b865e8fc --- /dev/null +++ b/Tests/openPMDTests/fields/Make.package @@ -0,0 +1 @@ +CEXE_sources += main.cpp diff --git a/Tests/openPMDTests/fields/inputs b/Tests/openPMDTests/fields/inputs new file mode 100644 index 00000000000..51fdba0631a --- /dev/null +++ b/Tests/openPMDTests/fields/inputs @@ -0,0 +1,38 @@ +# Domain size +ncells = 64 + +# Maximum allowable size of each subdomain in the problem domain; +# this is used to decompose the domain for parallel calculations. +max_grid_size = 8 + +# Number of levels +nlevs = 1 + +# Number of components in the multifabs +ncomp = 6 + +# Number of particles per cell +nppc = 2 + +# Number of plot files to write +nplotfile = 2 + +# Number of plot files to write +nparticlefile = 2 + +# Time to sleep before each write +sleeptime = 2 + +# Whether to check the correctness of Checkpoint / Restart +restart_check = 1 + +## Uncomment to read grids from file +#nlevs = 3 +#grids_from_file = 1 +#ref_ratio_file = 4 2 + +## Uncomment to enable compression +#hdf5compression=ZFP_ACCURACY@0.001 +#hdf5compression=SZ@sz.config + +directory=. diff --git a/Tests/openPMDTests/fields/main.cpp b/Tests/openPMDTests/fields/main.cpp new file mode 100644 index 00000000000..04c9a2f4476 --- /dev/null +++ b/Tests/openPMDTests/fields/main.cpp @@ -0,0 +1,193 @@ +#include +#include +#include + +#include +#include + +using namespace amrex; + +struct InputParams { + int ncells; + int max_grid_size; + int ncomp; + int nlevs; + + //int restart_check = 0; + int nplotfile = 1; + //int sleeptime = 0; +}; + + +void loadParameters(InputParams& input) +{ + ParmParse pp; + pp.get("ncells", input.ncells); + pp.get("max_grid_size", input.max_grid_size); + pp.get("ncomp", input.ncomp); + pp.get("nlevs", input.nlevs); + + pp.query("nplotfile", input.nplotfile); + //pp.query("sleeptime", input.sleeptime); + //pp.query("restart_check", input.restart_check); + //pp.query("grids_from_file", input.grids_from_file); +} + +void set_grids_nested (const InputParams& input, + Vector& domains, + Vector& grids, + Vector& ref_ratio); + +void test (); + +int main(int argc, char* argv[]) +{ + amrex::Initialize(argc,argv); + test(); + amrex::Finalize(); + std::cout<<"Finalized "< domains; + Vector ba; + set_grids_nested(inputs, domains, ba, m_Ref_ratio); + + RealBox real_box; + for (int n = 0; n < AMREX_SPACEDIM; n++) { + real_box.setLo(n, 0.0); + real_box.setHi(n, 1.0); + } + + // This sets the boundary conditions to be doubly or triply periodic + int is_per[AMREX_SPACEDIM] = {AMREX_D_DECL(1,1,1)}; + + // This defines a Geometry object for each level + m_Geom.resize(inputs.nlevs); + + m_Geom[0].define(domains[0], &real_box, CoordSys::cartesian, is_per); + for (int lev = 1; lev < inputs.nlevs; lev++) { + m_Geom[lev].define(domains[lev], &real_box, CoordSys::cartesian, is_per); + } + + Vector dmap(inputs.nlevs); + + m_mf.resize(inputs.nlevs); + for (int lev = 0; lev < inputs.nlevs; lev++) { + std::cout<(ba[lev], dmap[lev], inputs.ncomp, nghost); + m_mf[lev]->setVal(lev); + } + + for (int i = 0; i < inputs.ncomp; ++i) { + m_Varnames.push_back("component_" + std::to_string(i)); + } + + //Vector level_steps(inputs.nlevs, 0); + m_Level_steps.assign(inputs.nlevs, 0); + } + + Vector m_Level_steps; + + Vector m_Varnames; + Real m_Time = Real(0.0); + Vector m_Ref_ratio; + Vector m_Geom; + Vector > m_mf; +}; + + +void saveFile(char const* fname, const InputParams& inputs, const TestField& testField) +{ +#ifdef AMREX_USE_OPENPMD + openpmd_api::InitHandler(fname); + + for (int ts = 0; ts < inputs.nplotfile; ts++) + { + openpmd_api::SetStep(ts); + openpmd_api::WriteMultiLevel(//fname, + amrex::GetVecOfConstPtrs(testField.m_mf), testField.m_Varnames, + testField.m_Geom, testField.m_Time, /*testField.m_Level_steps,*/ testField.m_Ref_ratio); + openpmd_api::CloseStep(ts); + + //saveFile(fname, ts, inputs, testField); + amrex::Print()<<"Timestep: "<& domains, + Vector& grids, + Vector& ref_ratio) +{ + //int ncells, max_grid_size, nlevs; + + AMREX_ALWAYS_ASSERT(input.nlevs < 2); // relax this later + + IntVect domain_lo(AMREX_D_DECL(0, 0, 0)); + IntVect domain_hi(AMREX_D_DECL(input.ncells-1, input.ncells-1, input.ncells-1)); + + domains.resize(input.nlevs); + domains[0].setSmall(domain_lo); + domains[0].setBig(domain_hi); + + ref_ratio.resize(input.nlevs-1); + for (int lev = 1; lev < input.nlevs; lev++) { + ref_ratio[lev-1] = IntVect(AMREX_D_DECL(2, 2, 2)); + } + + grids.resize(input.nlevs); + grids[0].define(domains[0]); + + // Now we make the refined level be the center eighth of the domain + if (input.nlevs > 1) { + int n_fine = input.ncells * ref_ratio[0][0]; + IntVect refined_lo(AMREX_D_DECL(n_fine/4,n_fine/4,n_fine/4)); + IntVect refined_hi(AMREX_D_DECL(3*n_fine/4-1,3*n_fine/4-1,3*n_fine/4-1)); + + // Build a box for the level 1 domain + Box refined_patch(refined_lo, refined_hi); + grids[1].define(refined_patch); + } + + // break the BoxArrays at both levels into max_grid_size^3 boxes + for (int lev = 0; lev < input.nlevs; lev++) { + grids[lev].maxSize(input.max_grid_size); + } + + for (int lev = 1; lev < input.nlevs; lev++) { + domains[lev] = amrex::refine(domains[lev-1], ref_ratio[lev-1]); + } +} diff --git a/Tests/openPMDTests/ptls/CMakeLists.txt b/Tests/openPMDTests/ptls/CMakeLists.txt new file mode 100644 index 00000000000..697e875c5be --- /dev/null +++ b/Tests/openPMDTests/ptls/CMakeLists.txt @@ -0,0 +1,10 @@ +if (AMReX_PARTICLES) + foreach(D IN LISTS AMReX_SPACEDIM) + set(_sources main.cpp mypc.H trilinear_deposition_K.H) + set(_input_files inputs ) + setup_test(${D} _sources _input_files) + + unset(_sources) + unset(_input_files) + endforeach() +endif() diff --git a/Tests/openPMDTests/ptls/GNUmakefile b/Tests/openPMDTests/ptls/GNUmakefile new file mode 100644 index 00000000000..0eaf664cf9d --- /dev/null +++ b/Tests/openPMDTests/ptls/GNUmakefile @@ -0,0 +1,34 @@ +AMREX_HOME = ../../../ + +DEBUG = TRUE +DEBUG = FALSE + +DIM = 3 + +COMP = gcc + +TINY_PROFILE = TRUE +USE_PARTICLES = TRUE + +PRECISION = DOUBLE + +USE_MPI = TRUE +USE_OMP = FALSE + +USE_OPENPMD = TRUE + +################################################### + +EBASE = main + +include $(AMREX_HOME)/Tools/GNUMake/Make.defs + +include ./Make.package + +Pdirs := Base Boundary Particle AmrCore + +Ppack += $(foreach dir, $(Pdirs), $(AMREX_HOME)/Src/$(dir)/Make.package) + +include $(Ppack) + +include $(AMREX_HOME)/Tools/GNUMake/Make.rules diff --git a/Tests/openPMDTests/ptls/Make.package b/Tests/openPMDTests/ptls/Make.package new file mode 100644 index 00000000000..cf5bcd54ae0 --- /dev/null +++ b/Tests/openPMDTests/ptls/Make.package @@ -0,0 +1,4 @@ +CEXE_headers += mypc.H +CEXE_headers += trilinear_deposition_K.H +CEXE_sources += main.cpp + diff --git a/Tests/openPMDTests/ptls/inputs b/Tests/openPMDTests/ptls/inputs new file mode 100644 index 00000000000..15791c46d4b --- /dev/null +++ b/Tests/openPMDTests/ptls/inputs @@ -0,0 +1,26 @@ + +# Domain size + +nx = 32 # number of grid points along the x axis +ny = 32 # number of grid points along the y axis +nz = 32 # number of grid points along the z axis + +#nx = 64 # number of grid points along the x axis +#ny = 64 # number of grid points along the y axis +#nz = 64 # number of grid points along the z axis + +#nx = 128 # number of grid points along the x axis +#ny = 128 # number of grid points along the y axis +#nz = 128 # number of grid points along the z axis + +# Maximum allowable size of each subdomain in the problem domain; +# this is used to decompose the domain for parallel calculations. +max_grid_size = 16 + +# Number of particles per cell +nppc = 2 + +# Verbosity +verbose = true # set to true to get more verbosity + +nlevs=2 diff --git a/Tests/openPMDTests/ptls/main.cpp b/Tests/openPMDTests/ptls/main.cpp new file mode 100644 index 00000000000..74fb7450e7d --- /dev/null +++ b/Tests/openPMDTests/ptls/main.cpp @@ -0,0 +1,464 @@ +#include + +#include +#include +#include +#include "AMReX_Particles.H" +#include "AMReX_PlotFileUtil.H" +#include + +#include "mypc.H" +#include "trilinear_deposition_K.H" + +#include "warpxBTD.H" +#include "warpxWriter.H" + +using namespace amrex; +#define TEST_BTD 1 + +struct TestParams { + int nx; + int ny; + int nz; + int nlevs; + int max_grid_size; + int nppc; + bool verbose; + + int nplotfile=1; +}; + +void checkMFBox(const TestParams& parms, + Vector outputMF) +{ + for (int lev=0; lev < parms.nlevs; lev++) + { + auto const* curr_mf = outputMF[lev]; + int const ncomp = curr_mf->nComp(); + amrex::Print()<<" checking boxes, lev="< rr(parms.nlevs-1); + for (int lev = 1; lev < parms.nlevs; lev++) { + rr[lev-1] = IntVect(AMREX_D_DECL(2,2,2)); + } + + RealBox real_box; + for (int n = 0; n < BL_SPACEDIM; n++) { + real_box.setLo(n, 0.0); + real_box.setHi(n, 1.0); + } + + IntVect domain_lo(AMREX_D_DECL(0, 0, 0)); + IntVect domain_hi(AMREX_D_DECL(parms.nx - 1, parms.ny - 1, parms.nz-1)); + const Box base_domain(domain_lo, domain_hi); + + // This sets the boundary conditions to be doubly or triply periodic + int is_per[BL_SPACEDIM] = {AMREX_D_DECL(1,1,1)}; + + Vector geom(parms.nlevs); + geom[0].define(base_domain, &real_box, CoordSys::cartesian, is_per); + for (int lev = 1; lev < parms.nlevs; lev++) { + geom[lev].define(amrex::refine(geom[lev-1].Domain(), rr[lev-1]), + &real_box, CoordSys::cartesian, is_per); + } + + Vector ba(parms.nlevs); + Vector dm(parms.nlevs); + + Box domain = base_domain; + IntVect size = IntVect(AMREX_D_DECL(parms.nx, parms.ny, parms.nz)); + for (int lev = 0; lev < parms.nlevs; ++lev) + { + ba[lev].define(domain); + ba[lev].maxSize(parms.max_grid_size); + dm[lev].define(ba[lev]); + domain.grow(-size/4); // fine level cover the middle of the coarse domain + domain.refine(2); + } + + Vector density1(parms.nlevs); + Vector density2(parms.nlevs); + + for (int lev = 0; lev < parms.nlevs; lev++) { + // one field comp for density1 + density1[lev].define(ba[lev], dm[lev], 1, nghost); + density1[lev].setVal(0.0); + // and two field comp for density2 + density2[lev].define(ba[lev], dm[lev], 2, nghost); + density2[lev].setVal(2.0); + } + + MyParticleContainer myPC(geom, dm, ba, rr); + myPC.SetVerbose(false); + + bool serialize = true; + if (ParallelDescriptor::NProcs() > 1) { + serialize = false; + } + + amrex::Long num_particles = (amrex::Long)parms.nppc * parms.nx * parms.ny * parms.nz; + amrex::Print() << serialize<< " Total number of particles :" << num_particles << '\n'; + + int iseed = 451; + double mass = 10.0; + + MyParticleContainer::ParticleInitData pdata = {{mass}, {}, {}, {}}; + myPC.InitRandom(num_particles, iseed, pdata, serialize); + + // + // Here we provide an example of one way to call ParticleToMesh + // + amrex::ParticleToMesh(myPC, GetVecOfPtrs(density1), 0, parms.nlevs-1, + [=] AMREX_GPU_DEVICE (const MyParticleContainer::ParticleType& p, + amrex::Array4 const& rho, + amrex::GpuArray const& plo, + amrex::GpuArray const& dxi) noexcept + { + ParticleInterpolator::Linear interp(p, plo, dxi); + + interp.ParticleToMesh(p, rho, 0, 0, 1, + [=] AMREX_GPU_DEVICE (const MyParticleContainer::ParticleType& part, int comp) + { + return part.rdata(comp); // no weighting + }); + }); + + // + // Here we provide an example of another way to call ParticleToMesh + // + int start_part_comp = 0; + int start_mesh_comp = 0; + int num_comp = 1; + + amrex::ParticleToMesh(myPC,GetVecOfPtrs(density2),0,parms.nlevs-1, + TrilinearDeposition{start_part_comp,start_mesh_comp,num_comp}); + + // + // Now write the output from each into separate plotfiles for comparison + // + + Vector varnames1, varnames2; + // varname1 is for density1 + varnames1.push_back("density"); // has openPMD component scalar + + // varname2 is for density2, has two components, so assigned two names + varnames2.push_back("velosity"); // has openPMD component scalar + varnames2.push_back("Ex"); // has openPMD component E/x + + Vector particle_varnames; + particle_varnames.push_back("mass"); + + int output_levs = parms.nlevs; + + Vector outputMF1(output_levs); + Vector outputMF2(output_levs); + Vector outputRR(output_levs); + for (int lev = 0; lev < output_levs; ++lev) { + outputMF1[lev] = &density1[lev]; + outputMF2[lev] = &density2[lev]; + outputRR[lev] = IntVect(AMREX_D_DECL(2, 2, 2)); + } + + //checkMFBox(parms, outputMF1); + //checkMFBox(parms, outputMF2); + + // call count ptls to prepare ahead of time + myPC.CountParticles(); + + std::string fname; + openpmd_api::InitHandler(fname); + + // one specie per particle container + std::string specieName="ptlSpecie"; + + int nsteps = 3; + for (int ts = 0; ts < nsteps; ts++) + { + //Vector level_steps; + //level_steps.push_back(ts); + //level_steps.push_back(ts); + + openpmd_api::SetStep(ts); + openpmd_api::WriteParticles(myPC, specieName /*ts*/); // with default component names + + char name[512]; + std::snprintf(name, sizeof name, "plotfile_%05d", ts); + + if ( 1 == parms.nlevs ) + { + // example to store coarse level + openpmd_api::WriteSingleLevel(*(outputMF1[0]), varnames1, geom[0], 0.0); + openpmd_api::WriteSingleLevel(*(outputMF2[0]), varnames2, geom[0], 0.0); + } + else + { + // store multi mesh levels + openpmd_api::WriteMultiLevel(//parms.nlevs, + outputMF1, + varnames1, + geom, + 0.0, + //level_steps, + outputRR + ); + // store multi mesh levels + openpmd_api::WriteMultiLevel(//parms.nlevs, + outputMF2, + varnames2, + geom, + 0.0, + //level_steps, + outputRR + ); + } + openpmd_api::CloseStep(ts); + + amrex::Print()<<"Timestep: "< rr(parms.nlevs-1); + for (int lev = 1; lev < parms.nlevs; lev++) { + rr[lev-1] = IntVect(AMREX_D_DECL(2,2,2)); + } + + RealBox real_box; + for (int n = 0; n < BL_SPACEDIM; n++) { + real_box.setLo(n, 0.0); + real_box.setHi(n, 1.0); + } + + IntVect domain_lo(AMREX_D_DECL(0, 0, 0)); + IntVect domain_hi(AMREX_D_DECL(parms.nx - 1, parms.ny - 1, parms.nz-1)); + const Box base_domain(domain_lo, domain_hi); + + // This sets the boundary conditions to be doubly or triply periodic + int is_per[BL_SPACEDIM] = {AMREX_D_DECL(1,1,1)}; + + Vector geom(parms.nlevs); + geom[0].define(base_domain, &real_box, CoordSys::cartesian, is_per); + for (int lev = 1; lev < parms.nlevs; lev++) { + geom[lev].define(amrex::refine(geom[lev-1].Domain(), rr[lev-1]), + &real_box, CoordSys::cartesian, is_per); + } + + Vector ba(parms.nlevs); + Vector dm(parms.nlevs); + + Box domain = base_domain; + IntVect size = IntVect(AMREX_D_DECL(parms.nx, parms.ny, parms.nz)); + for (int lev = 0; lev < parms.nlevs; ++lev) + { + ba[lev].define(domain); + ba[lev].maxSize(parms.max_grid_size); + dm[lev].define(ba[lev]); + domain.grow(-size/4); // fine level cover the middle of the coarse domain + domain.refine(2); + } + + Vector density1(parms.nlevs); + + for (int lev = 0; lev < parms.nlevs; lev++) { + // one field comp for density1 + density1[lev].define(ba[lev], dm[lev], 1, nghost); + density1[lev].setVal(0.0); + } + + MyParticleContainer myPC(geom, dm, ba, rr); + myPC.SetVerbose(false); + + bool serialize = true; + if (ParallelDescriptor::NProcs() > 1) { + serialize = false; + } + + amrex::Long num_particles = (amrex::Long)parms.nppc * parms.nx * parms.ny * parms.nz; + amrex::Print() << serialize<< " Total number of particles :" << num_particles << '\n'; + + int iseed = 451; + double mass = 10.0; + + MyParticleContainer::ParticleInitData pdata = {{mass}, {}, {}, {}}; + myPC.InitRandom(num_particles, iseed, pdata, serialize); + + // + // Here we provide an example of one way to call ParticleToMesh + // + amrex::ParticleToMesh(myPC, GetVecOfPtrs(density1), 0, parms.nlevs-1, + [=] AMREX_GPU_DEVICE (const MyParticleContainer::ParticleType& p, + amrex::Array4 const& rho, + amrex::GpuArray const& plo, + amrex::GpuArray const& dxi) noexcept + { + ParticleInterpolator::Linear interp(p, plo, dxi); + + interp.ParticleToMesh(p, rho, 0, 0, 1, + [=] AMREX_GPU_DEVICE (const MyParticleContainer::ParticleType& part, int comp) + { + return part.rdata(comp); // no weighting + }); + }); + + // + // Now write the output from each into separate plotfiles for comparison + // + + Vector varnames1; + // varname1 is for density1 + varnames1.push_back("density"); // has openPMD component scalar + + Vector particle_varnames; + particle_varnames.push_back("mass"); + + int output_levs = parms.nlevs; + + Vector outputMF1(output_levs); + Vector outputRR(output_levs); + for (int lev = 0; lev < output_levs; ++lev) { + outputMF1[lev] = &density1[lev]; + outputRR[lev] = IntVect(AMREX_D_DECL(2, 2, 2)); + } + + //checkMFBox(parms, outputMF1); + //checkMFBox(parms, outputMF2); + + // call count ptls to prepare ahead of time + myPC.CountParticles(); + + std::string fname; + openpmd_api::InitHandler(fname); + + std::vector warpxPMLs(10); + auto* testWriter = new AMReX_warpxBTDWriter(warpxPMLs); // xxxxx is there a memory leak? + amrex::openpmd_api::UseCustomWriter(testWriter); + + // one specie per particle container + std::string specieName="ptlSpecie"; + + int nsteps = 4; + // + // mimicking BTD behavior. Total is 2 actual steps, each steps is written twice + // step writing order is 0 1 0 1, the last two writes are the final flushes + // To make things simple, at the end of the second write, we should see double the ptls. + // AsssignPtlOffsets(num_ptls) makes sure the second write starts off correctly + // + for (int its = 0; its < nsteps; its++) + { + int ts = its % 2; + openpmd_api::SetStep(ts); + + if ( (its - 2) == ts ) + { + //call AssignPtloffset() to assign the right starting offset of ptl batch + testWriter->AssignPtlOffset(num_particles); + testWriter->SetLastFlush(); + } +#if 0 + if (0) + { // test with default component names + openpmd_api::WriteParticles(myPC, specieName); + } + else +#endif + { // test with RZ style pos id + openpmd_api::WriteParticles(myPC, + specieName, + [=] (auto& /*pc*/, openPMD::ParticleSpecies& currSpecies, unsigned long long localTotal) + { + amrex::ParticleReal lcharge = 0.01; // warpx: pc->getCharge() + amrex::ParticleReal lmass = 0.5; // warpx: pc->getMass(); + + testWriter->SetConstantMassCharge(currSpecies, localTotal, lcharge, lmass); + }, + [=] (auto& pti, openPMD::ParticleSpecies& currSpecies, unsigned long long offset) + { + testWriter->SavePosId_RZ(pti, currSpecies, offset); // also supports RZ + }); + } + + { + // store multi mesh levels + openpmd_api::WriteMultiLevel(//parms.nlevs, + outputMF1, + varnames1, + geom, + 0.0, + outputRR + ); + } + + openpmd_api::CloseStep(ts); + + amrex::Print()<<"Timestep: "< + +using MyParticleContainer = amrex::ParticleContainer<20,3,7,9>; +// change to two parameters so i can check with plotfile output +//typedef amrex::ParticleContainer<3,4> MyParticleContainer; +//typedef amrex::ParticleContainer<1> MyParticleContainer; diff --git a/Tests/openPMDTests/ptls/pseudo_warpx.H b/Tests/openPMDTests/ptls/pseudo_warpx.H new file mode 100644 index 00000000000..ba8ee79d9c8 --- /dev/null +++ b/Tests/openPMDTests/ptls/pseudo_warpx.H @@ -0,0 +1,54 @@ +#ifndef JOKE_WARPX +#define JOKE_WARPX + +namespace WarpX +{ + static short electromagnetic_solver_id = 2; + static bool do_dive_cleaning = true; + static bool use_filter = true; + static amrex::IntVect filter_npass_each_dir = {11,22,33}; + static int n_rz_azimuthal_modes = 999; + static int nox = 77; + static int noy = 88; + static int noz = 99; + static short particle_pusher_algo = 0; + static short field_gathering_algo = 1; + static short current_deposition_algo = 0; + + struct CurrentDepositionAlgo { + enum { + Esirkepov = 0, + Direct = 1, + Vay = 2 + }; + }; + + struct GatheringAlgo { + enum { + EnergyConserving = 0, + MomentumConserving + }; + }; + + struct ParticlePusherAlgo { + enum { + Boris = 0, + Vay = 1, + HigueraCary = 2 + }; + }; + + struct ElectromagneticSolverAlgo { + enum { + None = 0, + Yee = 1, + CKC = 2, + PSATD = 3, + Other = 4 + }; + }; + +} + + +#endif diff --git a/Tests/openPMDTests/ptls/trilinear_deposition_K.H b/Tests/openPMDTests/ptls/trilinear_deposition_K.H new file mode 100644 index 00000000000..fba78c9fc3e --- /dev/null +++ b/Tests/openPMDTests/ptls/trilinear_deposition_K.H @@ -0,0 +1,29 @@ + +#include +#include +#include +#include +#include +#include + +struct TrilinearDeposition +{ + int start_part_comp; + int start_mesh_comp; + int num_comp; + + AMREX_GPU_DEVICE AMREX_FORCE_INLINE + void operator() (const MyParticleContainer::ParticleType& p, + amrex::Array4 const& rho, + amrex::GpuArray const& plo, + amrex::GpuArray const& dxi) const noexcept + { + amrex::ParticleInterpolator::Linear interp(p, plo, dxi); + + interp.ParticleToMesh(p, rho, start_part_comp, start_mesh_comp, num_comp, + [=] AMREX_GPU_DEVICE (const MyParticleContainer::ParticleType& part, int comp) + { + return part.rdata(comp); // no weighting + }); + } +}; diff --git a/Tests/openPMDTests/ptls/warpxBTD.H b/Tests/openPMDTests/ptls/warpxBTD.H new file mode 100644 index 00000000000..8d7c47a5133 --- /dev/null +++ b/Tests/openPMDTests/ptls/warpxBTD.H @@ -0,0 +1,259 @@ +#include + +#include +#include +#include +#include "AMReX_Particles.H" +#include "AMReX_PlotFileUtil.H" +#include + +#include "warpxWriter.H" + + +//////////////////////////// +//class AMReX_warpxBTDWriter final : public amrex::openpmd_api::AMReX_openPMDWriter +class AMReX_warpxBTDWriter final : public AMReX_warpxWriter +{ + public: + AMReX_warpxBTDWriter(std::vector fieldPMLdirections, + openPMD::IterationEncoding ie = openPMD::IterationEncoding::groupBased, + std::string filetype = "default", + std::string openpmdOptions = "{}"); + + // BTD does not support streams + [[nodiscard]] openPMD::Iteration GetIteration (int const iteration) const override + { + return m_Series->iterations[iteration]; + } + + [[nodiscard]] unsigned long long GetGrandOffset() const override + { + return m_AssignedPtlOffset; + } + + /* + */ + bool AllocatePtlProperties(openPMD::ParticleSpecies& currSpecies, + const amrex::Vector& write_real_comp, + const amrex::Vector& real_comp_names, + const amrex::Vector& write_int_comp, + const amrex::Vector& int_comp_names, + unsigned long long np) const override; + + + void SetupConstant(openPMD::ParticleSpecies& currSpecies, + const unsigned long long& np) const override; + + void CloseStep(int ts) override; + + void Init(openPMD::Access access) override; + + // + // in warpx, btd ignores "geom" and takes in a single geometry "full_BTD_snapshot" + // and loops over all levels in the vector "geom" + // might as well supply equivalent vector of full_BTD_snapshot instead of both "geom" and "full_BTD_snapshot" + // + // TODO -- need to find out whether BTD is only for levels if all levels are using full_BTD_snapshot in warpX + // + void WriteMesh(const std::vector& varnames, + const amrex::Vector& mf, + const amrex::Vector& geom, + //const int iteration, + double time ) const override; + // BTD specific: + void SetLastFlush() {m_LastFlush = true;} + [[nodiscard]] bool IsLastFlush() const {return m_LastFlush;} + void AssignPtlOffset(unsigned long long m) { m_AssignedPtlOffset = m; } + + template + void SavePosId_RZ(PIt& pti, + openPMD::ParticleSpecies& currSpecies, + unsigned long long offset) const; + + void SetConstantMassCharge(openPMD::ParticleSpecies& currSpecies, + const unsigned long long& np, + amrex::ParticleReal charge, + amrex::ParticleReal mass) const override; + + private: + bool m_LastFlush = false; + unsigned long long m_AssignedPtlOffset = 0; + bool m_FirstSightOfStep = false; +}; // class warpxBTDwriter + + + + +AMReX_warpxBTDWriter::AMReX_warpxBTDWriter( + std::vector fieldPMLdirections, + openPMD::IterationEncoding ie, + std::string filetype, + std::string openpmdOptions) +// :AMReX_openPMDWriter("{}", ie, filetype, openpmdOptions) + :AMReX_warpxWriter(std::move(fieldPMLdirections), ie, std::move(filetype), std::move(openpmdOptions)) +{ + m_openPMDDatasetOptions="{ \"resizable\": true }"; + + if (ie == openPMD::IterationEncoding::variableBased) + { + std::string warnMsg = " Unable to support BTD with streaming. Using GroupBased "; + m_openPMDEncoding = openPMD::IterationEncoding::groupBased; + } +} + +void AMReX_warpxBTDWriter::CloseStep(int ts) +{ + if (!m_LastFlush) { + return; + } + + AMReX_openPMDWriter::CloseStep(ts); + m_LastFlush = false; // this step is over +} + +void AMReX_warpxBTDWriter::Init(openPMD::Access access) +{ + if (m_Series != nullptr) { + m_FirstSightOfStep = ! m_Series->iterations.contains( m_CurrentStep ); + return; + } + m_FirstSightOfStep = true; + AMReX_openPMDWriter::Init(access); +} + +bool AMReX_warpxBTDWriter::AllocatePtlProperties(openPMD::ParticleSpecies& currSpecies, + const amrex::Vector& write_real_comp, + const amrex::Vector& real_comp_names, + const amrex::Vector& write_int_comp, + const amrex::Vector& int_comp_names, + const unsigned long long np) const +{ + //amrex::Print() <<"allocate np="< 0) + { + //AMReX_openPMDWriter::SetupPos(currSpecies, np + m_AssignedPtlOffset); + return AMReX_openPMDWriter::AllocatePtlProperties(currSpecies, write_real_comp, real_comp_names, + write_int_comp, int_comp_names, np + GetGrandOffset()); + } + // np == 0 + if ( m_LastFlush && ( m_AssignedPtlOffset == 0 ) ) + { + // properties were never allocated + //AMReX_openPMDWriter::SetupPos(currSpecies, np); + auto r = AMReX_openPMDWriter::AllocatePtlProperties(currSpecies, write_real_comp, real_comp_names, write_int_comp, int_comp_names, np); + amrex::ignore_unused(r); + } + return false; +} + +/* +*/ +void AMReX_warpxBTDWriter::SetupConstant(openPMD::ParticleSpecies& currSpecies, + const unsigned long long& np) const +{ + if (m_LastFlush) + { + AMReX_openPMDWriter::SetupConstant(currSpecies, np + GetGrandOffset()); + } +} + +void AMReX_warpxBTDWriter::SetConstantMassCharge(openPMD::ParticleSpecies& currSpecies, + const unsigned long long& np, + amrex::ParticleReal const charge, + amrex::ParticleReal const mass) const +{ + if (m_LastFlush) + { + BL_PROFILE("SetChargeMass(BTD)"); + AMReX_warpxWriter::SetConstantMassCharge(currSpecies, np + GetGrandOffset(), + charge, mass); + } +} + + + +void AMReX_warpxBTDWriter::WriteMesh(const std::vector& varnames, + const amrex::Vector& mf, + const amrex::Vector& geom, + //const int iteration, + const double time ) const +{ + BL_PROFILE("WriteMesh(BTD)"); + if (m_FirstSightOfStep) { + AMReX_openPMDWriter::WriteMesh(varnames, mf, geom, /*iteration,*/ time); + return; + } + + // field and meshcomp were already setup for this iteration + // only need to + openPMD::Iteration series_iteration = GetIteration(m_CurrentStep); + series_iteration.open(); + + auto meshes = series_iteration.meshes; + for (int lev=0; lev < geom.size(); lev++) + { + amrex::Geometry full_geom = geom[lev]; + //CompSetup(lev, meshes, full_geom, varnames, mf[lev]); + CompStorage(lev, meshes, full_geom, varnames, mf[lev]); +#ifdef AMREX_USE_GPU + amrex::Gpu::streamSynchronize(); +#endif + m_Series->flush(); + } // for lev loop +} + + + +template +void +AMReX_warpxBTDWriter::SavePosId_RZ(PIt& pti, + openPMD::ParticleSpecies& currSpecies, + unsigned long long offset) const +{ + BL_PROFILE("amrex::openpmd_api::SavePosId( RZ ..)"); + +#if defined(WARPX_DIM_RZ) + const auto& aos = pti.GetArrayOfStructs(); // size = numParticlesOnTile + { + // Save positions + auto const positionComponents = /*helper::*/getParticlePositionComponentLabels(); + { + std::shared_ptr z( + new amrex::ParticleReal[numParticleOnTile], + [](amrex::ParticleReal const *p) { delete[] p; } + ); + for (auto i = 0; i < numParticleOnTile; i++) + z.get()[i] = aos[i].pos(1); // {0: "r", 1: "z"} + + std::string const positionComponent = "z"; + currSpecies["position"]["z"].storeChunk(z, {offset}, {numParticleOnTile64}); + } + + // reconstruct x and y from polar coordinates r, theta + auto const& soa = pti.GetStructOfArrays(); + amrex::ParticleReal const* theta = soa.GetRealData(PIdx::theta).dataPtr(); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(theta != nullptr, "openPMD: invalid theta pointer."); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(int(soa.GetRealData(PIdx::theta).size()) == numParticleOnTile, + "openPMD: theta and tile size do not match"); + { + std::shared_ptr< amrex::ParticleReal > x( + new amrex::ParticleReal[numParticleOnTile], + [](amrex::ParticleReal const *p){ delete[] p; } + ); + std::shared_ptr< amrex::ParticleReal > y( + new amrex::ParticleReal[numParticleOnTile], + [](amrex::ParticleReal const *p){ delete[] p; } + ); + for (auto i=0; i + +#include +#include +#include +#include "AMReX_Particles.H" +#include "AMReX_PlotFileUtil.H" +#include + +#include "pseudo_warpx.H" +#include "AMReX_PlotFileUtilOPENPMD.H" + +class AMReX_warpxWriter: public amrex::openpmd_api::AMReX_openPMDWriter +{ +public: + AMReX_warpxWriter(std::vector fieldPMLdirections, + openPMD::IterationEncoding ie = openPMD::IterationEncoding::variableBased, + std::string filetype = "default", + std::string openpmdOptions = "{}") + :AMReX_openPMDWriter("{}", ie, std::move(filetype), std::move(openpmdOptions)), + m_fieldPMLdirections(std::move(fieldPMLdirections)) + {} + + // + // fields related funcitions + // + void SetupMeshComp (openPMD::Mesh& mesh, + const amrex::Geometry& full_geom, + amrex::MultiFab const& mf, + const amrex::openpmd_api::AMReX_VarNameParser& varName) const override; + + void SetupFields (openPMD::Container< openPMD::Mesh >& meshes, + amrex::Geometry& full_geom) const override; + + [[nodiscard]] inline std::vector< std::string > getParticlePositionComponentLabels() const override; + + virtual void SetConstantMassCharge(openPMD::ParticleSpecies& currSpecies, + const unsigned long long& np, + amrex::ParticleReal charge, + amrex::ParticleReal mass) const; + + void SetParticleSpecieAttributes(openPMD::ParticleSpecies& currSpecies) const override; + +private: + std::vector m_fieldPMLdirections; + + [[nodiscard]] static std::vector< std::string > getFieldAxisLabels ( const amrex::openpmd_api::AMReX_VarNameParser& varName ); +}; + + +void AMReX_warpxWriter::SetupMeshComp (openPMD::Mesh& mesh, + const amrex::Geometry& full_geom, + amrex::MultiFab const& mf, + const amrex::openpmd_api::AMReX_VarNameParser& varName) const +{ + BL_PROFILE("SetupMeshComp(WarpX)"); + + auto mesh_comp = mesh[varName.m_CompName]; + amrex::Box const & global_box = full_geom.Domain(); + auto global_size = amrex::openpmd_api::helper::getReversedVec(global_box.size() ); + + // - Grid spacing + std::vector const grid_spacing = amrex::openpmd_api::helper::getReversedVec(full_geom.CellSize()); + mesh.setGridSpacing(grid_spacing); + + // - Global offset + std::vector const global_offset = amrex::openpmd_api::helper::getReversedVec(full_geom.ProbLo()); + mesh.setGridGlobalOffset(global_offset); + + // - AxisLabels + std::vector axis_labels = getFieldAxisLabels(varName); + mesh.setAxisLabels(axis_labels); + +#if defined(WARPX_DIM_RZ) + auto & warpx = WarpX::GetInstance(); + if (curr.m_ThetaMode) { + global_size.emplace(global_size.begin(), warpx.ncomps); + } +#endif + + // Prepare the type of dataset that will be written + openPMD::Datatype const datatype = openPMD::determineDatatype(); + auto const dataset = openPMD::Dataset(datatype, global_size); + mesh.setDataOrder(openPMD::Mesh::DataOrder::C); + + if (varName.m_ThetaMode) { + mesh.setGeometry("thetaMode"); + mesh.setGeometryParameters("m=" + std::to_string(WarpX::n_rz_azimuthal_modes) + ";imag=+"); + } + + mesh.setAttribute("fieldSmoothing", "none"); + mesh_comp.resetDataset(dataset); + + //amrex::openpmd_api::helper::setOpenPMDUnit( mesh, varName.m_FieldName ); // removed this function. define in app if needed + + auto relative_cell_pos = amrex::openpmd_api::helper::getRelativeCellPosition(mf); // AMReX Fortran index order + std::reverse( relative_cell_pos.begin(), relative_cell_pos.end() ); // now in C order + mesh_comp.setPosition( relative_cell_pos ); +} + +void AMReX_warpxWriter::SetupFields (openPMD::Container< openPMD::Mesh >& meshes, + amrex::Geometry& full_geom) const +{ + // meta data for ED-PIC extension + auto const period = full_geom.periodicity(); + + std::vector fieldBoundary(6, "reflecting"); + std::vector particleBoundary(6, "absorbing"); + fieldBoundary.resize(AMREX_SPACEDIM * 2); + particleBoundary.resize(AMREX_SPACEDIM * 2); + +#if AMREX_SPACEDIM != 3 + fieldBoundary.resize(4); + particleBoundary.resize(4); +#endif + if (m_fieldPMLdirections.size() >= fieldBoundary.size()) + { + for (int i = 0; i < int(fieldBoundary.size() / 2); ++i) { + if (m_fieldPMLdirections.at(i)) { + fieldBoundary.at(i) = "open"; + } + } + } + + for (int i = 0; i < int(fieldBoundary.size() / 2); ++i) { + if (period.isPeriodic(i)) { + fieldBoundary.at(2 * i) = "periodic"; + fieldBoundary.at(2 * i + 1) = "periodic"; + particleBoundary.at(2 * i) = "periodic"; + particleBoundary.at(2 * i + 1) = "periodic"; + } + } + + meshes.setAttribute("fieldSolver", []() { + switch (WarpX::electromagnetic_solver_id) { + case WarpX::ElectromagneticSolverAlgo::Yee : + { return "Yee"; } + case WarpX::ElectromagneticSolverAlgo::CKC : + { return "CK"; } + case WarpX::ElectromagneticSolverAlgo::PSATD : + { return "PSATD"; } + default: + { return "other"; } + } + }()) + ; + + meshes.setAttribute("fieldBoundary", fieldBoundary); + meshes.setAttribute("particleBoundary", particleBoundary); + + meshes.setAttribute("currentSmoothing", []() { + if (WarpX::use_filter) { return "Binomial"; } + else { return "none"; } + }()); + + if (WarpX::use_filter) { + meshes.setAttribute("currentSmoothingParameters", []() { + std::stringstream ss; + ss << "period=1;compensator=false"; +#if (AMREX_SPACEDIM >= 2) + ss << ";numPasses_x=" << WarpX::filter_npass_each_dir[0]; +#endif +#if (AMREX_SPACEDIM == 3) + ss << ";numPasses_y=" << WarpX::filter_npass_each_dir[1]; + ss << ";numPasses_z=" << WarpX::filter_npass_each_dir[2]; +#endif +#if defined(WARPX_DIM_3D) + ss << ";numPasses_y=" << WarpX::filter_npass_each_dir[1]; + ss << ";numPasses_z=" << WarpX::filter_npass_each_dir[2]; +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + ss << ";numPasses_z=" << WarpX::filter_npass_each_dir[1]; +#elif defined(WARPX_DIM_1D_Z) + ss << ";numPasses_z=" << WarpX::filter_npass_each_dir[0]; +#endif + std::string currentSmoothingParameters = ss.str(); + return currentSmoothingParameters; + }()); + } + + meshes.setAttribute("chargeCorrection", []() { + if (WarpX::do_dive_cleaning) { return "hyperbolic"; } // TODO or "spectral" or something? double-check + else { return "none"; } + }()); + + if (WarpX::do_dive_cleaning) { + meshes.setAttribute("chargeCorrectionParameters", "period=1"); + } +} + + + +std::vector< std::string > AMReX_warpxWriter::getFieldAxisLabels ( const amrex::openpmd_api::AMReX_VarNameParser& /*varName*/ ) + { + using vs = std::vector< std::string >; + // Fortran order of the index labels for the AMReX FArrayBox +#if defined(WARPX_DIM_1D_Z) + vs const axisLabels{"z"}; // z varies fastest in memory +#elif defined(WARPX_DIM_XZ) + vs const axisLabels{"x", "z"}; // x varies fastest in memory +#elif defined(WARPX_DIM_RZ) + // when we write individual modes of a field (default) + vs const circAxisLabels{"r", "z"}; // r varies fastest in memory + // if we just write reconstructed 2D fields at theta=0 + vs const cartAxisLabels{"x", "z"}; // x varies fastest in memory + vs const axisLabels = varName.m_ThetaMode ? circAxisLabels : cartAxisLabels; +#elif defined(WARPX_DIM_3D) + vs const axisLabels{"x", "y", "z"}; // x varies fastest in memory +#else + //error Unknown WarpX dimensionality. + vs const axisLabels{"w", "k", "i"}; // default as this is just a test +#endif + // revert to C order (fastest varying index last) + return {axisLabels.rbegin(), axisLabels.rend()}; + } + + // + // ptl related warpx functions + // +inline std::vector< std::string > +AMReX_warpxWriter::getParticlePositionComponentLabels() const + { + using vs = std::vector< std::string >; +#if defined(WARPX_DIM_1D_Z) + vs positionComponents{"z"}; +#elif defined(WARPX_DIM_XZ) + vs positionComponents{"x", "z"}; +#elif defined(WARPX_DIM_RZ) + // note: although we internally store particle positions + // for AMReX in r,z and a theta attribute, we + // actually need them for algorithms (e.g. push) + // and I/O in Cartesian. + // Other attributes like momentum are consequently + // stored in x,y,z internally. + vs positionComponents{"x", "y", "z"}; +#elif (WARPX_DIM_3D) + vs positionComponents{"x", "y", "z"}; +#else + //# error Unknown WarpX dimensionality. + // default: + vs positionComponents{"kk", "jj", "ii"}; +#endif + return positionComponents; + } + + // TODO should it be here? +void AMReX_warpxWriter::SetConstantMassCharge(openPMD::ParticleSpecies& currSpecies, + const unsigned long long& np, + amrex::ParticleReal const charge, + amrex::ParticleReal const mass) const + { + auto realType = openPMD::Dataset(openPMD::determineDatatype(), {np}); + auto const * const scalar = openPMD::RecordComponent::SCALAR; + + currSpecies["charge"][scalar].resetDataset( realType ); + currSpecies["mass"][scalar].resetDataset( realType ); + + currSpecies["charge"][scalar].makeConstant( charge ); + currSpecies["mass"][scalar].makeConstant( mass ); + + currSpecies["charge"].setUnitDimension( amrex::openpmd_api::helper::getUnitDimension("charge") ); + currSpecies["mass"].setUnitDimension( amrex::openpmd_api::helper::getUnitDimension("mass") ); + + currSpecies["charge"].setAttribute( "macroWeighted", 0U ); + currSpecies["charge"].setAttribute( "weightingPower", 1.0 ); + currSpecies["mass"].setAttribute( "macroWeighted", 0U ); + currSpecies["mass"].setAttribute( "weightingPower", 1.0 ); + } + + + + // should be called after + // - setupPos so records "position" and "id" are defined + // - setConstants so record "positonOffset" is defined +void AMReX_warpxWriter::SetParticleSpecieAttributes(openPMD::ParticleSpecies& currSpecies) const + { + + // ED-PIC extension + currSpecies["id"].setAttribute( "macroWeighted", 0U ); + currSpecies["id"].setAttribute( "weightingPower", 0.0 ); + currSpecies["position"].setAttribute( "macroWeighted", 0U ); + currSpecies["position"].setAttribute( "weightingPower", 0.0 ); + + currSpecies["positionOffset"].setAttribute( "macroWeighted", 0U ); + currSpecies["positionOffset"].setAttribute( "weightingPower", 0.0 ); + + // meta data for ED-PIC extension + currSpecies.setAttribute( "particleShape", double( WarpX::noz ) ); + + currSpecies.setAttribute( "particleShapes", [](){ + return std::vector< double >{ + double(WarpX::nox), +#if AMREX_SPACEDIM==3 + double(WarpX::noy), +#endif + double(WarpX::noz) + }; + }() ); + + currSpecies.setAttribute( "particlePush", [](){ + switch( WarpX::particle_pusher_algo ) { + case WarpX::ParticlePusherAlgo::Boris : { return "Boris"; } + case WarpX::ParticlePusherAlgo::Vay : { return "Vay"; } + case WarpX::ParticlePusherAlgo::HigueraCary : { return "HigueraCary"; } + default: { return "other"; } + } + }() ); + currSpecies.setAttribute( "particleInterpolation", [](){ + switch( WarpX::field_gathering_algo ) { + case WarpX::GatheringAlgo::EnergyConserving : { return "energyConserving"; } + case WarpX::GatheringAlgo::MomentumConserving : { return "momentumConserving"; } + default: { return "other"; } + } + }() ); + currSpecies.setAttribute( "particleSmoothing", "none" ); + currSpecies.setAttribute( "currentDeposition", [](){ + switch( WarpX::current_deposition_algo ) { + case WarpX::CurrentDepositionAlgo::Esirkepov : { return "Esirkepov"; } + default: { return "directMorseNielson"; } + } + }() ); + } + + +#endif diff --git a/Tools/CMake/AMReXConfig.cmake.in b/Tools/CMake/AMReXConfig.cmake.in index 6ee60df4990..e2792c77d48 100644 --- a/Tools/CMake/AMReXConfig.cmake.in +++ b/Tools/CMake/AMReXConfig.cmake.in @@ -85,6 +85,7 @@ set(AMReX_PETSC_FOUND @AMReX_PETSC@) set(AMReX_SUNDIALS_FOUND @AMReX_SUNDIALS@) set(AMReX_HDF5_FOUND @AMReX_HDF5@) set(AMReX_HDF5_ZFP_FOUND @AMReX_HDF5_ZFP@) +set(AMReX_OPENPMD_FOUND @AMReX_OPENPMD@) # Compilation options set(AMReX_FPE_FOUND @AMReX_FPE@) @@ -137,6 +138,7 @@ set(AMReX_HYPRE @AMReX_HYPRE@) set(AMReX_PETSC @AMReX_PETSC@) set(AMReX_HDF5 @AMReX_HDF5@) set(AMReX_HDF5_ZFP @AMReX_HDF5_ZFP@) +set(AMReX_OPENPMD @AMReX_OPENPMD@) # Compilation options set(AMReX_FPE @AMReX_FPE@) @@ -215,6 +217,10 @@ if (@AMReX_HYPRE@) find_dependency(HYPRE 2.20.0 REQUIRED) endif () +if (@AMReX_OPENPMD@) + find_dependency(openpmd REQUIRED) +endif () + if (@AMReX_PETSC@) find_dependency(PETSc 2.13 REQUIRED) endif () diff --git a/Tools/CMake/AMReXOptions.cmake b/Tools/CMake/AMReXOptions.cmake index d19477e7698..fc221f7be64 100644 --- a/Tools/CMake/AMReXOptions.cmake +++ b/Tools/CMake/AMReXOptions.cmake @@ -352,6 +352,10 @@ cmake_dependent_option(AMReX_HDF5_ZFP "Enable ZFP compression in HDF5-based IO" "AMReX_HDF5" OFF ) print_option(AMReX_HDF5_ZFP) +# openPMD-api +option(AMReX_OPENPMD "Enable I/O through openPMD-api" OFF) +print_option(AMReX_OPENPMD) + # SUNDIALS option( AMReX_SUNDIALS "Enable SUNDIALS interfaces" OFF ) print_option( AMReX_SUNDIALS ) diff --git a/Tools/CMake/AMReXSetDefines.cmake b/Tools/CMake/AMReXSetDefines.cmake index fb428070809..3d6c348443d 100644 --- a/Tools/CMake/AMReXSetDefines.cmake +++ b/Tools/CMake/AMReXSetDefines.cmake @@ -170,6 +170,11 @@ add_amrex_define(AMREX_USE_HDF5 NO_LEGACY IF AMReX_HDF5) add_amrex_define(AMREX_USE_HDF5_ASYNC NO_LEGACY IF AMReX_HDF5_ASYNC) add_amrex_define(AMREX_USE_HDF5_ZFP NO_LEGACY IF AMReX_HDF5_ZFP) +# +# openPMD-api +# +add_amrex_define(AMREX_USE_OPENPMD NO_LEGACY IF AMReX_OPENPMD) + # # SUNDIALS diff --git a/Tools/CMake/AMReXThirdPartyLibraries.cmake b/Tools/CMake/AMReXThirdPartyLibraries.cmake index f8f49e9c478..a73ace6a5f8 100644 --- a/Tools/CMake/AMReXThirdPartyLibraries.cmake +++ b/Tools/CMake/AMReXThirdPartyLibraries.cmake @@ -1,4 +1,17 @@ # +# openPMD-api related targets +# +if (AMReX_OPENPMD) + if (AMReX_MPI) + find_package(openPMD 0.14.2 CONFIG REQUIRED COMPONENTS MPI) + else () + find_package(openPMD 0.14.2 CONFIG REQUIRED COMPONENTS NOMPI) + endif() + foreach(D IN LISTS AMReX_SPACEDIM) + target_link_libraries(amrex_${D}d PUBLIC openPMD::openPMD) + endforeach() +endif() +# # HDF5 -- here it would be best to create an imported target # if (AMReX_HDF5) diff --git a/Tools/CMake/AMReX_Config_ND.H.in b/Tools/CMake/AMReX_Config_ND.H.in index 5e499df1af1..fe2ae734c88 100644 --- a/Tools/CMake/AMReX_Config_ND.H.in +++ b/Tools/CMake/AMReX_Config_ND.H.in @@ -59,6 +59,7 @@ #cmakedefine AMREX_USE_HDF5 #cmakedefine AMREX_USE_HDF5_ASYNC #cmakedefine AMREX_USE_HDF5_ZFP +#cmakedefine AMREX_USE_OPENPMD #cmakedefine AMREX_USE_HYPRE #cmakedefine AMREX_USE_PETSC #cmakedefine AMREX_USE_SUNDIALS diff --git a/Tools/GNUMake/Make.defs b/Tools/GNUMake/Make.defs index 3394ce0d9ad..7b6578b3ec2 100644 --- a/Tools/GNUMake/Make.defs +++ b/Tools/GNUMake/Make.defs @@ -299,6 +299,12 @@ else USE_HDF5 := FALSE endif +ifdef USE_OPENPMD + USE_OPENPMD := $(strip $(USE_OPENPMD)) +else + USE_OPENPMD := FALSE +endif + ifdef EBASE EBASE := $(strip $(EBASE)) else @@ -1079,6 +1085,11 @@ ifeq ($(USE_HDF5),TRUE) include $(AMREX_HOME)/Tools/GNUMake/packages/Make.hdf5 endif +ifeq ($(USE_OPENPMD),TRUE) + $(info Loading $(AMREX_HOME)/Tools/GNUMake/packages/Make.openpmd...) + include $(AMREX_HOME)/Tools/GNUMake/packages/Make.openpmd +endif + ifeq ($(USE_BITTREE),TRUE) $(info Loading $(AMREX_HOME)/Tools/GNUMake/packages/Make.bittree...) include $(AMREX_HOME)/Tools/GNUMake/packages/Make.bittree diff --git a/Tools/GNUMake/packages/Make.openpmd b/Tools/GNUMake/packages/Make.openpmd new file mode 100644 index 00000000000..f512e697130 --- /dev/null +++ b/Tools/GNUMake/packages/Make.openpmd @@ -0,0 +1,20 @@ +CPPFLAGS += -DAMREX_USE_OPENPMD +include $(AMREX_HOME)/Src/Extern/openPMD-api/Make.package + +ifndef AMREX_OPENPMD_HOME +ifdef OPENPMD_DIR + AMREX_OPENPMD_HOME = $(OPENPMD_DIR) +endif +ifdef OPENPMD_HOME + AMREX_OPENPMD_HOME = $(OPENPMD_HOME) +endif +endif + +LIBRARIES += -lopenPMD + +ifdef AMREX_OPENPMD_HOME + OPENPMD_ABSPATH = $(abspath $(AMREX_OPENPMD_HOME)) + SYSTEM_INCLUDE_LOCATIONS += $(OPENPMD_ABSPATH)/include + LIBRARY_LOCATIONS += $(OPENPMD_ABSPATH)/lib + LDFLAGS += -Xlinker -rpath -Xlinker $(OPENPMD_ABSPATH)/lib +endif