From d7fd47ffd5fc160d8e975aabb3854d456f172cec Mon Sep 17 00:00:00 2001 From: Piet Jarmatz Date: Mon, 10 Aug 2026 19:50:37 +0200 Subject: [PATCH 01/16] SimpleMD host/device KOKKOS view separation compiles --- simplemd/LinkedCell.h | 45 +++++-- simplemd/MoleculeContainer.cpp | 123 ++++++++++++------ simplemd/MoleculeContainer.h | 82 ++++++++++-- .../ComputeMeanVelocityMapping.h | 1 + .../cell-mappings/ComputeTemperatureMapping.h | 1 + .../cell-mappings/DeleteMoleculesMapping.h | 1 + .../cell-mappings/LennardJonesForceMapping.h | 1 + .../LennardJonesPotentialEnergyMapping.h | 1 + .../ParallelBoundaryEmptyCellsMapping.h | 1 + ...iodicAndParallelBoundaryFillCellsMapping.h | 1 + .../PeriodicBoundaryEmptyCellsMapping.h | 1 + .../cell-mappings/ProfilePlotterMapping.h | 1 + simplemd/cell-mappings/RDFMapping.h | 1 + .../ResetPotentialEnergyMapping.h | 1 + .../cell-mappings/VaryCheckpointMapping.h | 1 + simplemd/molecule-mappings/Adios2Writer.h | 1 + .../ComputeMeanVelocityMapping.h | 1 + .../ConvertForcesFixedToFloatMapping.h | 1 + .../ConvertForcesFloatToFixedMapping.h | 1 + .../InitialPositionAndForceUpdate.h | 1 + .../SetMeanVelocityMapping.h | 1 + .../molecule-mappings/VTKMoleculeWriter.cpp | 2 +- .../molecule-mappings/VTKMoleculeWriter.h | 5 +- .../VelocityStoermerVerletMapping.h | 1 + .../WriteCheckPointMapping.h | 1 + .../LennardJonesForceMapping.h | 1 + test/integration/SimpleMDBench.h | 1 + 27 files changed, 218 insertions(+), 61 deletions(-) diff --git a/simplemd/LinkedCell.h b/simplemd/LinkedCell.h index 35d71fbbd..66b2733c8 100644 --- a/simplemd/LinkedCell.h +++ b/simplemd/LinkedCell.h @@ -17,8 +17,9 @@ class LinkedCell; class simplemd::LinkedCell { public: KOKKOS_FUNCTION LinkedCell() : _moleculeData(nullptr), _linkedCellNumMolecules(nullptr), _cellIndex(0), _isGhostCell(false) {} - KOKKOS_FUNCTION LinkedCell(const Kokkos::View* moleculeData, - const Kokkos::View* nMolecules, unsigned int cellIndex, bool isGhostCell) + + KOKKOS_FUNCTION LinkedCell(const void* moleculeData, + const void* nMolecules, unsigned int cellIndex, bool isGhostCell) : _moleculeData(moleculeData), _linkedCellNumMolecules(nMolecules), _cellIndex(cellIndex), _isGhostCell(isGhostCell) {} class Iterator { @@ -75,7 +76,14 @@ class simplemd::LinkedCell { *getMolecule(moleculeIdx) = *getMolecule(numMolecules() - 1); changeMoleculeCount(-1); } - KOKKOS_FUNCTION void clear() { (*_linkedCellNumMolecules)(_cellIndex) = 0; } + KOKKOS_FUNCTION void clear() { + KOKKOS_IF_ON_HOST(( + (*static_cast::host_mirror_type*>(_linkedCellNumMolecules))(_cellIndex) = 0; + )) + KOKKOS_IF_ON_DEVICE(( + (*static_cast*>(_linkedCellNumMolecules))(_cellIndex) = 0; + )) + } std::string to_string() const { std::stringstream to_ret; @@ -83,19 +91,40 @@ class simplemd::LinkedCell { return to_ret.str(); } - KOKKOS_INLINE_FUNCTION unsigned int numMolecules() const { return (*_linkedCellNumMolecules)(_cellIndex); } + KOKKOS_INLINE_FUNCTION unsigned int numMolecules() const { + KOKKOS_IF_ON_HOST(( + return (*static_cast::host_mirror_type*>(_linkedCellNumMolecules))(_cellIndex); + )) + KOKKOS_IF_ON_DEVICE(( + return (*static_cast*>(_linkedCellNumMolecules))(_cellIndex); + )) + } KOKKOS_INLINE_FUNCTION bool isGhostCell() const { return _isGhostCell; } KOKKOS_INLINE_FUNCTION size_t getIndex() const { return _cellIndex; } private: - KOKKOS_INLINE_FUNCTION void changeMoleculeCount(int by) { (*_linkedCellNumMolecules)(_cellIndex) += by; } + KOKKOS_INLINE_FUNCTION void changeMoleculeCount(int by) { + KOKKOS_IF_ON_HOST(( + (*static_cast::host_mirror_type*>(_linkedCellNumMolecules))(_cellIndex) += by; + )) + KOKKOS_IF_ON_DEVICE(( + (*static_cast*>(_linkedCellNumMolecules))(_cellIndex) += by; + )) + } - KOKKOS_INLINE_FUNCTION Molecule* getMolecule(unsigned int moleculeIndex) const { return &(*_moleculeData)(_cellIndex, moleculeIndex); } + KOKKOS_INLINE_FUNCTION Molecule* getMolecule(unsigned int moleculeIndex) const { + KOKKOS_IF_ON_HOST(( + return &(*static_cast::host_mirror_type*>(_moleculeData))(_cellIndex, moleculeIndex); + )) + KOKKOS_IF_ON_DEVICE(( + return &(*static_cast*>(_moleculeData))(_cellIndex, moleculeIndex); + )) + } - const Kokkos::View* _moleculeData; - const Kokkos::View* _linkedCellNumMolecules; + const void* _moleculeData; + const void* _linkedCellNumMolecules; const unsigned int _cellIndex; const bool _isGhostCell; }; diff --git a/simplemd/MoleculeContainer.cpp b/simplemd/MoleculeContainer.cpp index ffbffec5e..37244c903 100644 --- a/simplemd/MoleculeContainer.cpp +++ b/simplemd/MoleculeContainer.cpp @@ -11,10 +11,13 @@ simplemd::MoleculeContainer::MoleculeContainer(simplemd::services::ParallelTopol _cellCapacity(cellCapacity), _domainSize(parallelTopologyService.getGlobalDomainSize()), _domainOffset(parallelTopologyService.getGlobalDomainOffset()), _meshWidth(parallelTopologyService.getMeshWidth()), _globalIndexOfFirstCell(parallelTopologyService.getGlobalIndexOfFirstCell()), _localIndexOfFirstCell(parallelTopologyService.getLocalIndexOfFirstCell()), - _moleculeData("moleculeData", parallelTopologyService.getLocalNumberOfCellsLinear(true), cellCapacity), - _linkedCellNumMolecules("linkedCellNumMolecules", parallelTopologyService.getLocalNumberOfCellsLinear(true)), - _linkedCellIsGhostCell("linkedCellIsGhostCell", _linkedCellNumMolecules.size()), - _neighborOffsets("neighborOffsets", 26) { + _moleculeData_d("moleculeData", parallelTopologyService.getLocalNumberOfCellsLinear(true), cellCapacity), + _moleculeData_h(Kokkos::create_mirror_view(_moleculeData_d)), + _linkedCellNumMolecules_d("linkedCellNumMolecules", parallelTopologyService.getLocalNumberOfCellsLinear(true)), + _linkedCellNumMolecules_h(Kokkos::create_mirror_view(_linkedCellNumMolecules_d)), + _linkedCellIsGhostCell("linkedCellIsGhostCell", _linkedCellNumMolecules_d.size()), + _neighborOffsets("neighborOffsets", 26), + _memoryState(MemoryState::SYNCED) { auto host_mirror1 = Kokkos::create_mirror_view(_linkedCellIsGhostCell); for (unsigned int i = 0; i < _linkedCellIsGhostCell.size(); i++) @@ -34,54 +37,58 @@ simplemd::MoleculeContainer::MoleculeContainer(simplemd::services::ParallelTopol } void simplemd::MoleculeContainer::insert(unsigned int cellIdx, const simplemd::Molecule& molecule) { + synchronizeMemory(HostOperation{}); #if (MD_ERROR == MD_YES) - if (_linkedCellNumMolecules(cellIdx) + 1 > _cellCapacity) { + if (_linkedCellNumMolecules_h(cellIdx) + 1 > _cellCapacity) { Kokkos::printf("Cell capacity=%d would be exceeded by an operation! Exiting...", _cellCapacity); Kokkos::abort("simplemd::MoleculeContainer::insert\n"); } #endif - _moleculeData(cellIdx, _linkedCellNumMolecules(cellIdx)) = molecule; - _linkedCellNumMolecules(cellIdx) += 1; + _moleculeData_h(cellIdx, _linkedCellNumMolecules_h(cellIdx)) = molecule; + _linkedCellNumMolecules_h(cellIdx) += 1; } void simplemd::MoleculeContainer::insert(const simplemd::Molecule& molecule) { insert(positionToCellIndex(molecule.getConstPosition()), molecule); } void simplemd::MoleculeContainer::remove(unsigned int cellIdx, unsigned int moleculeIdx) { + synchronizeMemory(HostOperation{}); #if (MD_ERROR == MD_YES) - if (moleculeIdx >= _linkedCellNumMolecules(cellIdx)) { + if (moleculeIdx >= _linkedCellNumMolecules_h(cellIdx)) { Kokkos::printf("Deleting particle that does not exist! moleculeIdx: %d, cellIdx: %d, num molecules: %d\n", moleculeIdx, cellIdx, - _linkedCellNumMolecules(cellIdx)); + _linkedCellNumMolecules_h(cellIdx)); Kokkos::abort("ERROR simplemd::MoleculeContainer::remove\n"); } - if (cellIdx >= _linkedCellNumMolecules(cellIdx)) { + if (cellIdx >= _linkedCellNumMolecules_h(cellIdx)) { Kokkos::printf("Deleting particle from cell that does not exist! moleculeIdx: %d, cellIdx: %d, num molecules: %d\n", moleculeIdx, cellIdx, - _linkedCellNumMolecules(cellIdx)); + _linkedCellNumMolecules_h(cellIdx)); Kokkos::abort("ERROR simplemd::MoleculeContainer::remove\n"); } #endif - _moleculeData(cellIdx, moleculeIdx) = _moleculeData(cellIdx, _linkedCellNumMolecules(cellIdx) - 1); - _linkedCellNumMolecules(cellIdx) -= 1; + _moleculeData_h(cellIdx, moleculeIdx) = _moleculeData_h(cellIdx, _linkedCellNumMolecules_h(cellIdx) - 1); + _linkedCellNumMolecules_h(cellIdx) -= 1; } void simplemd::MoleculeContainer::clearLinkedCell(unsigned int cellIdx) { + synchronizeMemory(HostOperation{}); #if (MD_ERROR == MD_YES) - if (cellIdx >= _linkedCellNumMolecules.size()) { + if (cellIdx >= _linkedCellNumMolecules_h.size()) { Kokkos::printf("ERROR simplemd::MoleculeContainer::clearLinkedCell: Index out of range: %i (Local cell count = %i )\n", cellIdx, - _linkedCellNumMolecules.size()); + _linkedCellNumMolecules_h.size()); Kokkos::abort("\n"); } #endif - _linkedCellNumMolecules(cellIdx) = 0; + _linkedCellNumMolecules_h(cellIdx) = 0; } void simplemd::MoleculeContainer::sort(unsigned int cellIdx) { // set all outgoing molecules - for (size_t i = 0; i < _linkedCellNumMolecules(cellIdx); i++) { - unsigned int curMolIdx = positionToCellIndex(_moleculeData(cellIdx, i).getPosition()); + synchronizeMemory(HostOperation{}); + for (size_t i = 0; i < _linkedCellNumMolecules_h(cellIdx); i++) { + unsigned int curMolIdx = positionToCellIndex(_moleculeData_h(cellIdx, i).getPosition()); if (curMolIdx != cellIdx) { // if molecule does not belong to current cell anymore // write data to target end - _moleculeData(curMolIdx, _linkedCellNumMolecules(curMolIdx)) = _moleculeData(cellIdx, i); + _moleculeData_h(curMolIdx, _linkedCellNumMolecules_h(curMolIdx)) = _moleculeData_h(cellIdx, i); // increment target end - _linkedCellNumMolecules(curMolIdx)++; + _linkedCellNumMolecules_h(curMolIdx)++; // delete molecule at own position remove(cellIdx, i); // decrement iterator as the molecule at position i is now new @@ -128,6 +135,7 @@ void simplemd::MoleculeContainer::sort() { // parallelise loop for all cells that are to be traversed in this way printNonGhostCells(true, "host start sort"); + synchronizeMemory(DeviceOperation{}); Kokkos::parallel_for( "simplemd::MoleculeContainer::sort", Kokkos::RangePolicy(0, length), KOKKOS_CLASS_LAMBDA(const unsigned int j) { printNonGhostCells(j == 0, "device start sort"); @@ -162,16 +170,16 @@ void simplemd::MoleculeContainer::sort() { Kokkos::printf("Handle cell %u\n", index); #endif #if (MD_ERROR == MD_YES) - if (index >= _linkedCellNumMolecules.size()) { + if (index >= _linkedCellNumMolecules_d.size()) { Kokkos::abort("simplemd::MoleculeContainer::sort() out-of-bounds access to linked cell\n"); } #endif - for (size_t i = 0; i < _linkedCellNumMolecules(index); i++) { - unsigned int curMolIdx = positionToCellIndex(_moleculeData(index, i).getPosition()); + for (size_t i = 0; i < _linkedCellNumMolecules_d(index); i++) { + unsigned int curMolIdx = positionToCellIndex(_moleculeData_d(index, i).getPosition()); if (curMolIdx != index) { // if molecule does not belong to current cell anymore // write data to target end #if (MD_ERROR == MD_YES) - if (_linkedCellNumMolecules(curMolIdx) + 1 > _cellCapacity) { + if (_linkedCellNumMolecules_d(curMolIdx) + 1 > _cellCapacity) { Kokkos::printf("Cell capacity=%d would be exceeded by an operation. " "(Increase \"capacity-factor\" of the \"domain-configuration\" in the config!) " "Exiting...", @@ -179,12 +187,12 @@ void simplemd::MoleculeContainer::sort() { Kokkos::abort("simplemd::MoleculeContainer::insert\n"); } #endif - _moleculeData(curMolIdx, _linkedCellNumMolecules(curMolIdx)) = _moleculeData(index, i); + _moleculeData_d(curMolIdx, _linkedCellNumMolecules_d(curMolIdx)) = _moleculeData_d(index, i); // increment target end - _linkedCellNumMolecules(curMolIdx)++; + _linkedCellNumMolecules_d(curMolIdx)++; // delete molecule at own position - _moleculeData(index, i) = _moleculeData(index, _linkedCellNumMolecules(index) - 1); - _linkedCellNumMolecules(index) -= 1; + _moleculeData_d(index, i) = _moleculeData_d(index, _linkedCellNumMolecules_d(index) - 1); + _linkedCellNumMolecules_d(index) -= 1; // decrement iterator as the molecule at position i is now new i--; } @@ -202,23 +210,37 @@ void simplemd::MoleculeContainer::sort() { #endif } -simplemd::Molecule& simplemd::MoleculeContainer::getMoleculeAt(size_t i, size_t j) const { return _moleculeData(i, j); } +simplemd::Molecule& simplemd::MoleculeContainer::getMoleculeAt(size_t i, size_t j) const { + synchronizationCheck(); + KOKKOS_IF_ON_HOST(( + return _moleculeData_h(i, j); + )) + KOKKOS_IF_ON_DEVICE(( + return _moleculeData_d(i, j); + )) +} simplemd::LinkedCell simplemd::MoleculeContainer::operator[](const size_t idx) const { + synchronizationCheck(); #if (MD_ERROR == MD_YES) - if (idx >= _linkedCellNumMolecules.size()) { - Kokkos::printf("ERROR simplemd::MoleculeContainer::operator[]: Index out of range: %i (Local cell count = %i )\n", idx, _linkedCellNumMolecules.size()); + if (idx >= _linkedCellNumMolecules_d.size()) { + Kokkos::printf("ERROR simplemd::MoleculeContainer::operator[]: Index out of range: %i (Local cell count = %i )\n", idx, _linkedCellNumMolecules_d.size()); Kokkos::abort("\n"); } #endif - return simplemd::LinkedCell(&_moleculeData, &_linkedCellNumMolecules, idx, isGhostCell(idx)); + KOKKOS_IF_ON_HOST(( + return simplemd::LinkedCell(&_moleculeData_h, &_linkedCellNumMolecules_h, idx, isGhostCell(idx)); + )) + KOKKOS_IF_ON_DEVICE(( + return simplemd::LinkedCell(&_moleculeData_d, &_linkedCellNumMolecules_d, idx, isGhostCell(idx)); + )) } simplemd::LinkedCell simplemd::MoleculeContainer::operator[](const tarch::la::Vector cellIdx) const { return (*this)[vectorIndexToLinear(cellIdx)]; } -size_t simplemd::MoleculeContainer::getLocalNumberOfCellsScalarWithGhost() const { return _linkedCellNumMolecules.size(); } +size_t simplemd::MoleculeContainer::getLocalNumberOfCellsScalarWithGhost() const { return _linkedCellNumMolecules_d.size(); } unsigned int simplemd::MoleculeContainer::positionToCellIndex(const tarch::la::Vector& position) const { for (unsigned int d = 0; d < MD_DIM; d++) { @@ -261,6 +283,21 @@ unsigned int simplemd::MoleculeContainer::positionToCellIndex(const tarch::la::V return vectorIndexToLinear(cellVectorIndex); } +void simplemd::MoleculeContainer::synchronizationCheck() const { +#if (MD_ERROR == MD_YES) + KOKKOS_IF_ON_HOST(( + if (_memoryState == MemoryState::H_BEHIND){ + Kokkos::abort("simplemd::MoleculeContainer::synchronizationCheck: HOST ATTEMPT TO ACCESS OUTDATED MEMORY WITHOUT SYNCHRONIZATION\n"); + } + )) + KOKKOS_IF_ON_DEVICE(( + if (_memoryState == MemoryState::D_BEHIND){ + Kokkos::abort("simplemd::MoleculeContainer::synchronizationCheck: DEVICE ATTEMPT TO ACCESS OUTDATED MEMORY WITHOUT SYNCHRONIZATION\n"); + } + )) +#endif +} + size_t simplemd::MoleculeContainer::vectorIndexToLinear(const tarch::la::Vector& vectorIndex) const { unsigned int cellLinearIndex = 0; unsigned int stepSize = 1; @@ -273,9 +310,9 @@ size_t simplemd::MoleculeContainer::vectorIndexToLinear(const tarch::la::Vector< tarch::la::Vector simplemd::MoleculeContainer::getLocalCellIndexVector(const size_t cellIndex) const { #if (MD_ERROR == MD_YES) - if (cellIndex >= _linkedCellNumMolecules.size()) { + if (cellIndex >= _linkedCellNumMolecules_d.size()) { Kokkos::printf("ERROR simplemd::MoleculeContainer::getLocalCellIndexVector: Index out of range: %i (Local cell count = %i )\n", cellIndex, - _linkedCellNumMolecules.size()); + _linkedCellNumMolecules_d.size()); Kokkos::abort("\n"); } #endif @@ -310,7 +347,12 @@ bool simplemd::MoleculeContainer::isGhostCell(const size_t cellIndex) const { void simplemd::MoleculeContainer::printCellMolecules(size_t cellIndex) const { #if (MD_DUMP_CELLS != 0) - auto cellMoleculeCount = _linkedCellNumMolecules(cellIndex); + KOKKOS_IF_ON_HOST(( + auto cellMoleculeCount = _linkedCellNumMolecules_h(cellIndex); + )) + KOKKOS_IF_ON_DEVICE(( + auto cellMoleculeCount = _linkedCellNumMolecules_d(cellIndex); + )) for (size_t i = 0; i < cellMoleculeCount; i++) { Molecule& molecule = getMoleculeAt(cellIndex, i); Kokkos::printf("%u\t", cellIndex); @@ -332,7 +374,7 @@ void simplemd::MoleculeContainer::printNonGhostCells(bool shouldPrintCells, cons Kokkos::printf("=== BEGIN DUMP MOLECULE CONTAINER ===\n"); Kokkos::printf("Label: %s\n", label); Kokkos::printf("cell\tpos_x\tpos_y\tpos_z\tvel_x\tvel_y\tvel_z\tforce_x\tforce_y\tforce_z\n"); - size_t linkedCellCount = _linkedCellNumMolecules.size(); + size_t linkedCellCount = _linkedCellNumMolecules_d.size(); size_t cellsRemaining = MD_DUMP_CELLS == 0 ? linkedCellCount : std::min(MD_DUMP_CELLS, linkedCellCount); KOKKOS_IF_ON_HOST(( @@ -354,11 +396,12 @@ void simplemd::MoleculeContainer::printNonGhostCells(bool shouldPrintCells, cons #endif } -size_t simplemd::MoleculeContainer::getLocalNumberOfMoleculesWithGhost() const { +size_t simplemd::MoleculeContainer::getLocalNumberOfMoleculesWithGhost() { Kokkos::fence(); // Ensure molecule count per cell is up to date + synchronizeMemory(HostOperation{}); size_t moleculeCount = 0; - for (unsigned int i = 0; i < _linkedCellNumMolecules.size(); i++) { - moleculeCount += _linkedCellNumMolecules(i); + for (unsigned int i = 0; i < _linkedCellNumMolecules_h.size(); i++) { + moleculeCount += _linkedCellNumMolecules_h(i); } return moleculeCount; } diff --git a/simplemd/MoleculeContainer.h b/simplemd/MoleculeContainer.h index bdf0fb325..1e50a2481 100644 --- a/simplemd/MoleculeContainer.h +++ b/simplemd/MoleculeContainer.h @@ -29,6 +29,9 @@ class MoleculeContainer; */ class simplemd::MoleculeContainer { public: + + enum class MemoryState { SYNCED, H_BEHIND, D_BEHIND }; + /** * @brief Construct a new MoleculeContainer object * @@ -141,7 +144,7 @@ class simplemd::MoleculeContainer { * * @return const size_t */ - size_t getLocalNumberOfMoleculesWithGhost() const; + size_t getLocalNumberOfMoleculesWithGhost(); /** * @brief returns the index of the first (non-ghost) cell along each dimension @@ -309,6 +312,21 @@ class simplemd::MoleculeContainer { KOKKOS_FUNCTION unsigned int positionToCellIndex(const tarch::la::Vector& position) const; private: + template void synchronizeMemory(const A& a); + void synchronizationCheck() const; + + class HostOperation { + public: + static const bool IsParallel = false; + static const bool IsReadonly = false; + }; + + class DeviceOperation { + public: + static const bool IsParallel = true; + static const bool IsReadonly = false; + }; + /** * @brief Converts a 3D local linked cell index into a 1D local linked cell index. * @@ -355,14 +373,21 @@ class simplemd::MoleculeContainer { /** local index of the first cell within this domain */ const tarch::la::Vector _localIndexOfFirstCell; - Kokkos::View _moleculeData; - Kokkos::View _linkedCellNumMolecules; + Kokkos::View _moleculeData_d; + typename Kokkos::View::host_mirror_type _moleculeData_h; + Kokkos::View _linkedCellNumMolecules_d; + typename Kokkos::View::host_mirror_type _linkedCellNumMolecules_h; + Kokkos::View _linkedCellIsGhostCell; /** index offsets of all 26 neighbor cell directions */ Kokkos::View _neighborOffsets; + + /** current status of _moleculeData and _linkedCellNumMolecules */ + MemoryState _memoryState; }; template void simplemd::MoleculeContainer::iterateMolecules(A& a) { + synchronizeMemory(a); if constexpr (A::IsParallel) { iterateMoleculesParallel(a); } else { @@ -372,8 +397,8 @@ template void simplemd::MoleculeContainer::iterateMolecules(A& a) { template void simplemd::MoleculeContainer::iterateMoleculesSerial(A& a) { a.beginMoleculeIteration(); - for (unsigned int i = 0; i < _linkedCellNumMolecules.size(); i++) { - for (unsigned int j = 0; j < _linkedCellNumMolecules(i); j++) { + for (unsigned int i = 0; i < _linkedCellNumMolecules_h.size(); i++) { + for (unsigned int j = 0; j < _linkedCellNumMolecules_h(i); j++) { #if (MD_DEBUG == MD_YES) std::cout << "Handle molecule " << j << " in cell #" << i << std::endl; #endif @@ -393,7 +418,7 @@ template void simplemd::MoleculeContainer::iterateMoleculesParallel(A& "simplemd::MoleculeContainer::iterateMoleculesParallel", Kokkos::RangePolicy(0, _linkedCellIsGhostCell.size()), KOKKOS_CLASS_LAMBDA(const unsigned int i) { printNonGhostCells(i == 0, "device start iterateMoleculesParallel"); - for (unsigned int j = 0; j < _linkedCellNumMolecules(i); j++) { + for (unsigned int j = 0; j < _linkedCellNumMolecules_d(i); j++) { a.handleMolecule(getMoleculeAt(i, j)); } printNonGhostCells(i == 0, "device end iterateMoleculesParallel"); @@ -405,6 +430,7 @@ template void simplemd::MoleculeContainer::iterateMoleculesParallel(A& } template void simplemd::MoleculeContainer::iterateMoleculesWithCell(A& a) { + synchronizeMemory(a); if constexpr (A::IsParallel) { iterateMoleculesWithCellParallel(a); } else { @@ -422,8 +448,8 @@ template void simplemd::MoleculeContainer::handleCellNeighbors(A& a, M template void simplemd::MoleculeContainer::iterateMoleculesWithCellSerial(A& a) { a.beginMoleculeIteration(); - for (unsigned int i = 0; i < _linkedCellNumMolecules.size(); i++) { - for (unsigned int j = 0; j < _linkedCellNumMolecules(i); j++) { + for (unsigned int i = 0; i < _linkedCellNumMolecules_h.size(); i++) { + for (unsigned int j = 0; j < _linkedCellNumMolecules_h(i); j++) { simplemd::LinkedCell cell = (*this)[i]; Molecule& m = getMoleculeAt(i, j); a.handleMolecule(m, cell); @@ -436,7 +462,7 @@ template void simplemd::MoleculeContainer::iterateMoleculesWithCellSer template void simplemd::MoleculeContainer::iterateMoleculesWithCellParallel(A& a) { a.beginMoleculeIteration(); const unsigned int threads_per_cell = 5; - const unsigned int length = _linkedCellNumMolecules.size() * threads_per_cell; + const unsigned int length = _linkedCellNumMolecules_d.size() * threads_per_cell; Kokkos::parallel_for( "simplemd::MoleculeContainer::iterateMoleculesWithCellParallel", Kokkos::RangePolicy(0, length), KOKKOS_CLASS_LAMBDA(const unsigned int i) { @@ -444,7 +470,7 @@ template void simplemd::MoleculeContainer::iterateMoleculesWithCellPar const unsigned int cellIndex = i / threads_per_cell; simplemd::LinkedCell cell = (*this)[cellIndex]; - for (unsigned int j = i % threads_per_cell; j < _linkedCellNumMolecules(cellIndex); j+=threads_per_cell) { + for (unsigned int j = i % threads_per_cell; j < _linkedCellNumMolecules_d(cellIndex); j+=threads_per_cell) { Molecule& m = getMoleculeAt(cellIndex, j); a.handleMolecule(m,cell); handleCellNeighbors(a, m, cell); @@ -872,6 +898,7 @@ template void simplemd::MoleculeContainer::iterateCellPairsParallel(A& template void simplemd::MoleculeContainer::iterateCellPairs(A& a, const tarch::la::Vector& lowerLeftFrontCell, const tarch::la::Vector& cellRange) { + synchronizeMemory(a); if constexpr (A::IsParallel) { iterateCellPairsParallel(a, lowerLeftFrontCell, cellRange); } else { @@ -888,6 +915,7 @@ template void simplemd::MoleculeContainer::iterateCellPairs(A& a) { template void simplemd::MoleculeContainer::iterateCells(A& a, const tarch::la::Vector& lowerLeftFrontCell, const tarch::la::Vector& cellRange) { + synchronizeMemory(a); if constexpr (A::IsParallel) { iterateCellsParallel(a, lowerLeftFrontCell, cellRange); } else { @@ -897,4 +925,38 @@ void simplemd::MoleculeContainer::iterateCells(A& a, const tarch::la::Vector void simplemd::MoleculeContainer::iterateCells(A& a) { iterateCells(a, _ghostCellLayerThickness, _numLocalCellsNoGhost); } +template +void simplemd::MoleculeContainer::synchronizeMemory(const A& a){ + if(_memoryState == MemoryState::SYNCED){ + if constexpr (A::IsReadonly) + return; + if constexpr (A::IsParallel) + _memoryState = MemoryState::H_BEHIND; + else + _memoryState = MemoryState::D_BEHIND; + } else if (_memoryState == MemoryState::H_BEHIND) { + if constexpr (A::IsParallel) + return; + // Serial mapping to be executed here, in host behind state ==> D->H memory transfer required! + Kokkos::deep_copy(_moleculeData_h, _moleculeData_d); + Kokkos::deep_copy(_linkedCellNumMolecules_h, _linkedCellNumMolecules_d); + _memoryState = MemoryState::SYNCED; + if constexpr (A::IsReadonly) + return; + _memoryState = MemoryState::D_BEHIND; + } else if (_memoryState == MemoryState::D_BEHIND) { + if constexpr (!A::IsParallel) + return; + // Parallel mapping to be executed here, device behind ==> H->D memory transfer required! + Kokkos::deep_copy(_moleculeData_d, _moleculeData_h); + Kokkos::deep_copy(_linkedCellNumMolecules_d, _linkedCellNumMolecules_h); + _memoryState = MemoryState::SYNCED; + if constexpr (A::IsReadonly) + return; + _memoryState = MemoryState::H_BEHIND; + } else { + Kokkos::abort("simplemd::MoleculeContainer::synchronizeMemory: UNKNOWN MEMORY STATE\n"); + } +} + #endif // _MOLECULARDYNAMICS_MOLECULARCONTAINER_H_ diff --git a/simplemd/cell-mappings/ComputeMeanVelocityMapping.h b/simplemd/cell-mappings/ComputeMeanVelocityMapping.h index c69bacb98..9770b0a3c 100644 --- a/simplemd/cell-mappings/ComputeMeanVelocityMapping.h +++ b/simplemd/cell-mappings/ComputeMeanVelocityMapping.h @@ -75,6 +75,7 @@ class simplemd::cellmappings::ComputeMeanVelocityMapping { const unsigned int& getGlobalNumberMolecules() const { return _particleCounter; } static const bool IsParallel = false; + static const bool IsReadonly = true; private: simplemd::services::ParallelTopologyService& _parallelTopologyService; diff --git a/simplemd/cell-mappings/ComputeTemperatureMapping.h b/simplemd/cell-mappings/ComputeTemperatureMapping.h index 7a66e5a91..4f75917f6 100644 --- a/simplemd/cell-mappings/ComputeTemperatureMapping.h +++ b/simplemd/cell-mappings/ComputeTemperatureMapping.h @@ -77,6 +77,7 @@ class simplemd::cellmappings::ComputeTemperatureMapping { const double& getTemperature() const { return _temperature; } static const bool IsParallel = false; + static const bool IsReadonly = true; private: simplemd::services::ParallelTopologyService& _parallelTopologyService; diff --git a/simplemd/cell-mappings/DeleteMoleculesMapping.h b/simplemd/cell-mappings/DeleteMoleculesMapping.h index 1b6b3c9fa..29d41fc7f 100644 --- a/simplemd/cell-mappings/DeleteMoleculesMapping.h +++ b/simplemd/cell-mappings/DeleteMoleculesMapping.h @@ -30,6 +30,7 @@ class simplemd::cellmappings::DeleteMoleculesMapping { KOKKOS_FUNCTION void handleCell(LinkedCell& cell) const { cell.clear(); } static const bool IsParallel = true; + static const bool IsReadonly = false; private: }; diff --git a/simplemd/cell-mappings/LennardJonesForceMapping.h b/simplemd/cell-mappings/LennardJonesForceMapping.h index a977dc768..9165b0ee4 100644 --- a/simplemd/cell-mappings/LennardJonesForceMapping.h +++ b/simplemd/cell-mappings/LennardJonesForceMapping.h @@ -40,6 +40,7 @@ class simplemd::cellmappings::LennardJonesForceMapping { const tarch::la::Vector& position2) const; static const bool IsParallel = true; + static const bool IsReadonly = false; private: /** epsilon */ diff --git a/simplemd/cell-mappings/LennardJonesPotentialEnergyMapping.h b/simplemd/cell-mappings/LennardJonesPotentialEnergyMapping.h index 6bf728899..361bd2672 100644 --- a/simplemd/cell-mappings/LennardJonesPotentialEnergyMapping.h +++ b/simplemd/cell-mappings/LennardJonesPotentialEnergyMapping.h @@ -30,6 +30,7 @@ class simplemd::cellmappings::LennardJonesPotentialEnergyMapping { KOKKOS_FUNCTION void handleCellPair(LinkedCell& cell1, LinkedCell& cell2, const unsigned int& cellIndex1, const unsigned int& cellIndex2) const; static const bool IsParallel = true; + static const bool IsReadonly = false; private: /** epsilon */ diff --git a/simplemd/cell-mappings/ParallelBoundaryEmptyCellsMapping.h b/simplemd/cell-mappings/ParallelBoundaryEmptyCellsMapping.h index 60b510ed4..46162203d 100644 --- a/simplemd/cell-mappings/ParallelBoundaryEmptyCellsMapping.h +++ b/simplemd/cell-mappings/ParallelBoundaryEmptyCellsMapping.h @@ -35,6 +35,7 @@ class simplemd::cellmappings::ParallelBoundaryEmptyCellsMapping { } } static const bool IsParallel = false; + static const bool IsReadonly = false; private: simplemd::services::ParallelTopologyService& _parallelTopologyService; diff --git a/simplemd/cell-mappings/PeriodicAndParallelBoundaryFillCellsMapping.h b/simplemd/cell-mappings/PeriodicAndParallelBoundaryFillCellsMapping.h index 5b6d30df7..fbc1a21ec 100644 --- a/simplemd/cell-mappings/PeriodicAndParallelBoundaryFillCellsMapping.h +++ b/simplemd/cell-mappings/PeriodicAndParallelBoundaryFillCellsMapping.h @@ -39,6 +39,7 @@ class simplemd::cellmappings::PeriodicAndParallelBoundaryFillCellsMapping { void handleCell(LinkedCell& cell); static const bool IsParallel = false; + static const bool IsReadonly = false; private: simplemd::services::ParallelTopologyService& _parallelTopologyService; diff --git a/simplemd/cell-mappings/PeriodicBoundaryEmptyCellsMapping.h b/simplemd/cell-mappings/PeriodicBoundaryEmptyCellsMapping.h index 7d88e6716..7e02a0dca 100644 --- a/simplemd/cell-mappings/PeriodicBoundaryEmptyCellsMapping.h +++ b/simplemd/cell-mappings/PeriodicBoundaryEmptyCellsMapping.h @@ -37,6 +37,7 @@ class simplemd::cellmappings::PeriodicBoundaryEmptyCellsMapping { void handleCell(LinkedCell& cell); static const bool IsParallel = false; + static const bool IsReadonly = false; private: simplemd::services::ParallelTopologyService& _parallelTopologyService; diff --git a/simplemd/cell-mappings/ProfilePlotterMapping.h b/simplemd/cell-mappings/ProfilePlotterMapping.h index c6b59c269..51302a4b4 100644 --- a/simplemd/cell-mappings/ProfilePlotterMapping.h +++ b/simplemd/cell-mappings/ProfilePlotterMapping.h @@ -135,6 +135,7 @@ class simplemd::cellmappings::ProfilePlotterMapping { } static const bool IsParallel = false; + static const bool IsReadonly = true; private: const simplemd::services::ParallelTopologyService& _parallelTopologyService; diff --git a/simplemd/cell-mappings/RDFMapping.h b/simplemd/cell-mappings/RDFMapping.h index c12a34f27..b7462be8f 100644 --- a/simplemd/cell-mappings/RDFMapping.h +++ b/simplemd/cell-mappings/RDFMapping.h @@ -41,6 +41,7 @@ class simplemd::cellmappings::RDFMapping { void handleCellPair(LinkedCell& cell1, LinkedCell& cell2, const unsigned int& cellIndex1, const unsigned int& cellIndex2); static const bool IsParallel = false; + static const bool IsReadonly = true; private: const simplemd::services::ParallelTopologyService& _parallelTopologyService; diff --git a/simplemd/cell-mappings/ResetPotentialEnergyMapping.h b/simplemd/cell-mappings/ResetPotentialEnergyMapping.h index 1a463a9b6..6ccb1eac3 100644 --- a/simplemd/cell-mappings/ResetPotentialEnergyMapping.h +++ b/simplemd/cell-mappings/ResetPotentialEnergyMapping.h @@ -31,6 +31,7 @@ class simplemd::cellmappings::ResetPotentialEnergyMapping { } static const bool IsParallel = true; + static const bool IsReadonly = false; private: const double _zero; diff --git a/simplemd/cell-mappings/VaryCheckpointMapping.h b/simplemd/cell-mappings/VaryCheckpointMapping.h index ef2615ae6..5eaaf950f 100644 --- a/simplemd/cell-mappings/VaryCheckpointMapping.h +++ b/simplemd/cell-mappings/VaryCheckpointMapping.h @@ -57,6 +57,7 @@ class simplemd::cellmappings::VaryCheckpointMapping { } static const bool IsParallel = false; + static const bool IsReadonly = false; private: const double _molecularMass; diff --git a/simplemd/molecule-mappings/Adios2Writer.h b/simplemd/molecule-mappings/Adios2Writer.h index 71cf43220..59d3c863d 100644 --- a/simplemd/molecule-mappings/Adios2Writer.h +++ b/simplemd/molecule-mappings/Adios2Writer.h @@ -80,6 +80,7 @@ class simplemd::moleculemappings::Adios2Writer { void handleMolecule(Molecule& molecule); static const bool IsParallel = false; + static const bool IsReadonly = true; private: const simplemd::services::ParallelTopologyService& _parallelTopologyService; diff --git a/simplemd/molecule-mappings/ComputeMeanVelocityMapping.h b/simplemd/molecule-mappings/ComputeMeanVelocityMapping.h index 740778092..206049b33 100644 --- a/simplemd/molecule-mappings/ComputeMeanVelocityMapping.h +++ b/simplemd/molecule-mappings/ComputeMeanVelocityMapping.h @@ -48,6 +48,7 @@ class simplemd::moleculemappings::ComputeMeanVelocityMapping { } static const bool IsParallel = false; + static const bool IsReadonly = true; private: static constexpr double maxV = 1e3; diff --git a/simplemd/molecule-mappings/ConvertForcesFixedToFloatMapping.h b/simplemd/molecule-mappings/ConvertForcesFixedToFloatMapping.h index 8606593dc..eb9155f05 100644 --- a/simplemd/molecule-mappings/ConvertForcesFixedToFloatMapping.h +++ b/simplemd/molecule-mappings/ConvertForcesFixedToFloatMapping.h @@ -28,6 +28,7 @@ class simplemd::moleculemappings::ConvertForcesFixedToFloatMapping { } void endMoleculeIteration() const {} static const bool IsParallel = true; + static const bool IsReadonly = false; }; #endif // _MOLECULARDYNAMICS_MOLECULEMAPPINGS_CONVERTFORCES_MAPPING_H_ diff --git a/simplemd/molecule-mappings/ConvertForcesFloatToFixedMapping.h b/simplemd/molecule-mappings/ConvertForcesFloatToFixedMapping.h index 894258b2e..09e8961e6 100644 --- a/simplemd/molecule-mappings/ConvertForcesFloatToFixedMapping.h +++ b/simplemd/molecule-mappings/ConvertForcesFloatToFixedMapping.h @@ -35,4 +35,5 @@ class simplemd::moleculemappings::ConvertForcesFloatToFixedMapping { void endMoleculeIteration() const {} static const bool IsParallel = true; + static const bool IsReadonly = false; }; diff --git a/simplemd/molecule-mappings/InitialPositionAndForceUpdate.h b/simplemd/molecule-mappings/InitialPositionAndForceUpdate.h index 3a9063034..ec91588b5 100644 --- a/simplemd/molecule-mappings/InitialPositionAndForceUpdate.h +++ b/simplemd/molecule-mappings/InitialPositionAndForceUpdate.h @@ -31,6 +31,7 @@ class simplemd::moleculemappings::InitialPositionAndForceUpdate { molecule.setForce(_zero); } static const bool IsParallel = false; + static const bool IsReadonly = false; private: const double _dt; diff --git a/simplemd/molecule-mappings/SetMeanVelocityMapping.h b/simplemd/molecule-mappings/SetMeanVelocityMapping.h index 308d93681..1c3b9b9d3 100644 --- a/simplemd/molecule-mappings/SetMeanVelocityMapping.h +++ b/simplemd/molecule-mappings/SetMeanVelocityMapping.h @@ -33,6 +33,7 @@ class simplemd::moleculemappings::SetMeanVelocityMapping { } static const bool IsParallel = false; + static const bool IsReadonly = false; private: const tarch::la::Vector _oldVelocity; diff --git a/simplemd/molecule-mappings/VTKMoleculeWriter.cpp b/simplemd/molecule-mappings/VTKMoleculeWriter.cpp index d9c4e2b3e..3d6c4592b 100644 --- a/simplemd/molecule-mappings/VTKMoleculeWriter.cpp +++ b/simplemd/molecule-mappings/VTKMoleculeWriter.cpp @@ -5,7 +5,7 @@ #include "simplemd/molecule-mappings/VTKMoleculeWriter.h" simplemd::moleculemappings::VTKMoleculeWriter::VTKMoleculeWriter(const simplemd::services::ParallelTopologyService& parallelTopologyService, - const simplemd::MoleculeContainer& moleculeContainer, const std::string& filename) + simplemd::MoleculeContainer& moleculeContainer, const std::string& filename) : _parallelTopologyService(parallelTopologyService), _moleculeContainer(moleculeContainer), _filename(filename), _timestep(0) {} simplemd::moleculemappings::VTKMoleculeWriter::~VTKMoleculeWriter() {} diff --git a/simplemd/molecule-mappings/VTKMoleculeWriter.h b/simplemd/molecule-mappings/VTKMoleculeWriter.h index 8ae7b1776..adbdaef1c 100644 --- a/simplemd/molecule-mappings/VTKMoleculeWriter.h +++ b/simplemd/molecule-mappings/VTKMoleculeWriter.h @@ -33,7 +33,7 @@ class simplemd::moleculemappings::VTKMoleculeWriter { * @param moleculeContainer molecule container * @param filename filename for .vtk output file */ - VTKMoleculeWriter(const simplemd::services::ParallelTopologyService& parallelTopologyService, const simplemd::MoleculeContainer& moleculeContainer, + VTKMoleculeWriter(const simplemd::services::ParallelTopologyService& parallelTopologyService, simplemd::MoleculeContainer& moleculeContainer, const std::string& filename); /** Destructor */ @@ -56,10 +56,11 @@ class simplemd::moleculemappings::VTKMoleculeWriter { void handleMolecule(Molecule& molecule); static const bool IsParallel = false; + static const bool IsReadonly = true; private: const simplemd::services::ParallelTopologyService& _parallelTopologyService; - const simplemd::MoleculeContainer& _moleculeContainer; + simplemd::MoleculeContainer& _moleculeContainer; /** filename */ std::string _filename; /** current timestep */ diff --git a/simplemd/molecule-mappings/VelocityStoermerVerletMapping.h b/simplemd/molecule-mappings/VelocityStoermerVerletMapping.h index 24ae9a66b..a116b6cce 100644 --- a/simplemd/molecule-mappings/VelocityStoermerVerletMapping.h +++ b/simplemd/molecule-mappings/VelocityStoermerVerletMapping.h @@ -31,6 +31,7 @@ class simplemd::moleculemappings::VelocityStoermerVerletMapping { void handleMolecule(Molecule& molecule); static const bool IsParallel = false; + static const bool IsReadonly = false; private: tarch::la::Vector<2 * MD_DIM, bool> initReflectingBoundary(const tarch::la::Vector& boundary) const; diff --git a/simplemd/molecule-mappings/WriteCheckPointMapping.h b/simplemd/molecule-mappings/WriteCheckPointMapping.h index cc9ab0d7f..9458fb816 100644 --- a/simplemd/molecule-mappings/WriteCheckPointMapping.h +++ b/simplemd/molecule-mappings/WriteCheckPointMapping.h @@ -34,6 +34,7 @@ class simplemd::moleculemappings::WriteCheckPointMapping { void handleMolecule(Molecule& molecule); static const bool IsParallel = false; + static const bool IsReadonly = true; private: const simplemd::services::ParallelTopologyService& _parallelTopologyService; diff --git a/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.h b/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.h index 1440da170..85f86c7b2 100644 --- a/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.h +++ b/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.h @@ -22,5 +22,6 @@ class simplemd::moleculewithcellmappings::LennardJonesForceMapping : public simp KOKKOS_FUNCTION void handleMolecule(Molecule& molecule, const LinkedCell& cell) const; static const bool IsParallel = true; + static const bool IsReadonly = false; }; diff --git a/test/integration/SimpleMDBench.h b/test/integration/SimpleMDBench.h index 59ca6aaed..b6374aec5 100644 --- a/test/integration/SimpleMDBench.h +++ b/test/integration/SimpleMDBench.h @@ -30,6 +30,7 @@ class BenchSim : public simplemd::MolecularDynamicsSimulation { unsigned long long checksum() const { return sum; } static const bool IsParallel = false; + static const bool IsReadonly = true; private: void process(const double& data) { sum ^= *((unsigned long long*)&data); } From 4895a8b86dfa6c993aa21b629745da44476b5f8a Mon Sep 17 00:00:00 2001 From: Piet Jarmatz Date: Tue, 11 Aug 2026 13:40:32 +0200 Subject: [PATCH 02/16] disabled synchronizationCheck on host-only code --- simplemd/MoleculeContainer.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/simplemd/MoleculeContainer.cpp b/simplemd/MoleculeContainer.cpp index 37244c903..62861ed78 100644 --- a/simplemd/MoleculeContainer.cpp +++ b/simplemd/MoleculeContainer.cpp @@ -285,6 +285,7 @@ unsigned int simplemd::MoleculeContainer::positionToCellIndex(const tarch::la::V void simplemd::MoleculeContainer::synchronizationCheck() const { #if (MD_ERROR == MD_YES) +#if defined(KOKKOS_TARGET_CUDA) KOKKOS_IF_ON_HOST(( if (_memoryState == MemoryState::H_BEHIND){ Kokkos::abort("simplemd::MoleculeContainer::synchronizationCheck: HOST ATTEMPT TO ACCESS OUTDATED MEMORY WITHOUT SYNCHRONIZATION\n"); @@ -296,6 +297,7 @@ void simplemd::MoleculeContainer::synchronizationCheck() const { } )) #endif +#endif } size_t simplemd::MoleculeContainer::vectorIndexToLinear(const tarch::la::Vector& vectorIndex) const { From d4ce71f3d151102cb755fbb6aba09c7470305661 Mon Sep 17 00:00:00 2001 From: Piet Jarmatz Date: Tue, 11 Aug 2026 14:20:10 +0200 Subject: [PATCH 03/16] MoleculeContainer::synchronizationCheck() is a KOKKOS_FUNCTION --- simplemd/MoleculeContainer.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/simplemd/MoleculeContainer.h b/simplemd/MoleculeContainer.h index 1e50a2481..4efb1acd1 100644 --- a/simplemd/MoleculeContainer.h +++ b/simplemd/MoleculeContainer.h @@ -313,7 +313,7 @@ class simplemd::MoleculeContainer { private: template void synchronizeMemory(const A& a); - void synchronizationCheck() const; + KOKKOS_FUNCTION void synchronizationCheck() const; class HostOperation { public: From 050c77c9032d6c53cf048d44bc174336b6967302 Mon Sep 17 00:00:00 2001 From: Piet Jarmatz Date: Tue, 11 Aug 2026 17:05:15 +0200 Subject: [PATCH 04/16] Update MoleculeContainerTest mappings --- test/unit/simplemd/MoleculeContainerTest.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/unit/simplemd/MoleculeContainerTest.cpp b/test/unit/simplemd/MoleculeContainerTest.cpp index da157dcb8..7d532bae9 100644 --- a/test/unit/simplemd/MoleculeContainerTest.cpp +++ b/test/unit/simplemd/MoleculeContainerTest.cpp @@ -22,6 +22,7 @@ template class MolPosIncrMapping { } } static const bool IsParallel = par; + static const bool IsReadonly = false; }; template class CellMolPosIncrMapping { @@ -38,6 +39,7 @@ template class CellMolPosIncrMapping { } } static const bool IsParallel = par; + static const bool IsReadonly = false; }; template class CellPairPosMaxMapping { @@ -62,6 +64,7 @@ template class CellPairPosMaxMapping { } } static const bool IsParallel = par; + static const bool IsReadonly = false; }; class MoleculeContainerTest : public CppUnit::TestFixture { From 7dcb11cd833bcc632f55c0eeebbb6bba731ce790 Mon Sep 17 00:00:00 2001 From: Piet Jarmatz Date: Tue, 11 Aug 2026 17:06:31 +0200 Subject: [PATCH 05/16] VelocityStoermerVerletMapping is now parallel --- simplemd/Molecule.h | 16 ++++++------ .../VelocityStoermerVerletMapping.cpp | 26 +++++++++---------- .../VelocityStoermerVerletMapping.h | 6 ++--- 3 files changed, 23 insertions(+), 25 deletions(-) diff --git a/simplemd/Molecule.h b/simplemd/Molecule.h index c82c6c02d..bbcc39522 100644 --- a/simplemd/Molecule.h +++ b/simplemd/Molecule.h @@ -36,27 +36,27 @@ class simplemd::Molecule { /** get/set ID */ KOKKOS_FUNCTION const unsigned int& getID() const { return _id; } - void setID(const unsigned int& id) { _id = id; } + KOKKOS_FUNCTION void setID(const unsigned int& id) { _id = id; } /** get/ set position */ KOKKOS_FUNCTION tarch::la::Vector& getPosition() { return _position; } KOKKOS_FUNCTION const tarch::la::Vector& getConstPosition() const { return _position; } - void setPosition(const tarch::la::Vector& position) { _position = position; } + KOKKOS_FUNCTION void setPosition(const tarch::la::Vector& position) { _position = position; } /** get/ set velocity */ - tarch::la::Vector& getVelocity() { return _velocity; } + KOKKOS_FUNCTION tarch::la::Vector& getVelocity() { return _velocity; } KOKKOS_FUNCTION const tarch::la::Vector& getConstVelocity() const { return _velocity; } - void setVelocity(const tarch::la::Vector& velocity) { _velocity = velocity; } + KOKKOS_FUNCTION void setVelocity(const tarch::la::Vector& velocity) { _velocity = velocity; } /** get/ set force */ KOKKOS_FUNCTION tarch::la::Vector& getForce() { return _force; } KOKKOS_FUNCTION const tarch::la::Vector& getConstForce() const { return _force; } - void setForce(const tarch::la::Vector& force) { _force = force; } + KOKKOS_FUNCTION void setForce(const tarch::la::Vector& force) { _force = force; } /** get/ set force of last timestep */ - tarch::la::Vector& getForceOld() { return _forceOld; } - const tarch::la::Vector& getConstForceOld() const { return _forceOld; } - void setForceOld(const tarch::la::Vector& force) { _forceOld = force; } + KOKKOS_FUNCTION tarch::la::Vector& getForceOld() { return _forceOld; } + KOKKOS_FUNCTION const tarch::la::Vector& getConstForceOld() const { return _forceOld; } + KOKKOS_FUNCTION void setForceOld(const tarch::la::Vector& force) { _forceOld = force; } KOKKOS_FUNCTION double& getPotentialEnergy() { return _potentialEnergy; } KOKKOS_FUNCTION const double& getConstPotentialEnergy() const { return _potentialEnergy; } diff --git a/simplemd/molecule-mappings/VelocityStoermerVerletMapping.cpp b/simplemd/molecule-mappings/VelocityStoermerVerletMapping.cpp index 406b67e1f..da44bb34b 100644 --- a/simplemd/molecule-mappings/VelocityStoermerVerletMapping.cpp +++ b/simplemd/molecule-mappings/VelocityStoermerVerletMapping.cpp @@ -15,7 +15,7 @@ void simplemd::moleculemappings::VelocityStoermerVerletMapping::beginMoleculeIte void simplemd::moleculemappings::VelocityStoermerVerletMapping::endMoleculeIteration() {} -void simplemd::moleculemappings::VelocityStoermerVerletMapping::handleMolecule(Molecule& molecule) { +void simplemd::moleculemappings::VelocityStoermerVerletMapping::handleMolecule(Molecule& molecule) const { // if the molecule is fixed in space, return immediately: if (molecule.isFixed()) return; @@ -35,21 +35,19 @@ void simplemd::moleculemappings::VelocityStoermerVerletMapping::handleMolecule(M #if (MD_ERROR == MD_YES) for (unsigned int d = 0; d < MD_DIM; d++) { if (std::isnan(position[d]) || std::isinf(position[d])) { - std::cout << "ERROR simplemd::moleculemappings::VelocityStoermerVerletMapping::handleMolecule: Position "; - std::cout << d << " is out of range" << std::endl; - std::cout << "Position: " << position << ", molecule: " << molecule.getID() << std::endl; - std::cout << "Velocity: " << velocity << ", molecule: " << molecule.getID() << std::endl; - std::cout << "OldVelocity: " << oldVelocity << ", molecule: " << molecule.getID() << std::endl; - std::cout << "Force: " << molecule.getConstForce() << ", old: " << molecule.getConstForceOld() << std::endl; - std::cout << "Old position: " << oldPosition << std::endl; - exit(EXIT_FAILURE); + Kokkos::printf("simplemd::moleculemappings::VelocityStoermerVerletMapping::handleMolecule: Position " + "%u is out of range!\n", d); + Kokkos::printf("Force: %lf %lf %lf; " + "Position: %lf %lf %lf; " + "ID: %u;" + "\n", + molecule.getConstForce()[0], molecule.getConstForce()[1], molecule.getConstForce()[2], + position[0], position[1], position[2], + molecule.getID()); + Kokkos::abort("ERROR Position out of range"); } if (std::isnan(velocity[d]) || std::isinf(velocity[d])) { - std::cout << "ERROR simplemd::moleculemappings::VelocityStoermerVerletMapping::handleMolecule: Velocity "; - std::cout << d << " is NaN or Inf" << std::endl; - std::cout << velocity << std::endl; - std::cout << molecule.getConstForce() << ", " << molecule.getConstForceOld() << std::endl; - exit(EXIT_FAILURE); + Kokkos::abort("ERROR simplemd::moleculemappings::VelocityStoermerVerletMapping::handleMolecule: Velocity is NaN or Inf"); } } #endif diff --git a/simplemd/molecule-mappings/VelocityStoermerVerletMapping.h b/simplemd/molecule-mappings/VelocityStoermerVerletMapping.h index a116b6cce..117d2d38d 100644 --- a/simplemd/molecule-mappings/VelocityStoermerVerletMapping.h +++ b/simplemd/molecule-mappings/VelocityStoermerVerletMapping.h @@ -23,14 +23,14 @@ class simplemd::moleculemappings::VelocityStoermerVerletMapping { VelocityStoermerVerletMapping(const double& kB, const double& dt, const double& mass, const tarch::la::Vector& boundary, const tarch::la::Vector& domainOffset, const tarch::la::Vector& domainSize); - ~VelocityStoermerVerletMapping(); + KOKKOS_FUNCTION ~VelocityStoermerVerletMapping(); void beginMoleculeIteration(); void endMoleculeIteration(); - void handleMolecule(Molecule& molecule); + KOKKOS_FUNCTION void handleMolecule(Molecule& molecule) const; - static const bool IsParallel = false; + static const bool IsParallel = true; static const bool IsReadonly = false; private: From d8e29684b702814b24c7f18bc79172c69822e31a Mon Sep 17 00:00:00 2001 From: Piet Jarmatz Date: Tue, 11 Aug 2026 17:07:09 +0200 Subject: [PATCH 06/16] Add missing Kokkos::finalize() in SimpleMDBench --- test/integration/SimpleMDBench.h | 1 + 1 file changed, 1 insertion(+) diff --git a/test/integration/SimpleMDBench.h b/test/integration/SimpleMDBench.h index b6374aec5..1d30cd25c 100644 --- a/test/integration/SimpleMDBench.h +++ b/test/integration/SimpleMDBench.h @@ -179,6 +179,7 @@ class SimpleMDBench : public Test { std::cout << "ERROR SimpleMDBench: ERROR Checksum is wrong!! " << std::endl; std::cout << "(Note that this is expected and OK in Release build mode,)" << std::endl; std::cout << "(as checksum is supposed to match only in Debug / DebugOptimized build mode.)" << std::endl; + Kokkos::finalize(); exit(EXIT_FAILURE); } } From 9a85de9363e3cb4f2a5c208cf63f1c95d4968bd0 Mon Sep 17 00:00:00 2001 From: Piet Jarmatz Date: Wed, 12 Aug 2026 15:41:21 +0200 Subject: [PATCH 07/16] Lots if Kokkos::fence() removed --- simplemd/MoleculeContainer.cpp | 2 -- simplemd/MoleculeContainer.h | 14 ++++++++++---- test/integration/SimpleMDBench.h | 6 ++++-- 3 files changed, 14 insertions(+), 8 deletions(-) diff --git a/simplemd/MoleculeContainer.cpp b/simplemd/MoleculeContainer.cpp index 62861ed78..7cae25df1 100644 --- a/simplemd/MoleculeContainer.cpp +++ b/simplemd/MoleculeContainer.cpp @@ -199,7 +199,6 @@ void simplemd::MoleculeContainer::sort() { } printNonGhostCells(j == 0, "device end sorf"); }); // j, Kokkos::parallel_for - Kokkos::fence(); // Ensure results are available on the host printNonGhostCells(true, "host end sort"); } // x #if (MD_DIM > 1) @@ -399,7 +398,6 @@ void simplemd::MoleculeContainer::printNonGhostCells(bool shouldPrintCells, cons } size_t simplemd::MoleculeContainer::getLocalNumberOfMoleculesWithGhost() { - Kokkos::fence(); // Ensure molecule count per cell is up to date synchronizeMemory(HostOperation{}); size_t moleculeCount = 0; for (unsigned int i = 0; i < _linkedCellNumMolecules_h.size(); i++) { diff --git a/simplemd/MoleculeContainer.h b/simplemd/MoleculeContainer.h index 4efb1acd1..c5247d72c 100644 --- a/simplemd/MoleculeContainer.h +++ b/simplemd/MoleculeContainer.h @@ -2,6 +2,7 @@ #define _MOLECULARDYNAMICS_MOLECULARCONTAINER_H_ #include +#include #include #include @@ -9,6 +10,8 @@ #include #include +#define MEMSYNC_DEBUG + namespace simplemd { namespace services { // forward declarations to remove circular dependencies @@ -423,7 +426,6 @@ template void simplemd::MoleculeContainer::iterateMoleculesParallel(A& } printNonGhostCells(i == 0, "device end iterateMoleculesParallel"); }); - Kokkos::fence(); // Ensure results are available on the host printNonGhostCells(true, "host end iterateMoleculesParallel"); a.endMoleculeIteration(); #endif @@ -476,7 +478,6 @@ template void simplemd::MoleculeContainer::iterateMoleculesWithCellPar handleCellNeighbors(a, m, cell); } }); - Kokkos::fence(); // Ensure results are available on the host a.endMoleculeIteration(); @@ -716,7 +717,6 @@ void simplemd::MoleculeContainer::iterateCellsParallel(A& a, const tarch::la::Ve a.handleCell(cell); printNonGhostCells(i == 0, "device end iterateCellsParallel"); }); // Kokkos::parallel_for - Kokkos::fence(); // Ensure results are available on the host printNonGhostCells(true, "host end iterateCellsParallel"); // end iteration(); a.endCellIteration(); @@ -876,7 +876,6 @@ void simplemd::MoleculeContainer::iterateCellPairsParallel(A& a, const tarch::la } printNonGhostCells(j == 0, "device end iterateCellPairsParallel"); }); // j, Kokkos::parallel_for - Kokkos::fence(); // Ensure results are available on the host printNonGhostCells(true, "host end iterateCellPairsParallel"); } // x #if (MD_DIM > 1) @@ -938,8 +937,12 @@ void simplemd::MoleculeContainer::synchronizeMemory(const A& a){ if constexpr (A::IsParallel) return; // Serial mapping to be executed here, in host behind state ==> D->H memory transfer required! + Kokkos::fence(); Kokkos::deep_copy(_moleculeData_h, _moleculeData_d); Kokkos::deep_copy(_linkedCellNumMolecules_h, _linkedCellNumMolecules_d); +#ifdef MEMSYNC_DEBUG + std::cout << "D->H memory transfer caused by " << typeid(A).name() << std::endl; +#endif _memoryState = MemoryState::SYNCED; if constexpr (A::IsReadonly) return; @@ -950,6 +953,9 @@ void simplemd::MoleculeContainer::synchronizeMemory(const A& a){ // Parallel mapping to be executed here, device behind ==> H->D memory transfer required! Kokkos::deep_copy(_moleculeData_d, _moleculeData_h); Kokkos::deep_copy(_linkedCellNumMolecules_d, _linkedCellNumMolecules_h); +#ifdef MEMSYNC_DEBUG + std::cout << "H->D memory transfer caused by " << typeid(A).name() << std::endl; +#endif _memoryState = MemoryState::SYNCED; if constexpr (A::IsReadonly) return; diff --git a/test/integration/SimpleMDBench.h b/test/integration/SimpleMDBench.h index 1d30cd25c..c56f332d8 100644 --- a/test/integration/SimpleMDBench.h +++ b/test/integration/SimpleMDBench.h @@ -129,7 +129,8 @@ class SimpleMDBench : public Test { _simulation = std::make_unique(_simpleMDConfig); _simulation->initServices(); - std::cout << "INFO SimpleMDBench: Initial Checksum is " << _simulation->getChecksum() << std::endl; + unsigned long long sum = _simulation->getChecksum(); + std::cout << "INFO SimpleMDBench: Initial Checksum is " << sum << std::endl; } void shutdown() { _simulation->shutdownServices(); } @@ -137,7 +138,8 @@ class SimpleMDBench : public Test { void bench() { // warm-up timestep, for more reliable benchmarking result _simulation->simulateOneTimestep(0); - std::cout << "INFO SimpleMDBench: Warmup Checksum is " << _simulation->getChecksum() << std::endl; + unsigned long long sum = _simulation->getChecksum(); + std::cout << "INFO SimpleMDBench: Warmup Checksum is " << sum << std::endl; timeval start, end; gettimeofday(&start, NULL); From 5c87c4cf140768ec10ddf7264119844f6f54235b Mon Sep 17 00:00:00 2001 From: Ruben Horn <> Date: Wed, 12 Aug 2026 15:47:40 +0200 Subject: [PATCH 08/16] Few fixes for nvcc --- CMakeLists.txt | 24 ++++++++++++++----- .../LennardJonesForceMapping.cpp | 2 +- .../ComputeMeanVelocityMapping.h | 2 +- .../ConvertForcesFixedToFloatMapping.h | 2 +- .../ConvertForcesFloatToFixedMapping.h | 2 +- 5 files changed, 22 insertions(+), 10 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index b6ec196f9..172f4537b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -100,6 +100,22 @@ message(STATUS ">> Build with PINT_DEBUG : " ${BUILD_WITH_PINT_DEBUG}) option(BUILD_WITH_ADIOS2 "Enable output in the Adios2 format." OFF) message(STATUS ">> Build with Adios2: " ${BUILD_WITH_ADIOS2}) +# Profiling +if(CMAKE_BUILD_TYPE STREQUAL "RelWithDebInfo" OR (CMAKE_BUILD_TYPE MATCHES "^Debug")) + # Support source code view in ncu-ui + add_compile_options( + $<$:-lineinfo> + ) +endif() + +# Force -O3 optimization in RelWithDebInfo builds +# (Avoids multiple optimization levels in the compiler flags) +if(CMAKE_BUILD_TYPE STREQUAL "RelWithDebInfo") + string(REPLACE "-O2" "-O3" + CMAKE_CXX_FLAGS_RELWITHDEBINFO + "${CMAKE_CXX_FLAGS_RELWITHDEBINFO}") +endif() + # # Find dependencies # @@ -146,10 +162,6 @@ else() # Debugging if(CMAKE_BUILD_TYPE MATCHES "^Debug" AND (NOT CMAKE_BUILD_TYPE STREQUAL "DebugOptimized")) set(KOKKOS_DEBUG "ON") - # Support source code view in ncu-ui - add_compile_options( - $<$:-lineinfo> - ) else() set(KOKKOS_DEBUG "OFF") endif() @@ -309,7 +321,7 @@ elseif(CMAKE_BUILD_TYPE STREQUAL "DebugOptimized") target_compile_options(simplemd_interface INTERFACE --fmad=false) endif() elseif(CMAKE_BUILD_TYPE STREQUAL "RelWithDebInfo") - target_compile_options(simplemd_interface INTERFACE -O3 -g3) + target_compile_options(simplemd_interface INTERFACE -g3) elseif(CMAKE_BUILD_TYPE STREQUAL "Release") if(KOKKOS_TARGET STREQUAL "CUDA") target_compile_options(simplemd_interface INTERFACE -O3) @@ -526,7 +538,7 @@ elseif(CMAKE_BUILD_TYPE STREQUAL "DebugOptimized") target_compile_definitions(mamico INTERFACE MDCoupledError) target_compile_options(mamico INTERFACE -O3 -g3) elseif(CMAKE_BUILD_TYPE STREQUAL "RelWithDebInfo") - target_compile_options(mamico INTERFACE -O3 -g3) + target_compile_options(mamico INTERFACE -g3) elseif(CMAKE_BUILD_TYPE STREQUAL "Release") if(KOKKOS_TARGET STREQUAL "CUDA") target_compile_options(mamico INTERFACE -O3) diff --git a/simplemd/cell-mappings/LennardJonesForceMapping.cpp b/simplemd/cell-mappings/LennardJonesForceMapping.cpp index 4df427c62..43c62dc48 100644 --- a/simplemd/cell-mappings/LennardJonesForceMapping.cpp +++ b/simplemd/cell-mappings/LennardJonesForceMapping.cpp @@ -139,7 +139,7 @@ simplemd::cellmappings::LennardJonesForceMapping::getLennardJonesForce(const tar #if (TARCH_DEBUG == TARCH_YES) tarch::la::Vector res{24.0 * _epsilon / rij2 * (_sigma6 / rij6) * (1.0 - 2.0 * (_sigma6 / rij6)) * rij}; constexpr double maxF = 1e6; - constexpr double stepF = (double)(std::numeric_limits::max()) / maxF; + constexpr double stepF = static_cast(LLONG_MAX) / maxF; res = stepF * res; long long fb0{(long long)(res[0])}; long long fb1{(long long)(res[1])}; diff --git a/simplemd/molecule-mappings/ComputeMeanVelocityMapping.h b/simplemd/molecule-mappings/ComputeMeanVelocityMapping.h index 206049b33..779e7020e 100644 --- a/simplemd/molecule-mappings/ComputeMeanVelocityMapping.h +++ b/simplemd/molecule-mappings/ComputeMeanVelocityMapping.h @@ -52,7 +52,7 @@ class simplemd::moleculemappings::ComputeMeanVelocityMapping { private: static constexpr double maxV = 1e3; - static constexpr double stepV = (double)(std::numeric_limits::max()) / maxV; + static constexpr double stepV = static_cast(LLONG_MAX) / maxV; static constexpr double minV = 1 / stepV; tarch::la::Vector _meanVelocity; long long _particleCounter; diff --git a/simplemd/molecule-mappings/ConvertForcesFixedToFloatMapping.h b/simplemd/molecule-mappings/ConvertForcesFixedToFloatMapping.h index eb9155f05..6d861061b 100644 --- a/simplemd/molecule-mappings/ConvertForcesFixedToFloatMapping.h +++ b/simplemd/molecule-mappings/ConvertForcesFixedToFloatMapping.h @@ -18,7 +18,7 @@ class simplemd::moleculemappings::ConvertForcesFixedToFloatMapping { void beginMoleculeIteration() const {} KOKKOS_FUNCTION void handleMolecule(simplemd::Molecule& molecule) const { constexpr double maxF = 1e6; - constexpr double stepF = (double)(std::numeric_limits::max()) / maxF; + constexpr double stepF = static_cast(LLONG_MAX)/ maxF; constexpr double minF = 1 / stepF; tarch::la::Vector& force = molecule.getForce(); diff --git a/simplemd/molecule-mappings/ConvertForcesFloatToFixedMapping.h b/simplemd/molecule-mappings/ConvertForcesFloatToFixedMapping.h index 09e8961e6..8a4407726 100644 --- a/simplemd/molecule-mappings/ConvertForcesFloatToFixedMapping.h +++ b/simplemd/molecule-mappings/ConvertForcesFloatToFixedMapping.h @@ -18,7 +18,7 @@ class simplemd::moleculemappings::ConvertForcesFloatToFixedMapping { KOKKOS_FUNCTION void handleMolecule(simplemd::Molecule& molecule) const { constexpr double maxF = 1e6; - constexpr double stepF = (double)(std::numeric_limits::max()) / maxF; + constexpr double stepF = static_cast(LLONG_MAX) / maxF; tarch::la::Vector force = stepF * molecule.getForce(); long long fb0{(long long)(force[0])}; From 8d8a0b24376aa9c9dd220de0ecbe5c8586464129 Mon Sep 17 00:00:00 2001 From: Piet Jarmatz Date: Wed, 12 Aug 2026 17:38:08 +0200 Subject: [PATCH 09/16] Added MD120 and MD240 options to SimpleMDBench --- simplemd/configurations/DomainConfiguration.h | 2 ++ test/integration/SimpleMDBench.h | 22 +++++++++++++++++-- test/integration/benchmarks.cpp | 17 +++++++++++--- 3 files changed, 36 insertions(+), 5 deletions(-) diff --git a/simplemd/configurations/DomainConfiguration.h b/simplemd/configurations/DomainConfiguration.h index 5d8f8ee75..2b300124a 100644 --- a/simplemd/configurations/DomainConfiguration.h +++ b/simplemd/configurations/DomainConfiguration.h @@ -48,7 +48,9 @@ class simplemd::configurations::DomainConfiguration : public tarch::configuratio /** getters for all parsed and computed quantities */ const tarch::la::Vector& getMoleculesPerDirection() const { return _moleculesPerDirection; } + tarch::la::Vector& getMoleculesPerDirectionNonConst() { return _moleculesPerDirection; } const tarch::la::Vector& getGlobalDomainSize() const { return _domainSize; } + tarch::la::Vector& getGlobalDomainSizeNonConst() { return _domainSize; } const tarch::la::Vector& getGlobalDomainOffset() const { return _domainOffset; } const double& getCutoffRadius() const { return _cutoffRadius; } const tarch::la::Vector& getMeshWidth() const { return _meshWidth; } diff --git a/test/integration/SimpleMDBench.h b/test/integration/SimpleMDBench.h index c56f332d8..e21f7e9fd 100644 --- a/test/integration/SimpleMDBench.h +++ b/test/integration/SimpleMDBench.h @@ -11,6 +11,8 @@ #include #include +enum class SimpleMDBenchSize { MD60, MD120, MD240 }; + class BenchSim : public simplemd::MolecularDynamicsSimulation { public: BenchSim(const simplemd::configurations::MolecularDynamicsConfiguration& configuration) : simplemd::MolecularDynamicsSimulation(configuration) {} @@ -48,7 +50,8 @@ class BenchSim : public simplemd::MolecularDynamicsSimulation { class SimpleMDBench : public Test { public: - SimpleMDBench() : Test("SimpleMDBench") {} + SimpleMDBench(SimpleMDBenchSize benchsize) : + Test("SimpleMDBench"), _MDSize(benchsize) {} virtual ~SimpleMDBench() {} virtual void run() { @@ -126,6 +129,17 @@ class SimpleMDBench : public Test { std::cout << "ERROR SimpleMDBench: Invalid SimpleMD config!" << std::endl; exit(EXIT_FAILURE); } + if(_MDSize == SimpleMDBenchSize::MD120){ + _simpleMDConfig.getDomainConfigurationNonConst().getMoleculesPerDirectionNonConst() + = tarch::la::Vector<3, unsigned int>{112,112,112}; + _simpleMDConfig.getDomainConfigurationNonConst().getGlobalDomainSizeNonConst() + = tarch::la::Vector<3, double>{120,120,120}; + } else if(_MDSize == SimpleMDBenchSize::MD240){ + _simpleMDConfig.getDomainConfigurationNonConst().getMoleculesPerDirectionNonConst() + = tarch::la::Vector<3, unsigned int>{224,224,224}; + _simpleMDConfig.getDomainConfigurationNonConst().getGlobalDomainSizeNonConst() + = tarch::la::Vector<3, double>{240,240,240}; + } _simulation = std::make_unique(_simpleMDConfig); _simulation->initServices(); @@ -174,7 +188,10 @@ class SimpleMDBench : public Test { unsigned long long sum = _simulation->getChecksum(); std::cout << "INFO SimpleMDBench: Final XOR Checksum is " << sum << std::endl; - unsigned long long correct = 9224833479527670225u; + unsigned long long correct; + if(_MDSize == SimpleMDBenchSize::MD60) correct = 9224833479527670225u; + if(_MDSize == SimpleMDBenchSize::MD120) correct = 18377339084083654394u; + if(_MDSize == SimpleMDBenchSize::MD240) correct = 18342705775936782749u; if (sum == correct) std::cout << "INFO SimpleMDBench: SUCCESS Checksum is correct :-)" << std::endl; else { @@ -189,4 +206,5 @@ class SimpleMDBench : public Test { int _rank; simplemd::configurations::MolecularDynamicsConfiguration _simpleMDConfig; std::unique_ptr _simulation; + SimpleMDBenchSize _MDSize; }; diff --git a/test/integration/benchmarks.cpp b/test/integration/benchmarks.cpp index d888c9a50..7eca9d99d 100644 --- a/test/integration/benchmarks.cpp +++ b/test/integration/benchmarks.cpp @@ -30,10 +30,21 @@ int main(int argc, char* argv[]) { std::cout << "Available concurrency: " << mainExecSpace.concurrency() << std::endl; + SimpleMDBenchSize benchsize; // run tests - runTest(new CellIdxIterBench()); - std::cout << std::endl << "==================== ==================== ====================" << std::endl << std::endl; - runTest(new SimpleMDBench()); + if(argc == 1){ + runTest(new CellIdxIterBench()); + std::cout << std::endl << "==================== ==================== ====================" << std::endl << std::endl; + benchsize = SimpleMDBenchSize::MD60; + } else if(argc == 2 && !strcmp(argv[1],"MD60")) benchsize = SimpleMDBenchSize::MD60; + else if(argc == 2 && !strcmp(argv[1],"MD120")) benchsize = SimpleMDBenchSize::MD120; + else if(argc == 2 && !strcmp(argv[1],"MD240")) benchsize = SimpleMDBenchSize::MD240; + else { + std::cout << "ERROR unknown parameter" << std::endl; + Kokkos::finalize(); + exit(EXIT_FAILURE); + } + runTest(new SimpleMDBench(benchsize)); } #if (COUPLING_MD_PARALLEL == COUPLING_MD_YES) MPI_Finalize(); From 10f982c3bceedf4c6bec2e27e34ebc96e50ada5e Mon Sep 17 00:00:00 2001 From: Piet Jarmatz Date: Wed, 12 Aug 2026 18:18:48 +0200 Subject: [PATCH 10/16] Added missing timing fence in SimpleMDBench --- test/integration/SimpleMDBench.h | 1 + 1 file changed, 1 insertion(+) diff --git a/test/integration/SimpleMDBench.h b/test/integration/SimpleMDBench.h index e21f7e9fd..8477bbfcb 100644 --- a/test/integration/SimpleMDBench.h +++ b/test/integration/SimpleMDBench.h @@ -160,6 +160,7 @@ class SimpleMDBench : public Test { for (unsigned int t = 1; t < _simpleMDConfig.getSimulationConfiguration().getNumberOfTimesteps() + 1; t++) { _simulation->simulateOneTimestep(t); } + Kokkos::fence(); gettimeofday(&end, NULL); double runtime = (end.tv_sec - start.tv_sec) * 1000000 + (end.tv_usec - start.tv_usec); std::cout << "Runtime: " << (int)(runtime / 1000) << "ms" << std::endl; From 52761dc3ff3f2fbb73820476a8b644e5e3b36a29 Mon Sep 17 00:00:00 2001 From: Piet Jarmatz Date: Thu, 13 Aug 2026 14:55:00 +0200 Subject: [PATCH 11/16] Temporary hack: _boundaryTreatment disabled for better performance (TODO fix) --- simplemd/MolecularDynamicsSimulation.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/simplemd/MolecularDynamicsSimulation.cpp b/simplemd/MolecularDynamicsSimulation.cpp index 96417e7d3..60982ee7c 100644 --- a/simplemd/MolecularDynamicsSimulation.cpp +++ b/simplemd/MolecularDynamicsSimulation.cpp @@ -531,7 +531,8 @@ void simplemd::MolecularDynamicsSimulation::simulateOneTimestep(const unsigned i // buffers to arrive, then compute forces near boundaries, all within // putBoundaryParticlesToInnerCellsFillBoundaryCellsAndOverlapWithForceComputations. if (!_configuration.getSimulationConfiguration().useOverlappingCommunicationWithForceComputation()) { - _boundaryTreatment->putBoundaryParticlesToInnerCellsAndFillBoundaryCells(_localBoundary, *_parallelTopologyService); + // TODO refactor boundaryTreatment + //_boundaryTreatment->putBoundaryParticlesToInnerCellsAndFillBoundaryCells(_localBoundary, *_parallelTopologyService); // compute forces between molecules. #if defined(KOKKOS_TARGET_CUDA) _moleculeService->getContainer().iterateMoleculesWithCell(*_lennardJonesForce); @@ -548,7 +549,8 @@ void simplemd::MolecularDynamicsSimulation::simulateOneTimestep(const unsigned i evaluateStatistics(t); - _boundaryTreatment->emptyGhostBoundaryCells(); + // TODO refactor boundaryTreatment + //_boundaryTreatment->emptyGhostBoundaryCells(); // plot VTK output if ((_configuration.getVTKConfiguration().getWriteEveryTimestep() > 0) && (t % _configuration.getVTKConfiguration().getWriteEveryTimestep() == 0)) { From 048ac2e3ae11272d8487d0eb4fb2b7aea7b8da01 Mon Sep 17 00:00:00 2001 From: Piet Jarmatz Date: Fri, 14 Aug 2026 11:11:33 +0200 Subject: [PATCH 12/16] Automatic threads_per_cell selection --- simplemd/MoleculeContainer.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/simplemd/MoleculeContainer.h b/simplemd/MoleculeContainer.h index 12b6e3b06..539104189 100644 --- a/simplemd/MoleculeContainer.h +++ b/simplemd/MoleculeContainer.h @@ -462,7 +462,7 @@ template void simplemd::MoleculeContainer::iterateMoleculesWithCellSer template void simplemd::MoleculeContainer::iterateMoleculesWithCellParallel(A& a) { a.beginMoleculeIteration(); - const unsigned int threads_per_cell = 5; + const unsigned int threads_per_cell = _cellCapacity; const unsigned int length = _linkedCellNumMolecules_d.size() * threads_per_cell; Kokkos::parallel_for( "simplemd::MoleculeContainer::iterateMoleculesWithCellParallel", Kokkos::RangePolicy(0, length), From fdedf743a975552e1707baf6151221bd9e83f853 Mon Sep 17 00:00:00 2001 From: Piet Jarmatz Date: Fri, 14 Aug 2026 14:24:43 +0200 Subject: [PATCH 13/16] Added moleculewithcellmappings::LennardJonesForceMapping::handleMoleculeVeryFast --- simplemd/MoleculeContainer.h | 29 ++++++++----- .../cell-mappings/LennardJonesForceMapping.h | 4 +- .../LennardJonesForceMapping.cpp | 41 ++++++++++++++++++- .../LennardJonesForceMapping.h | 8 +++- 4 files changed, 66 insertions(+), 16 deletions(-) diff --git a/simplemd/MoleculeContainer.h b/simplemd/MoleculeContainer.h index 539104189..baa1c6c97 100644 --- a/simplemd/MoleculeContainer.h +++ b/simplemd/MoleculeContainer.h @@ -235,7 +235,7 @@ class simplemd::MoleculeContainer { /** * @brief applies molecule-with-cell mapping to all neighbors of cell */ - template KOKKOS_FUNCTION void handleCellNeighbors(A& a, Molecule& m, const LinkedCell& cell) const; + template KOKKOS_FUNCTION void handleCellNeighbors(A& a, Molecule& m, unsigned int index) const; /** * @brief applies molecule-with-cell mapping without any node-level parallelisation @@ -439,11 +439,9 @@ template void simplemd::MoleculeContainer::iterateMoleculesWithCell(A& } } -template void simplemd::MoleculeContainer::handleCellNeighbors(A& a, Molecule& m, const LinkedCell& cell) const { - unsigned int index = cell.getIndex(); +template void simplemd::MoleculeContainer::handleCellNeighbors(A& a, Molecule& m, unsigned int index) const { for (unsigned int i = 0; i < 26; i++) { - auto cell2 = (*this)[index + _neighborOffsets(i)]; - a.handleMolecule(m, cell2); + a.handleMoleculeVeryFast(m, index + _neighborOffsets(i)); } } @@ -461,19 +459,30 @@ template void simplemd::MoleculeContainer::iterateMoleculesWithCellSer } template void simplemd::MoleculeContainer::iterateMoleculesWithCellParallel(A& a) { - a.beginMoleculeIteration(); + Kokkos::View posData("posData",_linkedCellNumMolecules_d.size(),_cellCapacity); + Kokkos::parallel_for("simplemd::MoleculeContainer fill posData", Kokkos::RangePolicy(0, _linkedCellNumMolecules_d.size() * _cellCapacity), + KOKKOS_CLASS_LAMBDA(const unsigned int i) { + const unsigned int cellIndex = i / _cellCapacity; + unsigned int molIndex = i % _cellCapacity; + if (molIndex < _linkedCellNumMolecules_d(cellIndex)) { + const Molecule& m = getMoleculeAt(cellIndex, molIndex); + posData(cellIndex,molIndex,0) = m.getConstPosition()[0]; + posData(cellIndex,molIndex,1) = m.getConstPosition()[1]; + posData(cellIndex,molIndex,2) = m.getConstPosition()[2]; + } + }); + a.beginMoleculeIteration(posData, _linkedCellNumMolecules_d); + const unsigned int threads_per_cell = _cellCapacity; const unsigned int length = _linkedCellNumMolecules_d.size() * threads_per_cell; Kokkos::parallel_for( "simplemd::MoleculeContainer::iterateMoleculesWithCellParallel", Kokkos::RangePolicy(0, length), KOKKOS_CLASS_LAMBDA(const unsigned int i) { const unsigned int cellIndex = i / threads_per_cell; - simplemd::LinkedCell cell = (*this)[cellIndex]; - for (unsigned int j = i % threads_per_cell; j < _linkedCellNumMolecules_d(cellIndex); j+=threads_per_cell) { Molecule& m = getMoleculeAt(cellIndex, j); - a.handleMolecule(m,cell); - handleCellNeighbors(a, m, cell); + a.handleMoleculeVeryFast(m,cellIndex); + handleCellNeighbors(a, m, cellIndex); } }); diff --git a/simplemd/cell-mappings/LennardJonesForceMapping.h b/simplemd/cell-mappings/LennardJonesForceMapping.h index aab940cb4..a5bb28719 100644 --- a/simplemd/cell-mappings/LennardJonesForceMapping.h +++ b/simplemd/cell-mappings/LennardJonesForceMapping.h @@ -42,15 +42,13 @@ class simplemd::cellmappings::LennardJonesForceMapping { static const bool IsParallel = true; static const bool IsReadonly = false; -private: +protected: /** epsilon */ const double _epsilon; /** sigma^6 */ const double _sigma6; /** cutOffRadius*cutOffRadius */ const double _cutOffRadiusSquared; - -protected: /** external forces*/ tarch::la::Vector _externalForce; }; diff --git a/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.cpp b/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.cpp index 0094e2966..d83c79ee5 100644 --- a/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.cpp +++ b/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.cpp @@ -1,13 +1,17 @@ #include "simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.h" +#include "tarch/utils/Utils.h" simplemd::moleculewithcellmappings::LennardJonesForceMapping::LennardJonesForceMapping( simplemd::services::ExternalForceService& externalForceService, const simplemd::services::MolecularPropertiesService& molecularPropertiesService) : simplemd::cellmappings::LennardJonesForceMapping(externalForceService, molecularPropertiesService) {} -void simplemd::moleculewithcellmappings::LennardJonesForceMapping::beginMoleculeIteration() { +void simplemd::moleculewithcellmappings::LennardJonesForceMapping::beginMoleculeIteration(const Kokkos::View& posData, + const Kokkos::View& linkedCellNumMolecules) { #if (MD_DEBUG == MD_YES) Kokkos::printf("simplemd::moleculewithcellmappings::LennardJonesForceMapping::beginMoleculeIteration()\n"); #endif + _posData = posData; + _linkedCellNumMolecules = linkedCellNumMolecules; } /* @@ -62,4 +66,37 @@ void simplemd::moleculewithcellmappings::LennardJonesForceMapping::handleMolecul #endif addForce(force1, forceBuffer); } -} \ No newline at end of file +} + +void simplemd::moleculewithcellmappings::LennardJonesForceMapping::handleMoleculeVeryFast(Molecule& m, unsigned int cellIndex) const { + tarch::la::Vector& target = m.getForce(); +#if (TARCH_DEBUG == TARCH_YES) + if (_externalForce != tarch::la::Vector{0.0}) { + Kokkos::abort( + "ERROR simplemd::moleculewithcellmappings::LennardJonesForceMapping::handleCell(): externalForce not implemented in fixed point math debug mode!\n"); + } +#else + target += _externalForce; +#endif + for(unsigned int molIndex=0; molIndex<_linkedCellNumMolecules(cellIndex);molIndex++){ + const double rijx = _posData(cellIndex,molIndex,0) - m.getConstPosition()[0]; + const double rijy = _posData(cellIndex,molIndex,1) - m.getConstPosition()[1]; + const double rijz = _posData(cellIndex,molIndex,2) - m.getConstPosition()[2]; + const double rij2 = rijx*rijx + rijy*rijy + rijz*rijz; + if (rij2 <= _cutOffRadiusSquared && rij2 > 0) { + const double rij6 = rij2 * rij2 * rij2; + const double val = 24.0 * _epsilon / rij2 * (_sigma6 / rij6) + * (1.0 - 2.0 * (_sigma6 / rij6)); +#if (TARCH_DEBUG == TARCH_YES) + DEFINE_DECIMAL_FP_LIMITS(6); + *(long long*)(&target[0]) += (long long)(stepFP6 * (val * rijx)); + *(long long*)(&target[1]) += (long long)(stepFP6 * (val * rijy)); + *(long long*)(&target[2]) += (long long)(stepFP6 * (val * rijz)); +#else + target[0] += val * rijx; + target[1] += val * rijy; + target[2] += val * rijz; +#endif + } + } +} diff --git a/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.h b/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.h index 1ae3db53a..3483d707b 100644 --- a/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.h +++ b/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.h @@ -16,11 +16,17 @@ class simplemd::moleculewithcellmappings::LennardJonesForceMapping : public simp const simplemd::services::MolecularPropertiesService& molecularPropertiesService); KOKKOS_FUNCTION virtual ~LennardJonesForceMapping() {} - void beginMoleculeIteration(); + void beginMoleculeIteration(const Kokkos::View& posData, + const Kokkos::View& linkedCellNumMolecules); void endMoleculeIteration() { } KOKKOS_FUNCTION void handleMolecule(Molecule& molecule, const LinkedCell& cell) const; + KOKKOS_FUNCTION void handleMoleculeVeryFast(Molecule& molecule, unsigned int cellIndex) const; static const bool IsParallel = true; static const bool IsReadonly = false; + +private: + Kokkos::View _posData; + Kokkos::View _linkedCellNumMolecules; }; From 1bca1e8524415f3984720edc7f92666f66d5ccef Mon Sep 17 00:00:00 2001 From: Piet Jarmatz Date: Fri, 14 Aug 2026 15:27:02 +0200 Subject: [PATCH 14/16] LJ: size_T replaced by int (better performance) --- simplemd/MoleculeContainer.h | 32 +++++++++++-------- .../LennardJonesForceMapping.cpp | 8 +++-- .../LennardJonesForceMapping.h | 6 ++-- 3 files changed, 26 insertions(+), 20 deletions(-) diff --git a/simplemd/MoleculeContainer.h b/simplemd/MoleculeContainer.h index baa1c6c97..ebc99cc29 100644 --- a/simplemd/MoleculeContainer.h +++ b/simplemd/MoleculeContainer.h @@ -235,7 +235,7 @@ class simplemd::MoleculeContainer { /** * @brief applies molecule-with-cell mapping to all neighbors of cell */ - template KOKKOS_FUNCTION void handleCellNeighbors(A& a, Molecule& m, unsigned int index) const; + template KOKKOS_FUNCTION void handleCellNeighbors(A& a, Molecule& m, int index) const; /** * @brief applies molecule-with-cell mapping without any node-level parallelisation @@ -439,8 +439,8 @@ template void simplemd::MoleculeContainer::iterateMoleculesWithCell(A& } } -template void simplemd::MoleculeContainer::handleCellNeighbors(A& a, Molecule& m, unsigned int index) const { - for (unsigned int i = 0; i < 26; i++) { +template void simplemd::MoleculeContainer::handleCellNeighbors(A& a, Molecule& m, int index) const { + for (int i = 0; i < 26; i++) { a.handleMoleculeVeryFast(m, index + _neighborOffsets(i)); } } @@ -459,33 +459,37 @@ template void simplemd::MoleculeContainer::iterateMoleculesWithCellSer } template void simplemd::MoleculeContainer::iterateMoleculesWithCellParallel(A& a) { + Kokkos::View linkedCellNumMolecules("linkedCellNumMolecules_int",_linkedCellNumMolecules_d.size()); + Kokkos::parallel_for("simplemd::MoleculeContainer fill linkedCellNumMolecules", Kokkos::RangePolicy(0, _linkedCellNumMolecules_d.size()), + KOKKOS_CLASS_LAMBDA(const int i) { + linkedCellNumMolecules(i) = (int)_linkedCellNumMolecules_d(i); + }); Kokkos::View posData("posData",_linkedCellNumMolecules_d.size(),_cellCapacity); Kokkos::parallel_for("simplemd::MoleculeContainer fill posData", Kokkos::RangePolicy(0, _linkedCellNumMolecules_d.size() * _cellCapacity), - KOKKOS_CLASS_LAMBDA(const unsigned int i) { - const unsigned int cellIndex = i / _cellCapacity; - unsigned int molIndex = i % _cellCapacity; - if (molIndex < _linkedCellNumMolecules_d(cellIndex)) { + KOKKOS_CLASS_LAMBDA(const int i) { + const int cellIndex = i / _cellCapacity; + int molIndex = i % _cellCapacity; + if (molIndex < (int)_linkedCellNumMolecules_d(cellIndex)) { const Molecule& m = getMoleculeAt(cellIndex, molIndex); posData(cellIndex,molIndex,0) = m.getConstPosition()[0]; posData(cellIndex,molIndex,1) = m.getConstPosition()[1]; posData(cellIndex,molIndex,2) = m.getConstPosition()[2]; } }); - a.beginMoleculeIteration(posData, _linkedCellNumMolecules_d); + a.beginMoleculeIteration(posData, linkedCellNumMolecules); - const unsigned int threads_per_cell = _cellCapacity; - const unsigned int length = _linkedCellNumMolecules_d.size() * threads_per_cell; + const int threads_per_cell = _cellCapacity; + const int length = _linkedCellNumMolecules_d.size() * threads_per_cell; Kokkos::parallel_for( "simplemd::MoleculeContainer::iterateMoleculesWithCellParallel", Kokkos::RangePolicy(0, length), - KOKKOS_CLASS_LAMBDA(const unsigned int i) { - const unsigned int cellIndex = i / threads_per_cell; - for (unsigned int j = i % threads_per_cell; j < _linkedCellNumMolecules_d(cellIndex); j+=threads_per_cell) { + KOKKOS_CLASS_LAMBDA(const int i) { + const int cellIndex = i / threads_per_cell; + for (int j = i % threads_per_cell; j < linkedCellNumMolecules(cellIndex); j+=threads_per_cell) { Molecule& m = getMoleculeAt(cellIndex, j); a.handleMoleculeVeryFast(m,cellIndex); handleCellNeighbors(a, m, cellIndex); } }); - a.endMoleculeIteration(); } diff --git a/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.cpp b/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.cpp index d83c79ee5..cebcabae6 100644 --- a/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.cpp +++ b/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.cpp @@ -6,7 +6,7 @@ simplemd::moleculewithcellmappings::LennardJonesForceMapping::LennardJonesForceM : simplemd::cellmappings::LennardJonesForceMapping(externalForceService, molecularPropertiesService) {} void simplemd::moleculewithcellmappings::LennardJonesForceMapping::beginMoleculeIteration(const Kokkos::View& posData, - const Kokkos::View& linkedCellNumMolecules) { + const Kokkos::View& linkedCellNumMolecules) { #if (MD_DEBUG == MD_YES) Kokkos::printf("simplemd::moleculewithcellmappings::LennardJonesForceMapping::beginMoleculeIteration()\n"); #endif @@ -68,17 +68,19 @@ void simplemd::moleculewithcellmappings::LennardJonesForceMapping::handleMolecul } } -void simplemd::moleculewithcellmappings::LennardJonesForceMapping::handleMoleculeVeryFast(Molecule& m, unsigned int cellIndex) const { +void simplemd::moleculewithcellmappings::LennardJonesForceMapping::handleMoleculeVeryFast(Molecule& m, int cellIndex) const { tarch::la::Vector& target = m.getForce(); #if (TARCH_DEBUG == TARCH_YES) +#if (MD_DEBUG == MD_YES) if (_externalForce != tarch::la::Vector{0.0}) { Kokkos::abort( "ERROR simplemd::moleculewithcellmappings::LennardJonesForceMapping::handleCell(): externalForce not implemented in fixed point math debug mode!\n"); } +#endif #else target += _externalForce; #endif - for(unsigned int molIndex=0; molIndex<_linkedCellNumMolecules(cellIndex);molIndex++){ + for(int molIndex=0; molIndex<_linkedCellNumMolecules(cellIndex);molIndex++){ const double rijx = _posData(cellIndex,molIndex,0) - m.getConstPosition()[0]; const double rijy = _posData(cellIndex,molIndex,1) - m.getConstPosition()[1]; const double rijz = _posData(cellIndex,molIndex,2) - m.getConstPosition()[2]; diff --git a/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.h b/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.h index 3483d707b..fcf4540f8 100644 --- a/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.h +++ b/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.h @@ -17,16 +17,16 @@ class simplemd::moleculewithcellmappings::LennardJonesForceMapping : public simp KOKKOS_FUNCTION virtual ~LennardJonesForceMapping() {} void beginMoleculeIteration(const Kokkos::View& posData, - const Kokkos::View& linkedCellNumMolecules); + const Kokkos::View& linkedCellNumMolecules); void endMoleculeIteration() { } KOKKOS_FUNCTION void handleMolecule(Molecule& molecule, const LinkedCell& cell) const; - KOKKOS_FUNCTION void handleMoleculeVeryFast(Molecule& molecule, unsigned int cellIndex) const; + KOKKOS_FUNCTION void handleMoleculeVeryFast(Molecule& molecule, int cellIndex) const; static const bool IsParallel = true; static const bool IsReadonly = false; private: Kokkos::View _posData; - Kokkos::View _linkedCellNumMolecules; + Kokkos::View _linkedCellNumMolecules; }; From 6f87b6a863611d7236b732ccebd662e751275b0f Mon Sep 17 00:00:00 2001 From: Piet Jarmatz Date: Fri, 14 Aug 2026 15:35:31 +0200 Subject: [PATCH 15/16] Add ES1 optimization to LJForce --- simplemd/MoleculeContainer.h | 40 ++++++++++++++----- .../cell-mappings/LennardJonesForceMapping.h | 2 + .../LennardJonesForceMapping.cpp | 35 ++++++++++++++++ .../LennardJonesForceMapping.h | 1 + 4 files changed, 68 insertions(+), 10 deletions(-) diff --git a/simplemd/MoleculeContainer.h b/simplemd/MoleculeContainer.h index ebc99cc29..2d9ef1cfa 100644 --- a/simplemd/MoleculeContainer.h +++ b/simplemd/MoleculeContainer.h @@ -236,6 +236,7 @@ class simplemd::MoleculeContainer { * @brief applies molecule-with-cell mapping to all neighbors of cell */ template KOKKOS_FUNCTION void handleCellNeighbors(A& a, Molecule& m, int index) const; + template KOKKOS_FUNCTION void handleCellNeighborsES1(A& a, Molecule& m, int index) const; /** * @brief applies molecule-with-cell mapping without any node-level parallelisation @@ -445,6 +446,12 @@ template void simplemd::MoleculeContainer::handleCellNeighbors(A& a, M } } +template void simplemd::MoleculeContainer::handleCellNeighborsES1(A& a, Molecule& m, int index) const { + for (int i = 0; i < 26; i++) { + a.handleMoleculeVeryFastES1(m, index + _neighborOffsets(i)); + } +} + template void simplemd::MoleculeContainer::iterateMoleculesWithCellSerial(A& a) { a.beginMoleculeIteration(); for (unsigned int i = 0; i < _linkedCellNumMolecules_h.size(); i++) { @@ -480,16 +487,29 @@ template void simplemd::MoleculeContainer::iterateMoleculesWithCellPar const int threads_per_cell = _cellCapacity; const int length = _linkedCellNumMolecules_d.size() * threads_per_cell; - Kokkos::parallel_for( - "simplemd::MoleculeContainer::iterateMoleculesWithCellParallel", Kokkos::RangePolicy(0, length), - KOKKOS_CLASS_LAMBDA(const int i) { - const int cellIndex = i / threads_per_cell; - for (int j = i % threads_per_cell; j < linkedCellNumMolecules(cellIndex); j+=threads_per_cell) { - Molecule& m = getMoleculeAt(cellIndex, j); - a.handleMoleculeVeryFast(m,cellIndex); - handleCellNeighbors(a, m, cellIndex); - } - }); + if(a.epsilon_sigma_is_one()){ + Kokkos::parallel_for( + "simplemd::MoleculeContainer::iterateMoleculesWithCellParallel", Kokkos::RangePolicy(0, length), + KOKKOS_CLASS_LAMBDA(const int i) { + const int cellIndex = i / threads_per_cell; + for (int j = i % threads_per_cell; j < linkedCellNumMolecules(cellIndex); j+=threads_per_cell) { + Molecule& m = getMoleculeAt(cellIndex, j); + a.handleMoleculeVeryFastES1(m,cellIndex); + handleCellNeighborsES1(a, m, cellIndex); + } + }); + } else { + Kokkos::parallel_for( + "simplemd::MoleculeContainer::iterateMoleculesWithCellParallel", Kokkos::RangePolicy(0, length), + KOKKOS_CLASS_LAMBDA(const int i) { + const int cellIndex = i / threads_per_cell; + for (int j = i % threads_per_cell; j < linkedCellNumMolecules(cellIndex); j+=threads_per_cell) { + Molecule& m = getMoleculeAt(cellIndex, j); + a.handleMoleculeVeryFast(m,cellIndex); + handleCellNeighbors(a, m, cellIndex); + } + }); + } a.endMoleculeIteration(); } diff --git a/simplemd/cell-mappings/LennardJonesForceMapping.h b/simplemd/cell-mappings/LennardJonesForceMapping.h index a5bb28719..9df22d33d 100644 --- a/simplemd/cell-mappings/LennardJonesForceMapping.h +++ b/simplemd/cell-mappings/LennardJonesForceMapping.h @@ -42,6 +42,8 @@ class simplemd::cellmappings::LennardJonesForceMapping { static const bool IsParallel = true; static const bool IsReadonly = false; + bool epsilon_sigma_is_one() const {return _epsilon == 1.0 && _sigma6 == 1.0;}; + protected: /** epsilon */ const double _epsilon; diff --git a/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.cpp b/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.cpp index cebcabae6..d78595595 100644 --- a/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.cpp +++ b/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.cpp @@ -102,3 +102,38 @@ void simplemd::moleculewithcellmappings::LennardJonesForceMapping::handleMolecul } } } + +void simplemd::moleculewithcellmappings::LennardJonesForceMapping::handleMoleculeVeryFastES1(Molecule& m, int cellIndex) const { + tarch::la::Vector& target = m.getForce(); +#if (TARCH_DEBUG == TARCH_YES) +#if (MD_DEBUG == MD_YES) + if (_externalForce != tarch::la::Vector{0.0}) { + Kokkos::abort( + "ERROR simplemd::moleculewithcellmappings::LennardJonesForceMapping::handleCell(): externalForce not implemented in fixed point math debug mode!\n"); + } +#endif +#else + target += _externalForce; +#endif + for(int molIndex=0; molIndex<_linkedCellNumMolecules(cellIndex);molIndex++){ + const double rijx = _posData(cellIndex,molIndex,0) - m.getConstPosition()[0]; + const double rijy = _posData(cellIndex,molIndex,1) - m.getConstPosition()[1]; + const double rijz = _posData(cellIndex,molIndex,2) - m.getConstPosition()[2]; + const double rij2 = rijx*rijx + rijy*rijy + rijz*rijz; + if (rij2 <= _cutOffRadiusSquared && rij2 > 0) { + const double rij6 = rij2 * rij2 * rij2; + const double val = 24.0 / rij2 * (1.0 / rij6) + * (1.0 - 2.0 * (1.0 / rij6)); +#if (TARCH_DEBUG == TARCH_YES) + DEFINE_DECIMAL_FP_LIMITS(6); + *(long long*)(&target[0]) += (long long)(stepFP6 * (val * rijx)); + *(long long*)(&target[1]) += (long long)(stepFP6 * (val * rijy)); + *(long long*)(&target[2]) += (long long)(stepFP6 * (val * rijz)); +#else + target[0] += val * rijx; + target[1] += val * rijy; + target[2] += val * rijz; +#endif + } + } +} diff --git a/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.h b/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.h index fcf4540f8..a67ce4aa8 100644 --- a/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.h +++ b/simplemd/molecule-with-cell-mappings/LennardJonesForceMapping.h @@ -22,6 +22,7 @@ class simplemd::moleculewithcellmappings::LennardJonesForceMapping : public simp KOKKOS_FUNCTION void handleMolecule(Molecule& molecule, const LinkedCell& cell) const; KOKKOS_FUNCTION void handleMoleculeVeryFast(Molecule& molecule, int cellIndex) const; + KOKKOS_FUNCTION void handleMoleculeVeryFastES1(Molecule& molecule, int cellIndex) const; static const bool IsParallel = true; static const bool IsReadonly = false; From 9c4e1a3529d480fb2156205306bf55113d201519 Mon Sep 17 00:00:00 2001 From: Piet Jarmatz Date: Fri, 14 Aug 2026 17:10:17 +0200 Subject: [PATCH 16/16] Added SimpleMDBenchSize::MD480 --- test/integration/SimpleMDBench.h | 9 ++++++++- test/integration/benchmarks.cpp | 1 + 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/test/integration/SimpleMDBench.h b/test/integration/SimpleMDBench.h index 8477bbfcb..606b9d1cc 100644 --- a/test/integration/SimpleMDBench.h +++ b/test/integration/SimpleMDBench.h @@ -11,7 +11,7 @@ #include #include -enum class SimpleMDBenchSize { MD60, MD120, MD240 }; +enum class SimpleMDBenchSize { MD60, MD120, MD240, MD480 }; class BenchSim : public simplemd::MolecularDynamicsSimulation { public: @@ -107,6 +107,7 @@ class SimpleMDBench : public Test { linked-cell-size="2.5 ; 2.5 ; 2.5" k_B="1.0" block-size="100" + capacity-factor="2.25" bottom-south-west="reflecting" bottom-south="reflecting" bottom-south-east="reflecting" bottom-west="reflecting" bottom="reflecting" bottom-east="reflecting" @@ -139,6 +140,11 @@ class SimpleMDBench : public Test { = tarch::la::Vector<3, unsigned int>{224,224,224}; _simpleMDConfig.getDomainConfigurationNonConst().getGlobalDomainSizeNonConst() = tarch::la::Vector<3, double>{240,240,240}; + } else if(_MDSize == SimpleMDBenchSize::MD480){ + _simpleMDConfig.getDomainConfigurationNonConst().getMoleculesPerDirectionNonConst() + = tarch::la::Vector<3, unsigned int>{448,448,448}; + _simpleMDConfig.getDomainConfigurationNonConst().getGlobalDomainSizeNonConst() + = tarch::la::Vector<3, double>{480,480,480}; } _simulation = std::make_unique(_simpleMDConfig); @@ -193,6 +199,7 @@ class SimpleMDBench : public Test { if(_MDSize == SimpleMDBenchSize::MD60) correct = 9224833479527670225u; if(_MDSize == SimpleMDBenchSize::MD120) correct = 18377339084083654394u; if(_MDSize == SimpleMDBenchSize::MD240) correct = 18342705775936782749u; + if(_MDSize == SimpleMDBenchSize::MD480) correct = 18349058309621237939u; if (sum == correct) std::cout << "INFO SimpleMDBench: SUCCESS Checksum is correct :-)" << std::endl; else { diff --git a/test/integration/benchmarks.cpp b/test/integration/benchmarks.cpp index 7eca9d99d..fb3e3c976 100644 --- a/test/integration/benchmarks.cpp +++ b/test/integration/benchmarks.cpp @@ -39,6 +39,7 @@ int main(int argc, char* argv[]) { } else if(argc == 2 && !strcmp(argv[1],"MD60")) benchsize = SimpleMDBenchSize::MD60; else if(argc == 2 && !strcmp(argv[1],"MD120")) benchsize = SimpleMDBenchSize::MD120; else if(argc == 2 && !strcmp(argv[1],"MD240")) benchsize = SimpleMDBenchSize::MD240; + else if(argc == 2 && !strcmp(argv[1],"MD480")) benchsize = SimpleMDBenchSize::MD480; else { std::cout << "ERROR unknown parameter" << std::endl; Kokkos::finalize();