diff options
| -rw-r--r-- | CONTRIBUTING.md | 7 | ||||
| -rw-r--r-- | Docs/source/running_cpp/parameters.rst | 6 | ||||
| -rw-r--r-- | Examples/Physics_applications/laser_acceleration/inputs_3d | 2 | ||||
| -rw-r--r-- | Source/Diagnostics/WarpXOpenPMD.H | 35 | ||||
| -rw-r--r-- | Source/Diagnostics/WarpXOpenPMD.cpp | 367 | ||||
| -rw-r--r-- | Source/Parallelization/InterpolateCurrentFineToCoarse.H | 4 | ||||
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 1 | ||||
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.cpp | 1 | ||||
| -rw-r--r-- | Source/WarpX.H | 9 | ||||
| -rw-r--r-- | Source/WarpX.cpp | 21 | ||||
| -rwxr-xr-x | run_test.sh | 16 |
11 files changed, 365 insertions, 104 deletions
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index b9de662f2..90d7bf95f 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -95,7 +95,7 @@ Once your new feature is ready, you can check that you did not break anything. WarpX has automated tests running for each Pull Request. For easier debugging, it can be convenient to run the tests on your local machine with ``` -./run_tests.sh +./run_test.sh ``` from WarpX root folder. The tests can be influenced by environment variables: - `export WARPX_TEST_DIM=3`, `export WARPX_TEST_DIM=2` or `export WARPX_TEST_DIM=RZ` @@ -104,6 +104,11 @@ in order to select only the tests that correspond to this dimensionality run the tests on CPU or GPU respectively. - `export WARPX_TEST_COMMIT=...` in order to test a specific commit. +The command above (without command line arguments) runs all the tests defined in [Regression/WarpX-tests.ini](./Regression/WarpX-tests.ini). In order to run single tests, pass the test names as command line arguments: +``` +./run_test.sh test1 test2 +``` + ### Submit a Pull Request A Pull Request (PR) is the way to efficiently visualize the changes you made diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 0ee7740fa..3cd7af96c 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -386,6 +386,10 @@ Particle initialization If `1` is given, both charge deposition and current deposition will not be done, thus that species does not contribute to the fields. +* ``<species_name>.do_not_push`` (`0` or `1` optional; default `0`) + If `1` is given, this species will not be pushed + by any pusher during the simulation. + * ``<species>.plot_species`` (`0` or `1` optional; default `1`) Whether to plot particle quantities for this species. @@ -1085,7 +1089,7 @@ Diagnostics and output time interval is expressed in the laboratory frame). * ``slice.particle_slice_width_lab`` (`float`, in meters) - Only used when ``warpx.do_boosted_frame_diagnostic`` is ``1`` and + Only used when ``warpx.do_back_transformed_diagnostics`` is ``1`` and ``slice.num_slice_snapshots_lab`` is non-zero. Particles are copied from the full back-transformed diagnostic to the reduced slice diagnostic if there are within the user-defined width from diff --git a/Examples/Physics_applications/laser_acceleration/inputs_3d b/Examples/Physics_applications/laser_acceleration/inputs_3d index 017b82962..3a0368e72 100644 --- a/Examples/Physics_applications/laser_acceleration/inputs_3d +++ b/Examples/Physics_applications/laser_acceleration/inputs_3d @@ -54,7 +54,7 @@ electrons.momentum_distribution_type = "constant" electrons.do_continuous_injection = 1 ################################# -############ PLASMA ############# +############ LASER ############# ################################# lasers.nlasers = 1 lasers.names = laser1 diff --git a/Source/Diagnostics/WarpXOpenPMD.H b/Source/Diagnostics/WarpXOpenPMD.H index c79a12066..e4b48d605 100644 --- a/Source/Diagnostics/WarpXOpenPMD.H +++ b/Source/Diagnostics/WarpXOpenPMD.H @@ -1,10 +1,20 @@ #ifndef WARPX_OPEN_PMD_H_ #define WARPX_OPEN_PMD_H_ -#include <MultiParticleContainer.H> // has AMReX Vector etc used below +#include "WarpXParticleContainer.H" +#include "MultiParticleContainer.H" // PIdx + +#include <AMReX_ParallelDescriptor.H> +#include <AMReX_REAL.H> +#include <AMReX_Utility.H> #include <openPMD/openPMD.hpp> +#include <memory> +#include <string> +#include <vector> + + // // helper class // @@ -61,10 +71,14 @@ private: class WarpXOpenPMDPlot { public: - // not using const string, to allow std::move to be effective - WarpXOpenPMDPlot(bool, std::string& filetype); + /** Initialize openPMD I/O routines + * + * @param oneFilePerTS write one file per timestep + * @param filetype file backend, e.g. "bp" or "h5" + * @param fieldPMLdirections PML field solver, @see WarpX::getPMLdirections() + */ + WarpXOpenPMDPlot(bool oneFilePerTS, std::string filetype, std::vector<bool> fieldPMLdirections); - //WarpXOpenPMDPlot(const std::string& dir, const std::string& fileType); ~WarpXOpenPMDPlot(); void SetStep(int ts); @@ -82,12 +96,14 @@ private: void Init(//const std::string& filename, openPMD::AccessType accessType); - /** This function sets up the entries for storing the particle positions in an openPMD species + /** This function sets up the entries for storing the particle positions, global IDs, and constant records (charge, mass) * - * @param[in] currSpecies The openPMD species - * @param[in] np Number of particles + * @param[in] pc WarpX particle container + * @param[in] currSpecies Corresponding openPMD species + * @param[in] np Number of particles */ - void SetupPos(openPMD::ParticleSpecies& currSpecies, + void SetupPos(const std::unique_ptr<WarpXParticleContainer>& pc, + openPMD::ParticleSpecies& currSpecies, const unsigned long long& np) const ; /** This function sets up the entries for particle properties @@ -149,6 +165,9 @@ private: bool m_OneFilePerTS = true; //! write in openPMD fileBased manner for individual time steps std::string m_OpenPMDFileType = "bp"; //! MPI-parallel openPMD backend: bp or h5 int m_CurrentStep = -1; + + // meta data + std::vector< bool > m_fieldPMLdirections; //! @see WarpX::getPMLdirections() }; diff --git a/Source/Diagnostics/WarpXOpenPMD.cpp b/Source/Diagnostics/WarpXOpenPMD.cpp index 8e0bca549..ed2bf8020 100644 --- a/Source/Diagnostics/WarpXOpenPMD.cpp +++ b/Source/Diagnostics/WarpXOpenPMD.cpp @@ -1,12 +1,28 @@ #include "WarpXOpenPMD.H" +#include "WarpXAlgorithmSelection.H" #include "FieldIO.H" // for getReversedVec +#include <algorithm> +#include <cstdint> +#include <map> +#include <set> +#include <string> +#include <sstream> #include <tuple> #include <utility> namespace detail { + + /** Convert AMReX AoS .id and .cpu to a globally unique particle ID + */ + union GlobalID { + struct { int id; int cpu; }; //! MPI-rank local ID (and rank/cpu) + uint64_t global_id; //! global ID that is unique in the whole simulation + }; + static_assert(sizeof(int) * 2u <= sizeof(uint64_t), "int size might cause collisions in global IDs"); + /** Unclutter a real_names to openPMD record * * @param fullName name as in real_names variable @@ -25,14 +41,55 @@ namespace detail } return make_pair(record_name, component_name); } + + /** Get the openPMD physical dimensionality of a record + * + * @param record_name name of the openPMD record + * @return map with base quantities and power scaling + */ + 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 {}; + } } WarpXOpenPMDPlot::WarpXOpenPMDPlot(bool oneFilePerTS, - std::string& openPMDFileType) + std::string openPMDFileType, std::vector<bool> fieldPMLdirections) :m_Series(nullptr), m_OneFilePerTS(oneFilePerTS), - m_OpenPMDFileType(std::move(openPMDFileType)) - //m_OpenPMDFileType(openPMDFileType) + m_OpenPMDFileType(std::move(openPMDFileType)), + m_fieldPMLdirections(std::move(fieldPMLdirections)) {} WarpXOpenPMDPlot::~WarpXOpenPMDPlot() @@ -104,13 +161,13 @@ WarpXOpenPMDPlot::Init(openPMD::AccessType accessType) if( WarpX::authors.size() > 0u ) m_Series->setAuthor( WarpX::authors ); // more natural naming for PIC - m_Series->setMeshesPath("fields"); - // TODO conform to ED-PIC extension of openPMD - // uint32_t const openPMD_ED_PIC = 1u; - // m_Series->setOpenPMDextension(openPMD_ED_PIC); + 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("WarpX"); - m_Series->setSoftwareVersion(WARPX_GIT_VERSION); + m_Series->setSoftware( "WarpX" ); + m_Series->setSoftwareVersion( WarpX::Version() ) ; } @@ -124,6 +181,7 @@ WarpXOpenPMDPlot::WriteOpenPMDParticles(const std::unique_ptr<MultiParticleConta auto& pc = mpc->GetUniqueContainer(i); if (pc->plot_species) { + // names of amrex::Real and int particle attributes in SoA data amrex::Vector<std::string> real_names; amrex::Vector<std::string> int_names; amrex::Vector<int> int_flags; @@ -197,21 +255,55 @@ WarpXOpenPMDPlot::SavePlotFile (const std::unique_ptr<WarpXParticleContainer>& p openPMD::Iteration currIteration = m_Series->iterations[iteration]; openPMD::ParticleSpecies currSpecies = currIteration.particles[name]; + // meta data for ED-PIC extension + currSpecies.setAttribute( "particleShape", double( WarpX::noz ) ); + // TODO allow this per direction in the openPMD standard, ED-PIC extension? + 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 ParticlePusherAlgo::Boris : return "Boris"; + case ParticlePusherAlgo::Vay : return "Vay"; + case ParticlePusherAlgo::HigueraCary : return "HigueraCary"; + default: return "other"; + } + }() ); + currSpecies.setAttribute( "particleInterpolation", [](){ + switch( WarpX::field_gathering_algo ) { + case GatheringAlgo::EnergyConserving : return "energyConserving"; + case GatheringAlgo::MomentumConserving : return "momentumConserving"; + default: return "other"; + } + }() ); + currSpecies.setAttribute( "particleSmoothing", "none" ); + currSpecies.setAttribute( "currentDeposition", [](){ + switch( WarpX::current_deposition_algo ) { + case CurrentDepositionAlgo::Esirkepov : return "Esirkepov"; + default: return "directMorseNielson"; + } + }() ); + // // define positions & offsets // - SetupPos(currSpecies, counter.GetTotalNumParticles()); + SetupPos(pc, currSpecies, counter.GetTotalNumParticles()); SetupRealProperties(currSpecies, write_real_comp, real_comp_names, counter.GetTotalNumParticles()); // forces the files created by all processors! this is the key to resolve RZ storage issue!! m_Series->flush(); for (auto currentLevel = 0; currentLevel <= pc->finestLevel(); currentLevel++) { - //long numParticles = counter.m_ParticleSizeAtRank[currentLevel] - unsigned long long const numParticles = counter.m_ParticleSizeAtRank[currentLevel]; - unsigned long long offset = counter.m_ParticleOffsetAtRank[currentLevel]; + uint64_t const numParticles = static_cast<uint64_t>( counter.m_ParticleSizeAtRank[currentLevel] ); + uint64_t offset = static_cast<uint64_t>( counter.m_ParticleOffsetAtRank[currentLevel] ); - if (0 == numParticles) + if (0u == numParticles) return; // pc->NumIntComp() & NumRealComp() are protected, @@ -219,22 +311,33 @@ WarpXOpenPMDPlot::SavePlotFile (const std::unique_ptr<WarpXParticleContainer>& p // soa num real attributes = PIdx::nattribs, and num int in soa is 0 for (WarpXParIter pti(*pc, currentLevel); pti.isValid(); ++pti) { - auto numParticleOnTile = pti.numParticles(); + auto const numParticleOnTile = pti.numParticles(); + uint64_t const numParticleOnTile64 = static_cast<uint64_t>( numParticleOnTile ); - // get position from aos + // get position and particle ID from aos + // note: this implementation iterates the AoS 4x but saves memory... iterating once could be beneficial, too const auto& aos = pti.GetArrayOfStructs(); // size = numParticlesOnTile - { // :: Save Postions, 1D-3D:: - // this code is tested with laser 3d, not tested with 2D examples... + { + // Save positions std::vector<std::string> axisNames={"x", "y", "z"}; - for (auto currDim = 0; currDim < AMREX_SPACEDIM; currDim++) { - std::vector<amrex::ParticleReal> curr(numParticleOnTile, 0); + std::vector<amrex::ParticleReal> curr(numParticleOnTile, 0.); for (auto i=0; i<numParticleOnTile; i++) { curr[i] = aos[i].m_rdata.pos[currDim]; } - currSpecies["position"][axisNames[currDim]].storeChunk(curr, {offset}, {static_cast<unsigned long long>(numParticleOnTile)}); + currSpecies["position"][axisNames[currDim]].storeChunk(curr, {offset}, {numParticleOnTile64}); m_Series->flush(); } + + // save particle ID after converting it to a globally unique ID + std::vector<uint64_t> ids(numParticleOnTile, 0.); + for (auto i=0; i<numParticleOnTile; i++) { + detail::GlobalID const nextID = { aos[i].m_idata.id, aos[i].m_idata.cpu }; + ids[i] = nextID.global_id; + } + auto const scalar = openPMD::RecordComponent::SCALAR; + currSpecies["id"][scalar].storeChunk(ids, {offset}, {numParticleOnTile64}); + m_Series->flush(); } // save properties SaveRealProperty(pti, @@ -242,7 +345,7 @@ WarpXOpenPMDPlot::SavePlotFile (const std::unique_ptr<WarpXParticleContainer>& p offset, write_real_comp, real_comp_names); - offset += numParticleOnTile; + offset += numParticleOnTile64; } } } @@ -274,39 +377,40 @@ WarpXOpenPMDPlot::SetupRealProperties(openPMD::ParticleSpecies& currSpecies, void WarpXOpenPMDPlot::SaveRealProperty(WarpXParIter& pti, openPMD::ParticleSpecies& currSpecies, - unsigned long long offset, - const amrex::Vector<int>& write_real_comp, - const amrex::Vector<std::string>& real_comp_names) const + unsigned long long const offset, + amrex::Vector<int> const& write_real_comp, + amrex::Vector<std::string> const& real_comp_names) const { int numOutputReal = 0; int const totalRealAttrs = m_NumAoSRealAttributes + m_NumSoARealAttributes; - for (int i = 0; i < totalRealAttrs; ++i) - if (write_real_comp[i]) + for( int i = 0; i < totalRealAttrs; ++i ) + if( write_real_comp[i] ) ++numOutputReal; - auto numParticleOnTile = pti.numParticles(); - const auto& aos = pti.GetArrayOfStructs(); // size = numParticlesOnTile - const auto& soa = pti.GetStructOfArrays(); + auto const numParticleOnTile = pti.numParticles(); + auto const& aos = pti.GetArrayOfStructs(); // size = numParticlesOnTile + auto const& soa = pti.GetStructOfArrays(); // properties are saved separately { - for (auto idx=0; idx<m_NumAoSRealAttributes; idx++) { - if (write_real_comp[idx]) { + for( auto idx=0; idx<m_NumAoSRealAttributes; idx++ ) { + if( write_real_comp[idx] ) { // handle scalar and non-scalar records by name std::string record_name, component_name; std::tie(record_name, component_name) = detail::name2openPMD(real_comp_names[idx]); + auto& currRecord = currSpecies[record_name]; + auto& currRecordComp = currRecord[component_name]; - auto& currVar = currSpecies[record_name][component_name]; - typename amrex::ParticleReal *d = - static_cast<typename amrex::ParticleReal*> (malloc(sizeof(typename amrex::ParticleReal) * numParticleOnTile)); + amrex::ParticleReal *d = + static_cast<amrex::ParticleReal*>( malloc(sizeof(amrex::ParticleReal) * numParticleOnTile) ); - for (auto kk=0; kk<numParticleOnTile; kk++) + for( auto kk=0; kk<numParticleOnTile; kk++ ) d[kk] = aos[kk].m_rdata.arr[AMREX_SPACEDIM+idx]; - std::shared_ptr <typename amrex::ParticleReal> data(d, free); - currVar.storeChunk(data, + std::shared_ptr<amrex::ParticleReal> data(d, free); + currRecordComp.storeChunk(data, {offset}, {static_cast<unsigned long long>(numParticleOnTile)}); m_Series->flush(); } @@ -314,16 +418,30 @@ WarpXOpenPMDPlot::SaveRealProperty(WarpXParIter& pti, } { + std::set< std::string > addedRecords; // add meta-data per record only once for (auto idx=0; idx<m_NumSoARealAttributes; idx++) { auto ii = m_NumAoSRealAttributes + idx; if (write_real_comp[ii]) { // handle scalar and non-scalar records by name std::string record_name, component_name; std::tie(record_name, component_name) = detail::name2openPMD(real_comp_names[ii]); - - auto& currVar = currSpecies[record_name][component_name]; - currVar.storeChunk(openPMD::shareRaw(soa.GetRealData(idx)), - {offset}, {static_cast<unsigned long long>(numParticleOnTile)}); + auto& currRecord = currSpecies[record_name]; + auto& currRecordComp = currRecord[component_name]; + + // meta data for ED-PIC extension + bool newRecord = false; + std::tie(std::ignore, newRecord) = addedRecords.insert(record_name); + if( newRecord ) { + currRecord.setAttribute( "macroWeighted", 0u ); + if( record_name == "momentum" ) + currRecord.setAttribute( "weightingPower", 1.0 ); + else + currRecord.setAttribute( "weightingPower", 0.0 ); + currRecord.setUnitDimension( detail::getUnitDimension(record_name) ); + } + + currRecordComp.storeChunk(openPMD::shareRaw(soa.GetRealData(idx)), + {offset}, {static_cast<unsigned long long>(numParticleOnTile)}); } } m_Series->flush(); @@ -333,25 +451,44 @@ WarpXOpenPMDPlot::SaveRealProperty(WarpXParIter& pti, void -WarpXOpenPMDPlot::SetupPos(openPMD::ParticleSpecies& currSpecies, - const unsigned long long& np) const +WarpXOpenPMDPlot::SetupPos(const std::unique_ptr<WarpXParticleContainer>& pc, + openPMD::ParticleSpecies& currSpecies, + const unsigned long long& np) const { - auto particleLineup = openPMD::Dataset(openPMD::determineDatatype<amrex::ParticleReal>(), {np}); - - currSpecies["positionOffset"]["x"].resetDataset(particleLineup); - currSpecies["positionOffset"]["x"].makeConstant(0); - currSpecies["positionOffset"]["y"].resetDataset(particleLineup); - currSpecies["positionOffset"]["y"].makeConstant(0); - currSpecies["positionOffset"]["z"].resetDataset(particleLineup); - currSpecies["positionOffset"]["z"].makeConstant(0); - - //auto positions = openPMD::Dataset(openPMD::determineDatatype<amrex::ParticleReal>(), {np}); - currSpecies["position"]["x"].resetDataset(particleLineup); - currSpecies["position"]["y"].resetDataset(particleLineup); - currSpecies["position"]["z"].resetDataset(particleLineup); -} + auto const realType = openPMD::Dataset(openPMD::determineDatatype<amrex::ParticleReal>(), {np}); + auto const idType = openPMD::Dataset(openPMD::determineDatatype< uint64_t >(), {np}); + for( auto const& comp : {"x", "y", "z"} ) { + currSpecies["positionOffset"][comp].resetDataset( realType ); + currSpecies["positionOffset"][comp].makeConstant( 0. ); + currSpecies["position"][comp].resetDataset( realType ); + } + auto const scalar = openPMD::RecordComponent::SCALAR; + currSpecies["id"][scalar].resetDataset( idType ); + currSpecies["charge"][scalar].resetDataset( realType ); + currSpecies["charge"][scalar].makeConstant( pc->getCharge() ); + currSpecies["mass"][scalar].resetDataset( realType ); + currSpecies["mass"][scalar].makeConstant( pc->getMass() ); + + // meta data + currSpecies["position"].setUnitDimension( detail::getUnitDimension("position") ); + currSpecies["positionOffset"].setUnitDimension( detail::getUnitDimension("positionOffset") ); + currSpecies["charge"].setUnitDimension( detail::getUnitDimension("charge") ); + currSpecies["mass"].setUnitDimension( detail::getUnitDimension("mass") ); + + // meta data for ED-PIC extension + currSpecies["position"].setAttribute( "macroWeighted", 0u ); + currSpecies["position"].setAttribute( "weightingPower", 0.0 ); + currSpecies["positionOffset"].setAttribute( "macroWeighted", 0u ); + currSpecies["positionOffset"].setAttribute( "weightingPower", 0.0 ); + currSpecies["id"].setAttribute( "macroWeighted", 0u ); + currSpecies["id"].setAttribute( "weightingPower", 0.0 ); + currSpecies["charge"].setAttribute( "macroWeighted", 0u ); + currSpecies["charge"].setAttribute( "weightingPower", 1.0 ); + currSpecies["mass"].setAttribute( "macroWeighted", 0u ); + currSpecies["mass"].setAttribute( "weightingPower", 1.0 ); +} // @@ -365,48 +502,109 @@ WarpXOpenPMDPlot::WriteOpenPMDFields( //const std::string& filename, const int iteration, const double time ) const { - //This is AmrEx's tiny profiler. Possibly will apply it later + //This is AMReX's tiny profiler. Possibly will apply it later BL_PROFILE("WarpXOpenPMDPlot::WriteOpenPMDFields()"); if ( nullptr == m_Series) return; - const int ncomp = mf.nComp(); + int const ncomp = mf.nComp(); // Create a few vectors that store info on the global domain // Swap the indices for each of them, since AMReX data is Fortran order // and since the openPMD API assumes contiguous C order // - Size of the box, in integer number of cells - const amrex::Box& global_box = geom.Domain(); - auto global_size = getReversedVec(global_box.size()); + amrex::Box const & global_box = geom.Domain(); + auto const global_size = getReversedVec(global_box.size()); // - Grid spacing - std::vector<double> grid_spacing = getReversedVec(geom.CellSize()); + std::vector<double> const grid_spacing = getReversedVec(geom.CellSize()); // - Global offset - std::vector<double> global_offset = getReversedVec(geom.ProbLo()); + std::vector<double> const global_offset = getReversedVec(geom.ProbLo()); // - AxisLabels #if AMREX_SPACEDIM==3 - std::vector<std::string> axis_labels{"x", "y", "z"}; + std::vector<std::string> const axis_labels{"x", "y", "z"}; #else - std::vector<std::string> axis_labels{"x", "z"}; + std::vector<std::string> const axis_labels{"x", "z"}; #endif // Prepare the type of dataset that will be written - openPMD::Datatype datatype = openPMD::determineDatatype<amrex::Real>(); - auto dataset = openPMD::Dataset(datatype, global_size); + openPMD::Datatype const datatype = openPMD::determineDatatype<amrex::Real>(); + auto const dataset = openPMD::Dataset(datatype, global_size); + // meta data auto series_iteration = m_Series->iterations[iteration]; series_iteration.setTime( time ); + // meta data for ED-PIC extension + auto const period = geom.periodicity(); // TODO double-check: is this the proper global bound or of some level? + std::vector< std::string > fieldBoundary( 6, "reflecting" ); + std::vector< std::string > particleBoundary( 6, "absorbing" ); +#if AMREX_SPACEDIM!=3 + fieldBoundary.resize(4); + particleBoundary.resize(4); +#endif + + for( auto i = 0u; i < fieldBoundary.size() / 2u; ++i ) + if( m_fieldPMLdirections.at( i ) ) + fieldBoundary.at( i ) = "open"; + + for( auto i = 0u; i < fieldBoundary.size() / 2u; ++i ) + if( period.isPeriodic( i ) ) { + fieldBoundary.at(2u*i ) = "periodic"; + fieldBoundary.at(2u*i + 1u) = "periodic"; + particleBoundary.at(2u*i ) = "periodic"; + particleBoundary.at(2u*i + 1u) = "periodic"; + } + + auto meshes = series_iteration.meshes; + meshes.setAttribute( "fieldSolver", [](){ +#ifdef WARPX_USE_PSATD + return "PSATD"; // TODO double-check if WARPX_USE_PSATD_HYBRID is covered +#else + switch( WarpX::particle_pusher_algo ) { + case MaxwellSolverAlgo::Yee : return "Yee"; + case MaxwellSolverAlgo::CKC : return "CK"; + default: return "other"; + } +#endif + }() ); + 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"; + ss << ";numPasses_x=" << WarpX::filter_npass_each_dir[0]; +#if (AMREX_SPACEDIM == 3) + ss << ";numPasses_y=" << WarpX::filter_npass_each_dir[1]; + ss << ";numPasses_z=" << WarpX::filter_npass_each_dir[2]; +#else + ss << ";numPasses_z=" << WarpX::filter_npass_each_dir[1]; +#endif + std::string currentSmoothingParameters = ss.str(); + return std::move(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"); + // Loop through the different components, i.e. different fields stored in mf for (int icomp=0; icomp<ncomp; icomp++){ // Check if this field is a vector or a scalar, and extract the field name - const std::string& varname = varnames[icomp]; + std::string const & varname = varnames[icomp]; std::string field_name = varname; std::string comp_name = openPMD::MeshRecordComponent::SCALAR; - for (const char* vector_field: {"E", "B", "j"}){ - for (const char* comp: {"x", "y", "z"}){ - if (varname[0] == *vector_field && varname[1] == *comp ){ + for( char const* vector_field: {"E", "B", "j"} ) { + for( char const* comp: {"x", "y", "z"} ) { + if( varname[0] == *vector_field && varname[1] == *comp ) { field_name = varname[0] + varname.substr(2); // Strip component comp_name = varname[1]; } @@ -414,35 +612,36 @@ WarpXOpenPMDPlot::WriteOpenPMDFields( //const std::string& filename, } // Setup the mesh record accordingly - auto mesh = series_iteration.meshes[field_name]; - mesh.setDataOrder(openPMD::Mesh::DataOrder::F); // MultiFab: Fortran order + auto mesh = meshes[field_name]; + mesh.setDataOrder( openPMD::Mesh::DataOrder::F ); // MultiFab: Fortran order of indices and axes mesh.setAxisLabels( axis_labels ); mesh.setGridSpacing( grid_spacing ); mesh.setGridGlobalOffset( global_offset ); + mesh.setAttribute( "fieldSmoothing", "none" ); setOpenPMDUnit( mesh, field_name ); // Create a new mesh record component, and store the associated metadata auto mesh_comp = mesh[comp_name]; mesh_comp.resetDataset( dataset ); // Cell-centered data: position is at 0.5 of a cell size. - mesh_comp.setPosition(std::vector<double>{AMREX_D_DECL(0.5, 0.5, 0.5)}); + mesh_comp.setPosition( std::vector<double>{AMREX_D_DECL(0.5, 0.5, 0.5)} ); // Loop through the multifab, and store each box as a chunk, // in the openPMD file. - for (amrex::MFIter mfi(mf); mfi.isValid(); ++mfi) { + for( amrex::MFIter mfi(mf); mfi.isValid(); ++mfi ) { - const amrex::FArrayBox& fab = mf[mfi]; - const amrex::Box& local_box = fab.box(); + amrex::FArrayBox const& fab = mf[mfi]; + amrex::Box const& local_box = fab.box(); // Determine the offset and size of this chunk - amrex::IntVect box_offset = local_box.smallEnd() - global_box.smallEnd(); - auto chunk_offset = getReversedVec(box_offset); - auto chunk_size = getReversedVec(local_box.size()); + amrex::IntVect const box_offset = local_box.smallEnd() - global_box.smallEnd(); + auto const chunk_offset = getReversedVec( box_offset ); + auto const chunk_size = getReversedVec( local_box.size() ); // Write local data - const double* local_data = fab.dataPtr(icomp); - mesh_comp.storeChunk(openPMD::shareRaw(local_data), - chunk_offset, chunk_size); + double const * local_data = fab.dataPtr( icomp ); + mesh_comp.storeChunk( openPMD::shareRaw(local_data), + chunk_offset, chunk_size ); } } // Flush data to disk after looping over all components diff --git a/Source/Parallelization/InterpolateCurrentFineToCoarse.H b/Source/Parallelization/InterpolateCurrentFineToCoarse.H index cbbcdfab5..58451c6b7 100644 --- a/Source/Parallelization/InterpolateCurrentFineToCoarse.H +++ b/Source/Parallelization/InterpolateCurrentFineToCoarse.H @@ -63,9 +63,9 @@ public: // return zero for out-of-bounds accesses during interpolation // this is efficiently used as a method to add neutral elements beyond guards in the average below - auto const fine = [fine_unsafe] AMREX_GPU_DEVICE (int const j, int const k, int const l) noexcept + auto const fine = [fine_unsafe] AMREX_GPU_DEVICE (int const jj, int const kk, int const ll) noexcept { - return fine_unsafe.contains(j, k, l) ? fine_unsafe(j, k, l) : amrex::Real{0.}; + return fine_unsafe.contains(jj, kk, ll) ? fine_unsafe(jj, kk, ll) : amrex::Real{0.}; }; int const ii = i * m_refinement_ratio; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 5d56a5c42..21c8ddd31 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -40,6 +40,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_splitting", do_splitting); pp.query("split_type", split_type); pp.query("do_not_deposit", do_not_deposit); + pp.query("do_not_push", do_not_push); pp.query("do_continuous_injection", do_continuous_injection); pp.query("initialize_self_fields", initialize_self_fields); diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 03e83c5f3..85c4dc0d5 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -130,7 +130,6 @@ WarpXParticleContainer::ReadParameters () do_tiling = true; #endif pp.query("do_tiling", do_tiling); - pp.query("do_not_push", do_not_push); initialized = true; } diff --git a/Source/WarpX.H b/Source/WarpX.H index 4620ca6de..1549dded2 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -119,6 +119,9 @@ public: static long particle_pusher_algo; static int maxwell_fdtd_solver_id; + // div E cleaning + static int do_dive_cleaning; + // Interpolation order static long nox; static long noy; @@ -177,6 +180,9 @@ public: const amrex::MultiFab& getEfield_fp (int lev, int direction) {return *Efield_fp[lev][direction];} const amrex::MultiFab& getBfield_fp (int lev, int direction) {return *Bfield_fp[lev][direction];} + /** get low-high-low-high-... vector for each direction indicating if mother grid PMLs are enabled */ + std::vector<bool> getPMLdirections() const; + static amrex::MultiFab* getCosts (int lev) { if (m_instance) { return m_instance->costs[lev].get(); @@ -573,9 +579,6 @@ private: amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_buf; amrex::Vector<std::unique_ptr<amrex::MultiFab> > charge_buf; - // div E cleaning - int do_dive_cleaning = 0; - // PML int do_pml = 1; int pml_ncell = 10; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 49f8b5eaa..5f6f2b2e5 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -55,6 +55,7 @@ long WarpX::charge_deposition_algo; long WarpX::field_gathering_algo; long WarpX::particle_pusher_algo; int WarpX::maxwell_fdtd_solver_id; +int WarpX::do_dive_cleaning = 0; long WarpX::n_rz_azimuthal_modes = 1; long WarpX::ncomps = 1; @@ -152,7 +153,7 @@ WarpX::WarpX () ReadParameters(); #ifdef WARPX_USE_OPENPMD - m_OpenPMDPlotWriter = new WarpXOpenPMDPlot(openpmd_tspf, openpmd_backend); + m_OpenPMDPlotWriter = new WarpXOpenPMDPlot(openpmd_tspf, openpmd_backend, WarpX::getPMLdirections()); #endif // Geometry on all levels has been defined already. @@ -1217,6 +1218,24 @@ WarpX::GetPML (int lev) } } +std::vector< bool > +WarpX::getPMLdirections() const +{ + std::vector< bool > dirsWithPML( 6, false ); +#if AMREX_SPACEDIM!=3 + dirsWithPML.resize( 4 ); +#endif + if( do_pml ) + { + for( auto i = 0u; i < dirsWithPML.size() / 2u; ++i ) + { + dirsWithPML.at( 2u*i ) = bool(do_pml_Lo[i]); + dirsWithPML.at( 2u*i + 1u ) = bool(do_pml_Hi[i]); + } + } + return dirsWithPML; +} + void WarpX::BuildBufferMasks () { diff --git a/run_test.sh b/run_test.sh index b6d4a644b..67f4b48cf 100755 --- a/run_test.sh +++ b/run_test.sh @@ -12,6 +12,12 @@ # Use `export WARPX_TEST_ARCH=CPU` or `export WARPX_TEST_ARCH=GPU` in order # to run the tests on CPU or GPU respectively. +# Parse command line arguments: if test names are given as command line arguments, +# store them in variable tests_arg and define new command line argument to call +# regtest.py with option --tests (works also for single test) +tests_arg=$* +tests_run=${tests_arg:+--tests=${tests_arg}} + # Remove contents and link to a previous test directory (intentionally two arguments) rm -rf test_dir/* test_dir # Create a temporary test directory @@ -34,7 +40,7 @@ cd test_dir # Clone PICSAR and AMReX git clone --branch development https://github.com/AMReX-Codes/amrex.git # Use QED brach for QED tests -if [ ${HAS_QED} = "TRUE" ]; then +if [ "${HAS_QED}" = "TRUE" ]; then git clone --branch QED https://bitbucket.org/berkeleylab/picsar.git else git clone --branch master https://bitbucket.org/berkeleylab/picsar.git @@ -49,6 +55,12 @@ cd warpx/Regression python prepare_file_travis.py cp travis-tests.ini ../../rt-WarpX -# Run the tests +# Run tests cd ../../regression_testing/ +# run only tests specified in variable tests_arg (single test or multiple tests) +if [[ ! -z "${tests_arg}" ]]; then +python regtest.py ../rt-WarpX/travis-tests.ini --no_update all --source_git_hash=${WARPX_TEST_COMMIT} "${tests_run}" +# run all tests (variables tests_arg and tests_run are empty) +else python regtest.py ../rt-WarpX/travis-tests.ini --no_update all --source_git_hash=${WARPX_TEST_COMMIT} +fi |