diff --git a/.gitignore b/.gitignore index f198a0bfc0..3a9cb7f5ac 100644 --- a/.gitignore +++ b/.gitignore @@ -6,6 +6,7 @@ build*/ .idea/ .project .settings +.vscode # Generated documentation doxygen_doc/ diff --git a/cmake/modules/mamico.cmake b/cmake/modules/mamico.cmake index 5d4a6ed549..6939b2c601 100644 --- a/cmake/modules/mamico.cmake +++ b/cmake/modules/mamico.cmake @@ -1,18 +1,15 @@ option(MAMICO_COUPLING "Couple with MaMiCo" OFF) if (MAMICO_COUPLING) - message(STATUS "MaMiCo coupling enabled. ls1 mardyn will compile as library. No executable will be created.") - set(MAMICO_COMPILE_DEFINITIONS MAMICO_COUPLING MDDim3) - option(MAMICO_ENABLE_FPIC "Enable -fPIC flag for MaMiCo python bindings" OFF) - set(MAMICO_SRC_DIR CACHE PATH "Root directory of the MaMiCo codebase") + message(STATUS "MaMiCo coupling enabled. ls1 mardyn will compile as library. No executable will be created.") + set(MAMICO_COMPILE_DEFINITIONS MAMICO_COUPLING MDDim3) + option(MAMICO_ENABLE_FPIC "Enable -fPIC flag for MaMiCo python bindings" OFF) + set(MAMICO_SRC_DIR CACHE PATH "Root directory of the MaMiCo codebase") if(NOT MAMICO_SRC_DIR) - message(FATAL_ERROR "MaMiCo source directory not specified.") - endif() - if(ENABLE_MPI) - set(MAMICO_MPI_DEFINITIONS MDCoupledParallel TarchParallel) - endif() - if(MAMICO_ENABLE_FPIC) - set(MAMICO_COMPILE_OPTIONS "${MAMICO_COMPILE_OPTIONS} -fPIC") - endif() + message(FATAL_ERROR "MaMiCo source directory not specified.") + endif() + if(ENABLE_MPI) + set(MAMICO_MPI_DEFINITIONS MDCoupledParallel TarchParallel) + endif() else() message(STATUS "MaMiCo coupling disabled.") endif() \ No newline at end of file diff --git a/examples/all-options.xml b/examples/all-options.xml index 389f25f141..855bad744b 100644 --- a/examples/all-options.xml +++ b/examples/all-options.xml @@ -184,7 +184,8 @@ 0.0 0.0 0.0 - + + true @@ -196,6 +197,16 @@ False 5 False + + + reflective + outflow + periodic + + + reflectiveoutflowreflective + + + + + 2.2 + + + 1.0e+10 + + + + + vtkOutput + 1 + + + particles.bp + BP4 + 1 + + + + + diff --git a/examples/simple-boundary-test/simple_checkpoint.inp b/examples/simple-boundary-test/simple_checkpoint.inp new file mode 100644 index 0000000000..1da67c4a9c --- /dev/null +++ b/examples/simple-boundary-test/simple_checkpoint.inp @@ -0,0 +1,21 @@ +mardyn trunk 20120726 +currentTime 0 +Length 10 10 10 +Temperature 1.1 +NumberOfComponents 1 +1 0 0 0 0 +0 0 0 1 1 1 0 +0 0 0 +1e+10 +NumberOfMolecules 10 +MoleculeFormat IRV +0 7.5 5 5 -3 0 0 +1 2.5 5 5 0 0 0 +2 2.5 7.5 5 0 0 0 +3 2.5 2.5 5 0 0 0 +4 2.5 5 2.5 0 0 0 +5 2.5 5 7.5 0 0 0 +6 2.5 7.5 7.5 0 0 0 +7 2.5 7.5 2.5 0 0 0 +8 2.5 2.5 7.5 0 0 0.69 +9 2.5 2.5 2.5 0 -1.1 -1 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index f204718945..d2a8cd8894 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -8,16 +8,16 @@ if(NOT ENABLE_UNIT_TESTS) list(FILTER MY_SRC EXCLUDE REGEX "/tests/") endif() -# if mpi is not enabled, remove the uneeded source files +# if mpi is not enabled, remove the unneeded source files if(NOT ENABLE_MPI) # duplicate the list set(MY_SRC_BACK ${MY_SRC}) # exclude everything from parallel list(FILTER MY_SRC EXCLUDE REGEX "/parallel/") - # but include DomainDecompBase* and LoadCalc* + # but include DomainDecompBase*, LoadCalc* and Boundary utilities list(FILTER MY_SRC_BACK INCLUDE REGEX "/parallel/") - list(FILTER MY_SRC_BACK INCLUDE REGEX "DomainDecompBase|LoadCalc|Zonal|ForceHelper") + list(FILTER MY_SRC_BACK INCLUDE REGEX "boundaries/|DomainDecompBase|LoadCalc|Zonal|ForceHelper") list(APPEND MY_SRC ${MY_SRC_BACK}) else() if(NOT ENABLE_ALLLBL) @@ -54,7 +54,9 @@ if (MAMICO_COUPLING) TARGET_COMPILE_DEFINITIONS(MarDyn PUBLIC ${MAMICO_COMPILE_DEFINITIONS} ${MAMICO_MPI_DEFINITIONS} ) - TARGET_COMPILE_OPTIONS(MarDyn PUBLIC ${MAMICO_COMPILE_OPTIONS}) + if(MAMICO_ENABLE_FPIC) + SET_PROPERTY(TARGET MarDyn PROPERTY POSITION_INDEPENDENT_CODE ON) + endif() else() ADD_EXECUTABLE(MarDyn ${MY_SRC} diff --git a/src/Domain.cpp b/src/Domain.cpp index 00a3a3bf8b..80657e2f36 100644 --- a/src/Domain.cpp +++ b/src/Domain.cpp @@ -145,9 +145,9 @@ double Domain::getGlobalPressure() return globalTemperature * _globalRho + _globalRho * getAverageGlobalVirial()/3.; } -double Domain::getAverageGlobalVirial() { return _globalVirial/_globalNumMolecules; } +double Domain::getAverageGlobalVirial() const { return _globalVirial/_globalNumMolecules; } -double Domain::getAverageGlobalUpot() { return getGlobalUpot()/_globalNumMolecules; } +double Domain::getAverageGlobalUpot() const { return getGlobalUpot()/_globalNumMolecules; } double Domain::getGlobalUpot() const { return _globalUpot; } Comp2Param& Domain::getComp2Params(){ @@ -791,9 +791,9 @@ void Domain::updateMaxMoleculeID(ParticleContainer* particleContainer, DomainDec #endif } -double Domain::getglobalRho(){ return _globalRho;} +double Domain::getglobalRho() const { return _globalRho;} -void Domain::setglobalRho(double grho){ _globalRho = grho;} +void Domain::setglobalRho(double grho) { _globalRho = grho;} unsigned long Domain::getglobalRotDOF() { @@ -832,10 +832,10 @@ double Domain::cv() //! methods implemented by Stefan Becker // the following two methods are used by the MmspdWriter (writing the output file in a format used by MegaMol) -double Domain::getSigma(unsigned cid, unsigned nthSigma){ +double Domain::getSigma(unsigned cid, unsigned nthSigma) const { return _simulation.getEnsemble()->getComponent(cid)->getSigma(nthSigma); } -unsigned Domain::getNumberOfComponents(){ +unsigned Domain::getNumberOfComponents() const { return _simulation.getEnsemble()->getComponents()->size(); } @@ -863,10 +863,10 @@ void Domain::submitDU(unsigned /*cid*/, double DU, double* r) void Domain::setLocalUpotCompSpecific(double UpotCspec){_localUpotCspecif = UpotCspec;} - double Domain::getLocalUpotCompSpecific(){return _localUpotCspecif;} + double Domain::getLocalUpotCompSpecific() const {return _localUpotCspecif;} -double Domain::getAverageGlobalUpotCSpec() { +double Domain::getAverageGlobalUpotCSpec() const { Log::global_log->debug() << "number of fluid molecules = " << getNumFluidMolecules() << "\n"; return _globalUpotCspecif / getNumFluidMolecules(); } @@ -876,7 +876,7 @@ void Domain::setNumFluidComponents(unsigned nc){_numFluidComponent = nc;} unsigned Domain::getNumFluidComponents(){return _numFluidComponent;} -unsigned long Domain::getNumFluidMolecules(){ +unsigned long Domain::getNumFluidMolecules() const { unsigned long numFluidMolecules = 0; for(unsigned i = 0; i < _numFluidComponent; i++){ Component& ci=*(global_simulation->getEnsemble()->getComponent(i)); diff --git a/src/Domain.h b/src/Domain.h index faabdca061..cabc179377 100644 --- a/src/Domain.h +++ b/src/Domain.h @@ -140,7 +140,7 @@ class Domain { unsigned getNumFluidComponents(); //! @brief get the fluid and fluid-solid potential of the local process - double getLocalUpotCompSpecific(); + double getLocalUpotCompSpecific() const; //! @brief set the virial of the local process void setLocalVirial(double Virial); @@ -214,28 +214,30 @@ class Domain { //! //! Before this method is called, it has to be sure that the //! global potential has been calculated (method calculateGlobalValues) - double getAverageGlobalUpot(); + double getAverageGlobalUpot() const; double getGlobalUpot() const; //! by Stefan Becker: return the average global potential of the fluid-fluid and fluid-solid interaction (but NOT solid-solid interaction) - double getAverageGlobalUpotCSpec(); + double getAverageGlobalUpotCSpec() const; //! @brief get the global kinetic energy //! - //! Before this method is called, it has to be sure that the - //! global energies has been calculated (method calculateGlobalValues) - double getGlobalUkinTrans() { return 0.5*_globalsummv2; } - double getGlobalUkinRot() { return 0.5*_globalsumIw2; } + //! Before this method is called, the user has to be sure that the + //! global energy (rot and trans) has been calculated via calculateGlobalValues() + //! Since variables _globalsummv2 and _globalsumIw2 store the sum of m_i*(v_i^2). + //! Therefore, the constant factor 0.5 has to be applied to yield the kinetic energies + double getGlobalUkinTrans() const { return 0.5*_globalsummv2; } + double getGlobalUkinRot() const { return 0.5*_globalsumIw2; } //! by Stefan Becker: determine and return the totel number of fluid molecules //! this method assumes all molecules with a component-ID less than _numFluidComponent to be fluid molecules - unsigned long getNumFluidMolecules(); + unsigned long getNumFluidMolecules() const; //! @brief get the global average virial per particle //! //! Before this method is called, it has to be sure that the //! global virial has been calculated (method calculateGlobalValues) - double getAverageGlobalVirial(); + double getAverageGlobalVirial() const; //! @brief sets _localSummv2 to the given value void setLocalSummv2(double summv2, int thermostat); @@ -250,7 +252,7 @@ class Domain { } //! @brief get globalRho - double getglobalRho(); + double getglobalRho() const; //! @brief set globalRho void setglobalRho(double grho); @@ -386,9 +388,9 @@ class Domain { // by Stefan Becker /* method returning the sigma parameter of a component => needed in the output of the MmspdWriter (specifying the particles' radii in a movie) */ - double getSigma(unsigned cid, unsigned nthSigma); + double getSigma(unsigned cid, unsigned nthSigma) const; // needed for the MmspdWriter (MegaMol) - unsigned getNumberOfComponents(); + unsigned getNumberOfComponents() const; void setUpotCorr(double upotcorr){ _UpotCorr = upotcorr; } void setVirialCorr(double virialcorr){ _VirialCorr = virialcorr; } diff --git a/src/Simulation.cpp b/src/Simulation.cpp index ac459db6f8..78b29b737b 100644 --- a/src/Simulation.cpp +++ b/src/Simulation.cpp @@ -86,6 +86,8 @@ #include #endif +#include "parallel/boundaries/BoundaryUtils.h" + Simulation* global_simulation; Simulation::Simulation() @@ -424,6 +426,36 @@ void Simulation::readXML(XMLfileUnits& xmlconfig) { } _lastTraversalTimeHistory.setCapacity(timerForLoadAveragingLength); + if(xmlconfig.changecurrentnode("boundaries")) { + std::string xBoundaryFromFile, yBoundaryFromFile, zBoundaryFromFile; + xmlconfig.getNodeValue("x", xBoundaryFromFile); + xmlconfig.getNodeValue("y", yBoundaryFromFile); + xmlconfig.getNodeValue("z", zBoundaryFromFile); + BoundaryUtils::BoundaryType xBoundary = BoundaryUtils::convertStringToBoundary(xBoundaryFromFile); + BoundaryUtils::BoundaryType yBoundary = BoundaryUtils::convertStringToBoundary(yBoundaryFromFile); + BoundaryUtils::BoundaryType zBoundary = BoundaryUtils::convertStringToBoundary(zBoundaryFromFile); + _domainDecomposition->setGlobalBoundaryType(DimensionUtils::DimensionType::POSX, xBoundary); + _domainDecomposition->setGlobalBoundaryType(DimensionUtils::DimensionType::NEGX, xBoundary); + _domainDecomposition->setGlobalBoundaryType(DimensionUtils::DimensionType::POSY, yBoundary); + _domainDecomposition->setGlobalBoundaryType(DimensionUtils::DimensionType::NEGY, yBoundary); + _domainDecomposition->setGlobalBoundaryType(DimensionUtils::DimensionType::POSZ, zBoundary); + _domainDecomposition->setGlobalBoundaryType(DimensionUtils::DimensionType::NEGZ, zBoundary); + if (_domainDecomposition->hasGlobalInvalidBoundary()) { + Log::global_log->error() << "Invalid boundary type! Please check the config file" << std::endl; + exit(1); + } + if(_overlappingP2P && _domainDecomposition->hasGlobalNonPeriodicBoundary()) { + Log::global_log->info() << "Non-periodic boundaries not supported with overlappingP2P enabled! Exiting..." << std::endl; + exit(1); + } + Log::global_log->info() << "Boundary conditions: x - " << BoundaryUtils::convertBoundaryToString(xBoundary) + << " y - " << BoundaryUtils::convertBoundaryToString(yBoundary) + << " z - " << BoundaryUtils::convertBoundaryToString(zBoundary) << std::endl; + // go over all local boundaries, determine which are global + _domainDecomposition->setLocalBoundariesFromGlobal(_domain, _ensemble); + xmlconfig.changecurrentnode(".."); + } + xmlconfig.changecurrentnode(".."); } else { @@ -1014,13 +1046,14 @@ void Simulation::preSimLoopSteps() global_simulation->timers()->setOutputString("SIMULATION_UPDATE_CACHES", "Cache update took:"); global_simulation->timers()->setOutputString("COMMUNICATION_PARTNER_INIT_SEND", "initSend() took:"); global_simulation->timers()->setOutputString("COMMUNICATION_PARTNER_TEST_RECV", "testRecv() took:"); + global_simulation->timers()->setOutputString("SIMULATION_BOUNDARY_TREATMENT", "Enforcing boundary conditions took:"); // all timers except the ioTimer measure inside the main loop //global_simulation->timers()->getTimer("SIMULATION_LOOP")->set_sync(true); //global_simulation->timers()->setSyncTimer("SIMULATION_LOOP", true); #ifdef WITH_PAPI - const char *papi_event_list[] = { "PAPI_TOT_CYC", "PAPI_TOT_INS" /*, "PAPI_VEC_DP", "PAPI_L2_DCM", "PAPI_L2_ICM", "PAPI_L1_ICM", "PAPI_DP_OPS", "PAPI_VEC_INS" }; */ + const char *papi_event_list[] = { "PAPI_TOT_CYC", "PAPI_TOT_INS" };/*, "PAPI_VEC_DP", "PAPI_L2_DCM", "PAPI_L2_ICM", "PAPI_L1_ICM", "PAPI_DP_OPS", "PAPI_VEC_INS" }; */ int num_papi_events = sizeof(papi_event_list) / sizeof(papi_event_list[0]); global_simulation->timers()->getTimer("SIMULATION_LOOP")->add_papi_counters(num_papi_events, (char**) papi_event_list); #endif @@ -1077,9 +1110,13 @@ void Simulation::simulateOneTimestep() global_simulation->timers()->stop(plugin->getPluginName()); } - _ensemble->beforeEventNewTimestep(_moleculeContainer, _domainDecomposition, _simstep); - - _integrator->eventNewTimestep(_moleculeContainer, _domain); + _ensemble->beforeEventNewTimestep(_moleculeContainer, _domainDecomposition, _simstep); + + global_simulation->timers()->start("SIMULATION_BOUNDARY_TREATMENT"); + _domainDecomposition->processBoundaryConditions(_moleculeContainer, _integrator->getTimestepLength()); + global_simulation->timers()->stop("SIMULATION_BOUNDARY_TREATMENT"); + + _integrator->eventNewTimestep(_moleculeContainer, _domain); // beforeForces Plugin Call Log::global_log -> debug() << "[BEFORE FORCES] Performing BeforeForces plugin call" << std::endl; @@ -1402,6 +1439,10 @@ void Simulation::updateParticleContainerAndDecomposition(double lastTraversalTim _domain); global_simulation->timers()->stop("SIMULATION_MPI_OMP_COMMUNICATION"); + global_simulation->timers()->start("SIMULATION_BOUNDARY_TREATMENT"); + _domainDecomposition->removeNonPeriodicHalos(_moleculeContainer); + global_simulation->timers()->stop("SIMULATION_BOUNDARY_TREATMENT"); + // The cache of the molecules must be updated/build after the exchange process, // as the cache itself isn't transferred global_simulation->timers()->start("SIMULATION_UPDATE_CACHES"); diff --git a/src/io/ResultWriter.cpp b/src/io/ResultWriter.cpp index 9212a2b021..4a8d7ff899 100644 --- a/src/io/ResultWriter.cpp +++ b/src/io/ResultWriter.cpp @@ -67,9 +67,8 @@ void ResultWriter::endStep(ParticleContainer *particleContainer, DomainDecompBas // Writing of cavities now handled by CavityWriter - unsigned long globalNumMolecules = domain->getglobalNumMolecules(true, particleContainer, domainDecomp); - double cv = domain->cv(); - double ekin = domain->getGlobalUkinTrans()+domain->getGlobalUkinRot(); + const unsigned long globalNumMolecules = domain->getglobalNumMolecules(true, particleContainer, domainDecomp); + const double ekin = domain->getGlobalUkinTrans()+domain->getGlobalUkinRot(); _U_pot_acc->addEntry(domain->getGlobalUpot()); _U_kin_acc->addEntry(ekin); @@ -91,7 +90,7 @@ void ResultWriter::endStep(ParticleContainer *particleContainer, DomainDecompBas printOutput(_p_acc->getAverage()); printOutput(domain->getGlobalBetaTrans()); printOutput(domain->getGlobalBetaRot()); - printOutput(cv); + printOutput(domain->cv()); printOutput(globalNumMolecules); resultStream << std::endl; resultStream.close(); diff --git a/src/io/TimerProfiler.cpp b/src/io/TimerProfiler.cpp index 520702e376..c8ec9d6803 100644 --- a/src/io/TimerProfiler.cpp +++ b/src/io/TimerProfiler.cpp @@ -201,6 +201,7 @@ void TimerProfiler::readInitialTimersFromFile(std::string fileName){ std::make_tuple("QUICKSCHED", std::vector{"SIMULATION_LOOP"}, true), #endif std::make_tuple("SIMULATION_PER_STEP_IO", std::vector{"SIMULATION_LOOP"}, true), + std::make_tuple("SIMULATION_BOUNDARY_TREATMENT", std::vector{"SIMULATION_LOOP"}, true), std::make_tuple("SIMULATION_IO", std::vector{"SIMULATION"}, true), std::make_tuple("SIMULATION_UPDATE_CONTAINER", std::vector{"SIMULATION_DECOMPOSITION"}, true), std::make_tuple("SIMULATION_MPI_OMP_COMMUNICATION", std::vector{"SIMULATION_DECOMPOSITION"}, true), diff --git a/src/parallel/DomainDecompBase.cpp b/src/parallel/DomainDecompBase.cpp index 6a3499e5a1..dd9b2af653 100644 --- a/src/parallel/DomainDecompBase.cpp +++ b/src/parallel/DomainDecompBase.cpp @@ -15,6 +15,8 @@ #include "utils/MPI_Info_object.h" #endif +#include "boundaries/BoundaryUtils.h" + DomainDecompBase::DomainDecompBase() : _rank(0), _numProcs(1) { } @@ -24,6 +26,32 @@ DomainDecompBase::~DomainDecompBase() { void DomainDecompBase::readXML(XMLfileUnits& /* xmlconfig */) { } +void DomainDecompBase::setGlobalBoundaryType(DimensionUtils::DimensionType dimension, BoundaryUtils::BoundaryType boundary) { + _boundaryHandler.setGlobalWallType(dimension, boundary); +} + +void DomainDecompBase::setLocalBoundariesFromGlobal(Domain* domain, Ensemble* ensemble) { + //find which walls to consider + double startRegion[3], endRegion[3]; + getBoundingBoxMinMax(domain, startRegion, endRegion); + const double* globStartRegion = ensemble->domain()->rmin(); + const double* globEndRegion = ensemble->domain()->rmax(); + + _boundaryHandler.setLocalRegion(startRegion, endRegion); + _boundaryHandler.setGlobalRegion(globStartRegion, globEndRegion); + _boundaryHandler.updateGlobalWallLookupTable(); +} + +void DomainDecompBase::processBoundaryConditions(ParticleContainer* moleculeContainer, double timestepLength) { + if(hasGlobalNonPeriodicBoundary()) + _boundaryHandler.processGlobalWallLeavingParticles(moleculeContainer, timestepLength); +} + +void DomainDecompBase::removeNonPeriodicHalos(ParticleContainer* moleculeContainer) { + if(hasGlobalNonPeriodicBoundary()) + _boundaryHandler.removeNonPeriodicHalos(moleculeContainer); +} + void DomainDecompBase::addLeavingMolecules(std::vector& invalidMolecules, ParticleContainer* moleculeContainer) { for (auto& molecule : invalidMolecules) { @@ -286,6 +314,11 @@ void DomainDecompBase::handleDomainLeavingParticlesDirect(const HaloRegion& halo } void DomainDecompBase::populateHaloLayerWithCopies(unsigned dim, ParticleContainer* moleculeContainer) const { + + //reflecting and outflow boundaries do not expect halo particles + if(_boundaryHandler.getGlobalWallType(dim) != BoundaryUtils::BoundaryType::PERIODIC) + return; + double shiftMagnitude = moleculeContainer->getBoundingBoxMax(dim) - moleculeContainer->getBoundingBoxMin(dim); // molecules that have crossed the lower boundary need a positive shift diff --git a/src/parallel/DomainDecompBase.h b/src/parallel/DomainDecompBase.h index b5fc83d871..b29ee19b1e 100644 --- a/src/parallel/DomainDecompBase.h +++ b/src/parallel/DomainDecompBase.h @@ -16,7 +16,8 @@ #include "molecules/MoleculeForwardDeclaration.h" #include "utils/Logger.h" // is this used? - +#include "boundaries/BoundaryHandler.h" +#include "ensemble/EnsembleBase.h" class Component; class Domain; @@ -289,6 +290,23 @@ class DomainDecompBase: public MemoryProfilable { virtual void printCommunicationPartners(std::string filename) const {}; + /* Set the global boundary type for the _boundaryHandler object. */ + void setGlobalBoundaryType(DimensionUtils::DimensionType dimension, BoundaryUtils::BoundaryType boundary); + + /* Find which boundaries of a subdomain are actually global boundaries, and update _boundaryHandler. */ + void setLocalBoundariesFromGlobal(Domain* domain, Ensemble* ensemble); + + /* Check if any of the global boundaries are invalid. */ + bool hasGlobalInvalidBoundary() const { return _boundaryHandler.hasGlobalInvalidBoundary();} + + bool hasGlobalNonPeriodicBoundary() const { return _boundaryHandler.hasGlobalNonPeriodicBoundary();} + + /* Processes leaving particles according to the boundary coundition of the wall the particles would be leaving. */ + void processBoundaryConditions(ParticleContainer* moleculeContainer, double timestepLength); + + /* Delete all halo particles outside global boundas that are non-periodic. */ + void removeNonPeriodicHalos(ParticleContainer* moleculeContainer); + protected: void addLeavingMolecules(std::vector& invalidMolecules, ParticleContainer* moleculeContainer); @@ -338,6 +356,8 @@ class DomainDecompBase: public MemoryProfilable { //! total number of processes in the simulation int _numProcs; + BoundaryHandler _boundaryHandler; + private: CollectiveCommBase _collCommBase; int _sendLeavingAndCopiesSeparately = 0; diff --git a/src/parallel/GeneralDomainDecomposition.cpp b/src/parallel/GeneralDomainDecomposition.cpp index ec656fbefe..46d14c60b1 100644 --- a/src/parallel/GeneralDomainDecomposition.cpp +++ b/src/parallel/GeneralDomainDecomposition.cpp @@ -134,6 +134,8 @@ void GeneralDomainDecomposition::balanceAndExchange(double lastTraversalTime, bo DomainDecompMPIBase::exchangeMoleculesMPI(moleculeContainer, domain, HALO_COPIES); } } + _boundaryHandler.setLocalRegion(_boxMin.data(),_boxMax.data()); + _boundaryHandler.updateGlobalWallLookupTable(); } ++_steps; } diff --git a/src/parallel/KDDecomposition.cpp b/src/parallel/KDDecomposition.cpp index 814cbd6fad..7a44ecd319 100644 --- a/src/parallel/KDDecomposition.cpp +++ b/src/parallel/KDDecomposition.cpp @@ -331,6 +331,11 @@ void KDDecomposition::balanceAndExchange(double lastTraversalTime, bool forceReb initCommunicationPartners(_cutoffRadius, domain, moleculeContainer); DomainDecompMPIBase::exchangeMoleculesMPI(moleculeContainer, domain, HALO_COPIES, true /*doHaloPositionCheck*/, removeRecvDuplicates); + + double startRegion[3], endRegion[3]; + getBoundingBoxMinMax(domain, startRegion, endRegion); + _boundaryHandler.setLocalRegion(startRegion,endRegion); + _boundaryHandler.updateGlobalWallLookupTable(); } } diff --git a/src/parallel/boundaries/BoundaryHandler.cpp b/src/parallel/boundaries/BoundaryHandler.cpp new file mode 100644 index 0000000000..2b02d6029e --- /dev/null +++ b/src/parallel/boundaries/BoundaryHandler.cpp @@ -0,0 +1,177 @@ +/* + * BoundaryHandler.cpp + * + * Created on: 24 March 2023 + * Author: amartyads + */ + +#include "BoundaryHandler.h" + +#include "integrators/Integrator.h" + +#include "utils/Logger.h" +#include "utils/Math.h" +#include "utils/mardyn_assert.h" + +#include +#include // for ostringstream + +BoundaryUtils::BoundaryType BoundaryHandler::getGlobalWallType(DimensionUtils::DimensionType dimension) const { + return _boundaries.at(dimension); +} + +BoundaryUtils::BoundaryType BoundaryHandler::getGlobalWallType(int dimension) const { + return getGlobalWallType(DimensionUtils::convertLS1DimIndexToEnumPositive(dimension)); +} + +void BoundaryHandler::setGlobalWallType(DimensionUtils::DimensionType dimension, BoundaryUtils::BoundaryType value) { + if (dimension != DimensionUtils::DimensionType::ERROR) + _boundaries[dimension] = value; + else { + std::ostringstream error_message; + error_message << "DimensionType::ERROR received in BoundaryHandler::setGlobalWallType!!" << std::endl; + MARDYN_EXIT(error_message.str()); + } +} + +void BoundaryHandler::setGlobalRegion(const double *start, const double *end) { + for (short int i = 0; i < 3; i++) { + _globalRegionStart[i] = start[i]; + _globalRegionEnd[i] = end[i]; + } +} + +void BoundaryHandler::setLocalRegion(const double *start, const double *end) { + for (short int i = 0; i < 3; i++) { + _localRegionStart[i] = start[i]; + _localRegionEnd[i] = end[i]; + } +} + +void BoundaryHandler::updateGlobalWallLookupTable() { + _isGlobalWall[DimensionUtils::DimensionType::POSX] = isNearRel(_localRegionEnd[0], _globalRegionEnd[0]); + _isGlobalWall[DimensionUtils::DimensionType::NEGX] = isNearRel(_localRegionStart[0], _globalRegionStart[0]); + _isGlobalWall[DimensionUtils::DimensionType::POSY] = isNearRel(_localRegionEnd[1], _globalRegionEnd[1]); + _isGlobalWall[DimensionUtils::DimensionType::NEGY] = isNearRel(_localRegionStart[1], _globalRegionStart[1]); + _isGlobalWall[DimensionUtils::DimensionType::POSZ] = isNearRel(_localRegionEnd[2], _globalRegionEnd[2]); + _isGlobalWall[DimensionUtils::DimensionType::NEGZ] = isNearRel(_localRegionStart[2], _globalRegionStart[2]); +} + +bool BoundaryHandler::hasGlobalInvalidBoundary() const { + return std::any_of(_boundaries.begin(), _boundaries.end(), [](const auto &keyVal) { + const auto [dim, boundaryType] = keyVal; + return boundaryType == BoundaryUtils::BoundaryType::ERROR; + }); +} + +bool BoundaryHandler::hasGlobalNonPeriodicBoundary() const { + return std::any_of(_boundaries.begin(), _boundaries.end(), [](const auto &keyVal) { + const auto [dim, boundaryType] = keyVal; + return boundaryType != BoundaryUtils::BoundaryType::PERIODIC; + }); +} + +bool BoundaryHandler::isGlobalWall(DimensionUtils::DimensionType dimension) const { + return _isGlobalWall.at(dimension); +} + +bool BoundaryHandler::isGlobalWall(int dimension) const { + return isGlobalWall(DimensionUtils::convertLS1DimIndexToEnumPositive(dimension)); +} + +void BoundaryHandler::processGlobalWallLeavingParticles(ParticleContainer *moleculeContainer, + double timestepLength) const { + const auto cutoff = moleculeContainer->getCutoff(); + for (auto const [currentDim, currentWallIsGlobalWall] : _isGlobalWall) { + if (!currentWallIsGlobalWall) + continue; + + switch (getGlobalWallType(currentDim)) { + case BoundaryUtils::BoundaryType::PERIODIC: + // nothing changes from normal ls1 behaviour, so leaving particles not touched by BoundaryHandler and + // are processed by DomainDecompBase::handleDomainLeavingParticles() + break; + + case BoundaryUtils::BoundaryType::OUTFLOW: + [[fallthrough]]; + case BoundaryUtils::BoundaryType::REFLECTING: { + // create region by using getInnerRegionSlab() + const auto [curWallRegionBegin, curWallRegionEnd] = + RegionUtils::getInnerRegionSlab(_localRegionStart, _localRegionEnd, currentDim, cutoff); + // grab an iterator from the converted coords + const auto particlesInRegion = moleculeContainer->regionIterator( + curWallRegionBegin.data(), curWallRegionEnd.data(), ParticleIterator::ONLY_INNER_AND_BOUNDARY); + + // iterate through all molecules + for (auto moleculeIter = particlesInRegion; moleculeIter.isValid(); ++moleculeIter) { + // Calculate the change in velocity, which the leapfrog method will + // apply in the next velocity update to the dimension of interest. + const int currentDimInt = DimensionUtils::convertEnumToLS1DimIndex(currentDim); + const double halfTimestep = .5 * timestepLength; + const double halfTimestepByMass = halfTimestep / moleculeIter->mass(); + const double force = moleculeIter->F(currentDimInt); + const double nextStepVelAdjustment = halfTimestepByMass * force; + + // check if the molecule would leave the bounds + if (RegionUtils::isMoleculeLeaving(*moleculeIter, curWallRegionBegin, curWallRegionEnd, currentDim, + timestepLength, nextStepVelAdjustment)) { + if (getGlobalWallType(currentDim) == BoundaryUtils::BoundaryType::REFLECTING) { + const double currentVel = moleculeIter->v(currentDimInt); + // change the velocity in the dimension of interest such that when + // the leapfrog integrator adds nextStepVelAdjustment in the next + // velocity update, the final result ends up being the intended, + // reversed velocity: -(currentVel+nextStepVelAdjustment) + moleculeIter->setv(currentDimInt, + -currentVel - nextStepVelAdjustment - nextStepVelAdjustment); + } else { // outflow, delete the particle if it would leave + moleculeContainer->deleteMolecule(moleculeIter, false); + } + } + } + break; + } + default: + std::ostringstream error_message; + error_message << "BoundaryType::ERROR received in BoundaryHandler::processGlobalWallLeavingParticles!" << std::endl; + MARDYN_EXIT(error_message.str()); + } + } +} + +void BoundaryHandler::removeNonPeriodicHalos(ParticleContainer *moleculeContainer) const { + // get halo lengths in each dimension + const std::array haloWidths = {moleculeContainer->getHaloWidthForDimension(0), + moleculeContainer->getHaloWidthForDimension(1), + moleculeContainer->getHaloWidthForDimension(2)}; + for (auto const [currentDim, currentWallIsGlobalWall] : _isGlobalWall) { + if (!currentWallIsGlobalWall) + continue; + + switch (getGlobalWallType(currentDim)) { + case BoundaryUtils::BoundaryType::PERIODIC: + // nothing changes from normal ls1 behaviour, so empty case, and halo particles left untouched + break; + + case BoundaryUtils::BoundaryType::OUTFLOW: + [[fallthrough]]; + case BoundaryUtils::BoundaryType::REFLECTING: { + // create region by using getOuterRegionSlab() + auto const [curWallRegionBegin, curWallRegionEnd] = + RegionUtils::getOuterRegionSlab(_localRegionStart, _localRegionEnd, currentDim, haloWidths); + + // grab an iterator from the converted coords + auto particlesInRegion = moleculeContainer->regionIterator( + curWallRegionBegin.data(), curWallRegionEnd.data(), ParticleIterator::ALL_CELLS); + for (auto moleculeIter = particlesInRegion; moleculeIter.isValid(); ++moleculeIter) { + // delete all halo particles + moleculeContainer->deleteMolecule(moleculeIter, false); + } + break; + } + default: + std::ostringstream error_message; + error_message << "BoundaryType::ERROR received in BoundaryHandler::removeNonPeriodicHalos!" << std::endl; + MARDYN_EXIT(error_message.str()); + } + } +} diff --git a/src/parallel/boundaries/BoundaryHandler.h b/src/parallel/boundaries/BoundaryHandler.h new file mode 100644 index 0000000000..08bcea4945 --- /dev/null +++ b/src/parallel/boundaries/BoundaryHandler.h @@ -0,0 +1,193 @@ +/* + * BoundaryHandler.h + * + * Created on: 24 March 2023 + * Author: amartyads + */ + +#pragma once + +#include "BoundaryUtils.h" +#include "DimensionUtils.h" +#include "RegionUtils.h" +#include "particleContainer/ParticleContainer.h" + +#include +#include +#include +#include + +/** + * @brief Class to handle boundary conditions, namely leaving and halo particles. + * @author Amartya Das Sharma + * + * The objects of this class store the local and global bounds of the subdomain in every process, and provide functions + * to deal with leaving particles, and delete halo particles. + * + * The internal walls of the subdomain, touching other subdomains are 'local' walls while the walls that are also the + * limits of the global domain are 'global' walls. + * + * All subdomains have a copy of the global wall types, and do not store their local boundary conditions; instead they + * check whether a particular local wall is also a global wall, and then use the global lookup table to ascertain what + * boundary effects to use. + * + * The default state for all global walls is 'PERIODIC', since that is the default ls1 behaviour. + */ +class BoundaryHandler { +public: + BoundaryHandler() = default; + + /** + * @brief Find the boundary type of a global wall for a particular dimension. + * + * Since this returns the global wall type, every subdomain on every rank will return + * the same value for the same input. + * + * @param dimension the dimension of interest, must be DimensionType enum member + * @return BoundaryUtils::BoundaryType enum member corresponding to the global wall type + */ + BoundaryUtils::BoundaryType getGlobalWallType(DimensionUtils::DimensionType dimension) const; + + /** + * @brief Find the boundary type of a global wall for a particular dimension. + * + * Since this returns the global wall type, every subdomain on every rank will return + * the same value for the same input. + * + * @param dimension tte dimension of interest, must be an ls1-style index (0 for x, 1 for y, 2 for z) + * @return BoundaryUtils::BoundaryType enum member corresponding to the global wall type + */ + BoundaryUtils::BoundaryType getGlobalWallType(int dimension) const; + + /** + * @brief Set the boundary type of a global wall for a particular dimension. + * + * @param dimension the dimension of interest, must be DimensionType enum member + * @param value the type of the boundary, must be BoundaryType enum member + */ + void setGlobalWallType(DimensionUtils::DimensionType dimension, BoundaryUtils::BoundaryType value); + + /** + * @brief Check if any of the global boundaries have invalid types. + * + * @return true if any of the global boundaries is BoundaryType::ERROR + * @return false if none of the global boundaries are BoundaryType::ERROR + */ + bool hasGlobalInvalidBoundary() const; + + /** + * @brief Check if any of the global boundaries are non-periodic. + * + * This check helps bypass all boundary-handling related code if default behaviour + * (all periodic boundaries) is found. + * + * @return true if any of the global boundaries are non-periodic + * @return false if all of the global boundaries are periodic + */ + bool hasGlobalNonPeriodicBoundary() const; + + /** + * @brief Set bounds for global domain. + * + * @param start double[3] array with coordinates of the start point + * @param end double[3] array with coordinates of the end point + */ + void setGlobalRegion(const double *start, const double *end); + + /** + * @brief Set bounds for local subdomain. + * + * @param start double[3] array with coordinates of the start point + * @param end double[3] array with coordinates of the end point + */ + void setLocalRegion(const double *start, const double *end); + + /** + * @brief Determine which walls in the local region are actually global walls. + * + * Should be called after changing global and local regions (typically after a rebalance). + */ + void updateGlobalWallLookupTable(); + + /** + * @brief Check if the local wall in a particular dimension is actually a global wall. + * + * @param dimension the dimension of interest, must be DimensionType enum member + * @return true if the local wall in the given dimension is a global wall + * @return false if the local wall in the given dimension is not a global wall + */ + bool isGlobalWall(DimensionUtils::DimensionType dimension) const; + + /** + * @brief Check if the local wall in a particular dimension is actually a global wall. + * + * @param dimension the dimension of interest, must be an ls1-style index (0 for x, 1 for y, 2 for z) + * @return true if the local wall in the given dimension is a global wall + * @return false if the local wall in the given dimension is not a global wall + */ + bool isGlobalWall(int dimension) const; + + /** + * @brief Processes all particles that would leave the global domain. + * + * If a subdomain has no global walls, this function does nothing. + * For every global wall, the function iterates through all particles that are + * within one cutoff distance away from the wall. If these particles would + * leave the global box in the next simulation, the following is done: + * + * PERIODIC - No actions taken (default behaviour). + * REFLECTING - The particle's velocity is reversed normal to the wall it's leaving. + * OUTFLOW - The particle is deleted. + * + * If a particle would hit multiple walls with multiple types, the following happens: + * - If any of the walls is an outflow wall, the particle is deleted immediately, otherwise the following happen + * - The particle's velocity is reversed for the components where the wall is reflecting + * - Periodic effects happen as normal + * + * @param moleculeContainer used to get the cutoff, region iterators, and to delete particles if needed + * @param timestepLength the full timestep length, used to determine the position of the molecule + * in the next timestep to check whether it is leaving + */ + void processGlobalWallLeavingParticles(ParticleContainer *moleculeContainer, double timestepLength) const; + + /** + * @brief Processes all halo particles outside the global domain. + * + * If a subdomain has no global walls, this function does nothing. + * For every global wall, the function iterates through all halo particles + * that are within one cutoff distance away from the wall. The following is + * done for each particle: + * + * PERIODIC - No actions taken (default behaviour). + * REFLECTING / OUTFLOW - The halo particle is deleted, so that particles + * approaching the boundary do not decelerate due to influence from the halo + * particles, and preserve their velocities before being bounced/deleted + * + * @param moleculeContainer used to get the cutoff, region iterators, and to delete particles if needed + */ + void removeNonPeriodicHalos(ParticleContainer *moleculeContainer) const; + +private: + /** + * @brief Lookup table for global boundary type in every dimension. + * + * Set as periodic by default, since the default behaviour of LS1 is all periodic boundaries. + * This allows the tag to be optional. + */ + std::map _boundaries{ + {DimensionUtils::DimensionType::POSX, BoundaryUtils::BoundaryType::PERIODIC}, + {DimensionUtils::DimensionType::POSY, BoundaryUtils::BoundaryType::PERIODIC}, + {DimensionUtils::DimensionType::POSZ, BoundaryUtils::BoundaryType::PERIODIC}, + {DimensionUtils::DimensionType::NEGX, BoundaryUtils::BoundaryType::PERIODIC}, + {DimensionUtils::DimensionType::NEGY, BoundaryUtils::BoundaryType::PERIODIC}, + {DimensionUtils::DimensionType::NEGZ, BoundaryUtils::BoundaryType::PERIODIC}}; + + /* Lookup table to check if a local wall is also global. */ + std::map _isGlobalWall; + + /* Global region start/end. */ + std::array _globalRegionStart, _globalRegionEnd; + + /* Local region start/end. */ + std::array _localRegionStart, _localRegionEnd; +}; diff --git a/src/parallel/boundaries/BoundaryUtils.cpp b/src/parallel/boundaries/BoundaryUtils.cpp new file mode 100644 index 0000000000..747f39d2ae --- /dev/null +++ b/src/parallel/boundaries/BoundaryUtils.cpp @@ -0,0 +1,48 @@ +/* + * BoundaryUtils.cpp + * + * Created on: 24 March 2023 + * Author: amartyads + */ + +#include "BoundaryUtils.h" + +#include "utils/Logger.h" +#include "utils/mardyn_assert.h" + +#include +#include +#include +#include // for ostringstream + +BoundaryUtils::BoundaryType BoundaryUtils::convertStringToBoundary(const std::string &boundary) { + std::string boundaryLowercase(boundary); + std::transform(boundaryLowercase.begin(), boundaryLowercase.end(), boundaryLowercase.begin(), + [](unsigned char c) { return std::tolower(c); }); + if (boundaryLowercase.find("per") != std::string::npos) + return BoundaryType::PERIODIC; + if (boundaryLowercase.find("ref") != std::string::npos) + return BoundaryType::REFLECTING; + if (boundaryLowercase.find("out") != std::string::npos) + return BoundaryType::OUTFLOW; + std::ostringstream error_message; + error_message << "Invalid boundary type passed to BoundaryUtils::convertStringToBoundary. Check your input file!" << std::endl; + MARDYN_EXIT(error_message.str()); + return BoundaryType::ERROR; // warning suppression +} + +std::string BoundaryUtils::convertBoundaryToString(BoundaryType boundary) { + switch (boundary) { + case BoundaryType::PERIODIC: + return "periodic"; + case BoundaryType::REFLECTING: + return "reflecting"; + case BoundaryType::OUTFLOW: + return "outflow"; + default: + std::ostringstream error_message; + error_message << "BoundaryType::ERROR received in BoundaryUtils::convertBoundaryToString!" << std::endl; + MARDYN_EXIT(error_message.str()); + } + return "error"; // warning suppression +} diff --git a/src/parallel/boundaries/BoundaryUtils.h b/src/parallel/boundaries/BoundaryUtils.h new file mode 100644 index 0000000000..f5f039b74b --- /dev/null +++ b/src/parallel/boundaries/BoundaryUtils.h @@ -0,0 +1,55 @@ +/* + * BoundaryUtils.h + * + * Created on: 24 March 2023 + * Author: amartyads + */ + +#pragma once + +#include + +/** + * @brief Includes the BoundaryType enum and helper functions for boundary types. + * @author Amartya Das Sharma + */ +namespace BoundaryUtils { + +/** + * @brief enum storing the types of global boundaries currently supported. + * + * The currently supported boundary types are + * PERIODIC - periodic boundaries, using the default (communicate particle transfers with neighbours) behaviour + * OUTFLOW - molecules exiting the boundary are deleted + * REFLECTING - molecules exiting the boundary have their velocities reversed in the direction they are leaving (i.e. + * reflected) ERROR - kept for sanity checks + * + * This can be extended if needed. + */ +enum class BoundaryType { PERIODIC, OUTFLOW, REFLECTING, ERROR }; + +/** + * @brief Convert strings read from the XML input file into BoundaryType enum members. + * + * The conversion logic is as follows: + * - Any string with "per" - BoundaryType::PERIODIC + * - Any string with "ref" - BoundaryType::REFLECTING + * - Any string with "out" - BoundaryType::OUTFLOW + * An error occurs if none of the substrings are found. + * + * @param boundary string read from xml file + * @return BoundaryType which is the corresponding BoundaryType enum member + */ +BoundaryType convertStringToBoundary(const std::string &boundary); + +/** + * @brief Converts BoundaryType enum members into strings. + * + * Used mainly for logging. + * + * @param boundary BoundaryType enum member + * @return std::string corresponding to the enum member + */ +std::string convertBoundaryToString(BoundaryType boundary); + +} // namespace BoundaryUtils diff --git a/src/parallel/boundaries/DimensionUtils.cpp b/src/parallel/boundaries/DimensionUtils.cpp new file mode 100644 index 0000000000..9788512f1a --- /dev/null +++ b/src/parallel/boundaries/DimensionUtils.cpp @@ -0,0 +1,113 @@ +/* + * DimensionUtils.cpp + * + * Created on: 18 Sep 2024 + * Author: amartyads + */ + +#include "DimensionUtils.h" +#include "utils/Logger.h" +#include "utils/mardyn_assert.h" // for MARDYN_EXIT() + +#include // for ostringstream + +bool DimensionUtils::isDimensionNumericPermissible(int dim) { + // permissible dimensions are {-1, -2, -3, 1, 2, 3} + return (dim >= -3 && dim <= 3 && dim != 0); +} + +DimensionUtils::DimensionType DimensionUtils::convertNumericToDimension(int dim) { + if (!isDimensionNumericPermissible(dim)) { + std::ostringstream error_message; + error_message << "Invalid dimension passed for enum conversion. Received value: " << dim + << std::endl; + MARDYN_EXIT(error_message.str()); + return DimensionType::ERROR; + } + switch (dim) { + case 1: + return DimensionType::POSX; + case 2: + return DimensionType::POSY; + case 3: + return DimensionType::POSZ; + case -1: + return DimensionType::NEGX; + case -2: + return DimensionType::NEGY; + case -3: + return DimensionType::NEGZ; + } + return DimensionType::ERROR; // warning suppression +} + +DimensionUtils::DimensionType DimensionUtils::convertLS1DimIndexToEnumPositive(int dim) { + switch (dim) { + case 0: + return DimensionType::POSX; + case 1: + return DimensionType::POSY; + case 2: + return DimensionType::POSZ; + default: + return DimensionType::ERROR; + } +} + +std::string DimensionUtils::convertDimensionToString(DimensionType dimension) { + switch (dimension) { + case DimensionType::POSX: + return "+x"; + case DimensionType::POSY: + return "+y"; + case DimensionType::POSZ: + return "+z"; + case DimensionType::NEGX: + return "-x"; + case DimensionType::NEGY: + return "-y"; + case DimensionType::NEGZ: + return "-z"; + default: // ERROR + std::ostringstream error_message; + error_message << "DimesionType::ERROR received in DimensionUtils::convertDimensionToString!" + << std::endl; + MARDYN_EXIT(error_message.str()); + return "error"; + } +} + +std::string DimensionUtils::convertDimensionToStringAbs(DimensionType dimension) { + return convertDimensionToString(dimension).substr(1, 1); +} + +int DimensionUtils::convertDimensionToNumeric(DimensionType dimension) { + switch (dimension) { + case DimensionType::POSX: + return 1; + case DimensionType::POSY: + return 2; + case DimensionType::POSZ: + return 3; + case DimensionType::NEGX: + return -1; + case DimensionType::NEGY: + return -2; + case DimensionType::NEGZ: + return -3; + default: + std::ostringstream error_message; + error_message << "DimesionType::ERROR received in DimensionUtils::convertDimensionToNumeric!" + << std::endl; + MARDYN_EXIT(error_message.str()); + return 0; + } +} + +int DimensionUtils::convertDimensionToNumericAbs(DimensionType dimension) { + return std::abs(convertDimensionToNumeric(dimension)); +} + +int DimensionUtils::convertEnumToLS1DimIndex(DimensionType dimension) { + return convertDimensionToNumericAbs(dimension) - 1; +} diff --git a/src/parallel/boundaries/DimensionUtils.h b/src/parallel/boundaries/DimensionUtils.h new file mode 100644 index 0000000000..2b6fb0ba0a --- /dev/null +++ b/src/parallel/boundaries/DimensionUtils.h @@ -0,0 +1,123 @@ +/* + * DimensionUtils.h + * + * Created on: 18 Sep 2024 + * Author: amartyads + */ + +#pragma once + +#include "utils/Math.h" +#include + +/** + * @brief Includes the DimensionType enum and helper functions for dimension types. + * @author Amartya Das Sharma + */ +namespace DimensionUtils { + +/** + * @brief enum storing the axes and direction. + * + * The dimensions are POSX, NEGX, POSY, NEGY, POSZ and NEGZ. + * + * This is hardcoded for 3D, and ERROR is included for sanity checks. + */ +enum class DimensionType { POSX, NEGX, POSY, NEGY, POSZ, NEGZ, ERROR }; + +/** + * @brief Checks if a numeric dimension is allowed. + * + * Allowed numeric dimensions are +-1 for x, +-2 for y, +-3 for z + * + * @param dim integer dimension + * @return true if `dim` is one of (1, 2, 3, -1, -2, -3) + * @return false if `dim` is not one of (1, 2, 3, -1, -2, -3) + */ +bool isDimensionNumericPermissible(int dim); + +/** + * @brief Converts a dimension from number to DimensionType, where x = +-1, y = +-2 and z = +-3. + * + * Throws an error and exits if the dimension is not permissible. + * + * @param dim integer dimension + * @return DimensionType enum member corresponding to the integer dimension + */ +DimensionType convertNumericToDimension(int dim); + +/** + * @brief Converts an LS1 internal dimension representation int to DimensionType, where + * x = 0, y = 1 and z = 2 Since this does not contain direction information, the + * positive direction is returned. + * + * Throws an error and exits if the dimension is not permissible. + * + * @param dim integer dimension (either 0, 1 or 2) + * @return DimensionType enum member corresponding to the integer dimension + */ +DimensionType convertLS1DimIndexToEnumPositive(int dim); + +/** + * @brief Converts a DimensionType into a string. + * + * The expected return values are +x, -x, +y, -y, +z, -z. Can be used for logging purposes. + * Throws an error and exits if DimensionType::ERROR is encountered. + * + * @param dimension DimensionType enum member + * @return std::string which can be +x, -x, +y, -y, +z, -z + */ +std::string convertDimensionToString(DimensionType dimension); + +/** + * @brief Converts a DimensionType into a string, and remove directional information. + * + * POSX and NEGX both return "x", POSY and NEGY return "y", and POSZ and NEGZ return "z". + * Can be used for logging purposes. Throws an error and exits if DimensionType::ERROR is encountered. + * + * @param dimension DimensionType enum member + * @return std::string which can be x, y, or z + */ +std::string convertDimensionToStringAbs(DimensionType dimension); + +/** + * @brief Converts a dimension from DimensionType to number, where x = +-1, y = +-2 and + * z = +-3. Throws an error and exits if DimensionType::ERROR is encountered. + * + * @param dimension DimensionType enum member + * @return int which is one of (1, 2, 3, -1, -2, -3) + */ +int convertDimensionToNumeric(DimensionType dimension); + +/** + * @brief Converts a DimensionType into an int, and remove directions. + * + * +x/-x returns 1, +y/-y returns 2, +z/-z returns 3. + * Throws an error and exits if DimensionType::ERROR is encountered. + * + * @param dimension DimensionType enum member + * @return int which is one of (1, 2, 3) + */ +int convertDimensionToNumericAbs(DimensionType dimension); + +/** + * @brief Converts a DimensionType to LS1 internal dimension representation int, where x + * = 0, y = 1 and z = 2 + * + * Since this does not contain direction information, both POSX and NEGX return 0, for example. + * Throws an error and exits if DimensionType::ERROR is encountered. + * + * @param dimension DimensionType enum member + * @return int which is one of (0, 1, 2) + */ +int convertEnumToLS1DimIndex(DimensionType dimension); + +/** + * @brief Returns the direction from a given DimensionType. + * + * @param dimension DimensionType enum member + * @return int which is 1 if the direction is positive, and -1 otherwise + */ +inline int findSign(DimensionType dimension) { return ::findSign(convertDimensionToNumeric(dimension)); } + +} // namespace DimensionUtils diff --git a/src/parallel/boundaries/RegionUtils.cpp b/src/parallel/boundaries/RegionUtils.cpp new file mode 100644 index 0000000000..337ffc864b --- /dev/null +++ b/src/parallel/boundaries/RegionUtils.cpp @@ -0,0 +1,117 @@ +/* + * RegionUtils.cpp + * + * Created on: 18 Sep 2024 + * Author: amartyads + */ + +#include "RegionUtils.h" +#include "utils/mardyn_assert.h" //for MARDYN_EXIT() + +#include // for ostringstream + +std::tuple, std::array> RegionUtils::getInnerRegionSlab( + const std::array &givenRegionBegin, const std::array &givenRegionEnd, + DimensionUtils::DimensionType dimension, double regionWidth) { + + std::array returnRegionBegin = givenRegionBegin; + std::array returnRegionEnd = givenRegionEnd; + + const int dimensionLS1 = convertEnumToLS1DimIndex(dimension); + switch (dimension) { + // in positive case, set the beginning to end-width, or whole domain if width + // too large + case DimensionUtils::DimensionType::POSX: + [[fallthrough]]; + case DimensionUtils::DimensionType::POSY: + [[fallthrough]]; + case DimensionUtils::DimensionType::POSZ: + returnRegionBegin[dimensionLS1] = + std::max(returnRegionEnd[dimensionLS1] - regionWidth, givenRegionBegin[dimensionLS1]); + break; + // in negative case, set the end to beginning+width, or whole domain if width + // too large + case DimensionUtils::DimensionType::NEGX: + [[fallthrough]]; + case DimensionUtils::DimensionType::NEGY: + [[fallthrough]]; + case DimensionUtils::DimensionType::NEGZ: + returnRegionEnd[dimensionLS1] = + std::min(returnRegionBegin[dimensionLS1] + regionWidth, givenRegionEnd[dimensionLS1]); + break; + + default: + std::ostringstream error_message; + error_message << "DimensionType::ERROR received in RegionUtils::getInnerRegionSlab" << std::endl; + MARDYN_EXIT(error_message.str()); + } + return {returnRegionBegin, returnRegionEnd}; +} + +bool RegionUtils::isMoleculeLeaving(const Molecule &molecule, const std::array ®ionBegin, + const std::array ®ionEnd, DimensionUtils::DimensionType dimension, + double timestepLength, double nextStepVelAdjustment) { + const int ls1dim = convertEnumToLS1DimIndex(dimension); + const int direction = findSign(dimension); + const double newPos = molecule.r(ls1dim) + (timestepLength * (molecule.v(ls1dim) + nextStepVelAdjustment)); + if (newPos <= regionBegin[ls1dim] && direction < 0) + return true; + if (newPos >= regionEnd[ls1dim] && direction > 0) + return true; + return false; +} + +std::tuple, std::array> RegionUtils::getOuterRegionSlab( + const std::array &givenRegionBegin, const std::array &givenRegionEnd, + DimensionUtils::DimensionType dimension, const std::array ®ionWidth) { + std::array returnRegionBegin = givenRegionBegin; + std::array returnRegionEnd = givenRegionEnd; + + for (int i = 0; i < 3; i++) { + returnRegionBegin[i] = givenRegionBegin[i]; + returnRegionEnd[i] = givenRegionEnd[i]; + } + + const int dimensionLS1 = convertEnumToLS1DimIndex(dimension); + + // find the two dimensions that are not being considered + const int extraDim1 = dimensionLS1 == 0 ? 1 : 0; + const int extraDim2 = dimensionLS1 == 2 ? 1 : 2; + + // extend the extra dimensions to cover all ghost areas + returnRegionBegin[extraDim1] -= regionWidth[extraDim1]; + returnRegionEnd[extraDim1] += regionWidth[extraDim1]; + + returnRegionBegin[extraDim2] -= regionWidth[extraDim2]; + returnRegionEnd[extraDim2] += regionWidth[extraDim2]; + + switch (dimension) { + // in positive case, move the box begin to edge of domain, and box end to + // beyond + case DimensionUtils::DimensionType::POSX: + [[fallthrough]]; + case DimensionUtils::DimensionType::POSY: + [[fallthrough]]; + case DimensionUtils::DimensionType::POSZ: + returnRegionBegin[dimensionLS1] = returnRegionEnd[dimensionLS1]; + returnRegionEnd[dimensionLS1] += regionWidth[dimensionLS1]; + break; + + // in negative case, move the box end to edge of domain, and box begin to + // beyond + case DimensionUtils::DimensionType::NEGX: + [[fallthrough]]; + case DimensionUtils::DimensionType::NEGY: + [[fallthrough]]; + case DimensionUtils::DimensionType::NEGZ: + returnRegionEnd[dimensionLS1] = returnRegionBegin[dimensionLS1]; + returnRegionBegin[dimensionLS1] -= regionWidth[dimensionLS1]; + break; + + default: + std::ostringstream error_message; + error_message << "DimensionType::ERROR received in RegionUtils::getOuterRegionSlab" << std::endl; + MARDYN_EXIT(error_message.str()); + } + return std::make_tuple(returnRegionBegin, returnRegionEnd); +} diff --git a/src/parallel/boundaries/RegionUtils.h b/src/parallel/boundaries/RegionUtils.h new file mode 100644 index 0000000000..6f613178c7 --- /dev/null +++ b/src/parallel/boundaries/RegionUtils.h @@ -0,0 +1,99 @@ +/* + * RegionUtils.h + * + * Created on: 18 Sep 2024 + * Author: amartyads + */ + +#pragma once + +#include "DimensionUtils.h" +#include "molecules/Molecule.h" +#include + +/** + * @brief Includes helper functions for region-related computation during boundary handling. + * @author Amartya Das Sharma + * + * Does not actually work on any molecule containers, hence can be used + * in a general sense for other calculations. + */ +namespace RegionUtils { + +/** + * @brief When given a domain delimited by givenRegionBegin and givenRegionEnd + * and a DimensionType, and a requested regionWidth, this function + * returns a cropped box from the domain. This box is of size regionWidth + * in the dimension specified, and the size of the domain otherwise. + * + * E.g. if the dimension given is +x, and the regionWidth is 3, this function + * returns a box, which would contain every particle that is within 3 + * units from the +x boundary wall, in the original domain + * demarcated by givenRegion{begin,end}. + * + * @param givenRegionBegin the start of the region to be cropped + * @param givenRegionEnd the end of the region to be cropped + * @param dimension the dimension which is cropped i.e. the resulting slab will be regionWidth wide in this dimension, + * and unchanged otherwise. + * @param regionWidth the width of the slab + * @return std::tuple, std::array> the corners of the resultant slab inside the given + * region + */ +std::tuple, std::array> getInnerRegionSlab( + const std::array &givenRegionBegin, const std::array &givenRegionEnd, + DimensionUtils::DimensionType dimension, double regionWidth); + +/** + * @brief Checks if a given molecule is leaving the given region in the next timestep, + * with added scalar velocity adjustment, but only in the given dimension. + * + * The molecule's position and velocity in the specified dimension is extracted + * from the molecule itself, and the next position is calculated using the + * velocity adjustment and the timestep length. + * + * If the molecule was already out of the box (and stay out of the box after position adjustment), this function returns + * true. If the molecule somehow ends up entering the box after the position adjustment, the function returns false. + * + * When this function is called, it is usually expected that the forces have + * been updated for the current timestep and the positions and velocities + * aren't. Hence the nextStepVelAdjustment is provided to the caller function, + * so that it can handle the change in velocity due to forces. This is not done + * by the util function itself, to keep it as generic as possible. + * + * @param molecule the molecule to be checked + * @param regionBegin the start of the region to be checked + * @param regionEnd the end of the region to be checked + * @param dimension the dimension in which the molecule's new position must be calculated, to check for leaving + * @param timestepLength the length of the whole time step, after which the molecule's position is calculated + * @param nextStepVelAdjustment a flat scalar velocity added to the velocity in the `dimension` component after + * calculations + * @return true if the molecule would leave the box after timestepLength + * @return false if the molecule would not leave the box after timestepLength + */ +bool isMoleculeLeaving(const Molecule &molecule, const std::array ®ionBegin, + const std::array ®ionEnd, DimensionUtils::DimensionType dimension, + double timestepLength, double nextStepVelAdjustment); + +/** + * @brief When given a domain delimited by givenRegionBegin and givenRegionEnd + * and a DimensionType, and a requested regionWidth, this function + * returns a box outside the domain. This box is of size regionWidth + * in the dimension specified, and the size of the domain plus regionwidth + * otherwise. + * + * E.g. if the dimension given is +x, and the regionWidth is 3, this function + * returns a box, which contains every ghost particle that is 3 units away + * from +x. + * + * @param givenRegionBegin the start of the region under consideration + * @param givenRegionEnd the end of the region under consideration + * @param dimension the dimension in which the box should be sharing a wall with the givenRegion + * @param regionWidth the width of the slab + * @return std::tuple, std::array> the corners of the resultant region outside the + * given region + */ +std::tuple, std::array> getOuterRegionSlab( + const std::array &givenRegionBegin, const std::array &givenRegionEnd, + DimensionUtils::DimensionType dimension, const std::array ®ionWidth); + +} // namespace RegionUtils diff --git a/src/particleContainer/AutoPasContainer.cpp b/src/particleContainer/AutoPasContainer.cpp index d445392748..b28bf3a668 100644 --- a/src/particleContainer/AutoPasContainer.cpp +++ b/src/particleContainer/AutoPasContainer.cpp @@ -633,7 +633,7 @@ void AutoPasContainer::deleteOuterParticles() { } } -double AutoPasContainer::get_halo_L(int /*index*/) const { return _cutoff; } +double AutoPasContainer::getHaloWidthForDimension(int /*index*/) const { return _cutoff; } double AutoPasContainer::getCutoff() const { return _cutoff; } diff --git a/src/particleContainer/AutoPasContainer.h b/src/particleContainer/AutoPasContainer.h index 5167ba5818..ebdddf686c 100644 --- a/src/particleContainer/AutoPasContainer.h +++ b/src/particleContainer/AutoPasContainer.h @@ -89,7 +89,7 @@ class AutoPasContainer : public ParticleContainer { void deleteOuterParticles() override; - double get_halo_L(int index) const override; + double getHaloWidthForDimension(int index) const override; double getCutoff() const override; diff --git a/src/particleContainer/LinkedCells.cpp b/src/particleContainer/LinkedCells.cpp index 9087b21d74..eab3591bea 100644 --- a/src/particleContainer/LinkedCells.cpp +++ b/src/particleContainer/LinkedCells.cpp @@ -629,7 +629,7 @@ void LinkedCells::deleteOuterParticles() { } } -double LinkedCells::get_halo_L(int index) const { +double LinkedCells::getHaloWidthForDimension(int index) const { return _haloLength[index]; } diff --git a/src/particleContainer/LinkedCells.h b/src/particleContainer/LinkedCells.h index a444fea3ae..add747275e 100644 --- a/src/particleContainer/LinkedCells.h +++ b/src/particleContainer/LinkedCells.h @@ -190,7 +190,7 @@ class LinkedCells : public ParticleContainer { //! @brief gets the width of the halo region in dimension index //! @todo remove this method, because a halo_L shouldn't be necessary for every ParticleContainer //! e.g. replace it by the cutoff-radius - double get_halo_L(int index) const override; + double getHaloWidthForDimension(int index) const override; double getCutoff() const override { return _cutoffRadius; } void setCutoff(double rc) override { _cutoffRadius = rc; } diff --git a/src/particleContainer/ParticleContainer.h b/src/particleContainer/ParticleContainer.h index eea8af21cc..012b68f2b3 100644 --- a/src/particleContainer/ParticleContainer.h +++ b/src/particleContainer/ParticleContainer.h @@ -192,7 +192,7 @@ class ParticleContainer: public MemoryProfilable { //! @brief returns the width of the halo stripe (for the given dimension index) //! @todo remove this method, because a halo_L shouldn't be necessary for every ParticleContainer //! e.g. replace it by the cutoff-radius - virtual double get_halo_L(int index) const = 0; + virtual double getHaloWidthForDimension(int index) const = 0; virtual double getCutoff() const = 0; diff --git a/src/plugins/MamicoCoupling.cpp b/src/plugins/MamicoCoupling.cpp index 4bb5739eef..c652a16be4 100644 --- a/src/plugins/MamicoCoupling.cpp +++ b/src/plugins/MamicoCoupling.cpp @@ -1,64 +1,32 @@ +#ifdef MAMICO_COUPLING + #include "MamicoCoupling.h" #include "Domain.h" -void MamicoCoupling::readXML(XMLfileUnits& xmlconfig) -{ - return; -} - -void MamicoCoupling::init(ParticleContainer* particleContainer, - DomainDecompBase* domainDecomp, Domain* domain) -{ -#ifdef MAMICO_COUPLING - //code to print to log that plugin is initialised +void MamicoCoupling::init(ParticleContainer *particleContainer, DomainDecompBase *domainDecomp, Domain *domain) { Log::global_log->info() << "MaMiCo coupling plugin initialized" << std::endl; -#endif -} - -void MamicoCoupling::beforeEventNewTimestep(ParticleContainer* particleContainer, - DomainDecompBase* domainDecomp, unsigned long simstep) -{ - } -void MamicoCoupling::beforeForces(ParticleContainer* particleContainer, - DomainDecompBase* domainDecomp, unsigned long simstep) -{ -#ifdef MAMICO_COUPLING - if(_couplingEnabled) - { - // This object should be set by MaMiCo after the plugins are created in the simulation readxml file +void MamicoCoupling::beforeForces(ParticleContainer *particleContainer, DomainDecompBase *domainDecomp, + unsigned long simstep) { + if (_couplingEnabled) { + // This object should be set by MaMiCo after the plugins are created in the simulation readxml file. // Even though this method is called before the object is set, at this point the coupling switch is always off - _macroscopicCellService->processInnerMacroscopicCellAfterMDTimestep(); - _macroscopicCellService->distributeMass(simstep); - _macroscopicCellService->applyTemperatureToMolecules(simstep); + _couplingCellService->processInnerCouplingCellAfterMDTimestep(); + _couplingCellService->distributeMass(simstep); + _couplingCellService->applyTemperatureToMolecules(simstep); #ifndef MARDYN_AUTOPAS particleContainer->deleteOuterParticles(); #endif } -#endif } -void MamicoCoupling::afterForces(ParticleContainer* particleContainer, - DomainDecompBase* domainDecomp, unsigned long simstep) -{ -#ifdef MAMICO_COUPLING - if(_couplingEnabled) - { - _macroscopicCellService->distributeMomentum(simstep); - _macroscopicCellService->applyBoundaryForce(simstep); +void MamicoCoupling::afterForces(ParticleContainer *particleContainer, DomainDecompBase *domainDecomp, + unsigned long simstep) { + if (_couplingEnabled) { + _couplingCellService->distributeMomentum(simstep); + _couplingCellService->applyBoundaryForce(simstep); } -#endif } -void MamicoCoupling::endStep(ParticleContainer* particleContainer, - DomainDecompBase* domainDecomp, Domain* domain, unsigned long simstep) -{ - -} - -void MamicoCoupling::finish(ParticleContainer* particleContainer, - DomainDecompBase* domainDecomp, Domain* domain) -{ - -} +#endif // MAMICO_COUPLING diff --git a/src/plugins/MamicoCoupling.h b/src/plugins/MamicoCoupling.h index c8a7c2f5c2..1512edb15c 100644 --- a/src/plugins/MamicoCoupling.h +++ b/src/plugins/MamicoCoupling.h @@ -1,115 +1,119 @@ -#pragma once - -#include "PluginBase.h" #ifdef MAMICO_COUPLING -#include -#include -#endif -/** @brief Allows execution of MaMiCo code to enable coupling with MaMiCo. - * - * When MaMiCo is coupled with any simulation software, it needs to control the simulation from both the outside and inside. - * From the outside, MaMiCo needs to be able to start and stop the simulation, and keep unique simulations distinct. - * From the inside, MaMiCo code needs to be executed at specific points within the same simulation step. This plugin enables the "inside" behaviour. - * - * MaMiCo coupling requires the MAMICO_COUPLING flag to be set, and the MAMICO_SRD_DIR flag to point to MaMiCo source files. With these enabled, ls1 compiles as a library. - * To make sure the program compiles and works with tests, the relevant MaMiCo portions for the code are put in #ifdef regions. - * - * */ +#pragma once -class MamicoCoupling: public PluginBase { +#include +#include +#include "PluginBase.h" +/** + * @brief Allows execution of MaMiCo code to enable coupling with MaMiCo. + * @author Amartya Das Sharma + * + * When MaMiCo is coupled with any simulation software, it needs to control the simulation from both the outside and + * inside. From the outside, MaMiCo needs to be able to start and stop the simulation, and keep unique simulations + * distinct. From the inside, MaMiCo code needs to be executed at specific points within the same simulation step. This + * plugin enables the "inside" behaviour. + * + * MaMiCo coupling requires the MAMICO_COUPLING flag to be set, and the MAMICO_SRD_DIR flag to point to MaMiCo source + * files. With these enabled, ls1 compiles as a library. To make sure the program compiles and works with tests, the + * relevant MaMiCo portions for the code are put in #ifdef regions. + * + * \code{.xml} + * + * + * \endcode + */ +class MamicoCoupling : public PluginBase { public: - MamicoCoupling() {} - virtual ~MamicoCoupling() {} - - /** @brief Prints to log that mamico coupling is initialized - * - * */ - void init(ParticleContainer* particleContainer, - DomainDecompBase* domainDecomp, Domain* domain) override; + MamicoCoupling() = default; + ~MamicoCoupling() override = default; + void init(ParticleContainer *particleContainer, DomainDecompBase *domainDecomp, Domain *domain) override; - void readXML(XMLfileUnits& xmlconfig) override; + /** + * No XML tags defined for this plugin, so does nothing + */ + void readXML(XMLfileUnits &) {} - void beforeEventNewTimestep( - ParticleContainer* particleContainer, DomainDecompBase* domainDecomp, - unsigned long simstep - ) override; + void beforeEventNewTimestep(ParticleContainer *, DomainDecompBase *, unsigned long) {} - /** @brief Takes coupling steps such as particle insertion, to make sure they are accounted for before forces are calculated. - * + /** + * @brief Takes coupling steps such as particle insertion, to make sure they are + * accounted for before forces are calculated. + * * Following steps are taken, if coupling is switched on: * - Iterate over cells to average values like momentum and mass, to pass to macroscopic solvers * - Distribute incoming mass from macroscopic solver by inserting perticles (if enabled) * - Run the MaMiCo thermostat cell by cell - * - * The distributeMass method calls the updateParticleContainerAndDecomposition() function at the end, so we end up with a container with halo particles present. - * Hence a manual halo clearance is done to restore the state of the container. - * */ - void beforeForces( - ParticleContainer* particleContainer, DomainDecompBase* domainDecomp, - unsigned long simstep - ) override; - - /** @brief Performs adjustments after force calculation - * + * + * The distributeMass method calls the updateParticleContainerAndDecomposition() function at the end, so we end up + * with a container with halo particles present. Hence a manual halo clearance is done to restore the state of the + * container. + */ + void beforeForces(ParticleContainer *particleContainer, DomainDecompBase *domainDecomp, + unsigned long simstep) override; + + /** + * @brief Performs adjustments after force calculation + * * Following steps are taken, if coupling is switched on: * - Distribute incoming momentum among affected cells * - Apply boundary force to molecules near microscopic domain boundary - * */ - void afterForces( - ParticleContainer* particleContainer, DomainDecompBase* domainDecomp, - unsigned long simstep - ) override; - - void endStep( - ParticleContainer* particleContainer, DomainDecompBase* domainDecomp, - Domain* domain, unsigned long simstep) override; - - void finish(ParticleContainer* particleContainer, - DomainDecompBase* domainDecomp, Domain* domain) override; - - std::string getPluginName() override { - return std::string("MamicoCoupling"); - } - - static PluginBase* createInstance() { return new MamicoCoupling(); } - -#ifdef MAMICO_COUPLING -/** @brief Sets the macroscopicCellService object that controls the inner coupling logic. - * - * MaMiCo extracts the MamicoCoupling plugin object from the simulation object after initialization and uses this function to set the macroscopicCellCervice. - * The code for this object can be found in https://github.com/HSU-HPC/MaMiCo/blob/master/coupling/services/MacroscopicCellService.cpph - * */ - void setMamicoMacroscopicCellService(coupling::services::MacroscopicCellService<3>* macroscopicCellService){ - _macroscopicCellService = static_cast*> - (macroscopicCellService); + */ + void afterForces(ParticleContainer *particleContainer, DomainDecompBase *domainDecomp, + unsigned long simstep) override; + + void endStep(ParticleContainer *, DomainDecompBase *, Domain *, unsigned long) {} + + void finish(ParticleContainer *, DomainDecompBase *, Domain *) {} + + std::string getPluginName() override { return "MamicoCoupling"; } + + static PluginBase *createInstance() { return new MamicoCoupling(); } + + /** + * @brief Sets the macroscopicCellService object that controls the inner coupling logic. + * + * MaMiCo extracts the MamicoCoupling plugin object from the simulation object after initialization and uses this + * function to set the macroscopicCellCervice. + * + * The code for this object can be found in + * https://github.com/HSU-HPC/MaMiCo/blob/master/coupling/services/CouplingCellService.cpph + */ + void setMamicoCouplingCellService(coupling::services::CouplingCellService<3> *couplingCellService) { + _couplingCellService = + static_cast *>(couplingCellService); } -#endif - -/** @brief Enables coupling logic, allowing coupling steps to run while simulation is running. - * - * Set from within MaMiCo, check https://github.com/HSU-HPC/MaMiCo/blob/master/coupling/interface/MDSimulationFactory.h, class LS1MDSimulation - * */ - void switchOnCoupling(){ _couplingEnabled = true; } - -/** @brief Disables coupling logic, not allowing coupling steps to run. Typically done when equilibrating the simulation initially. - * - * Set from within MaMiCo, check https://github.com/HSU-HPC/MaMiCo/blob/master/coupling/interface/MDSimulationFactory.h, class LS1MDSimulation - * */ - void switchOffCoupling(){ _couplingEnabled = false; } - -/** @brief Getter method for the coupling state. - * - * Not used currently, but may be used in future - * */ - bool getCouplingState() { return _couplingEnabled;} + /** + * @brief Enables coupling logic, allowing coupling steps to run while simulation is running. + * + * Set from within MaMiCo, check + * https://github.com/HSU-HPC/MaMiCo/blob/master/coupling/interface/MDSimulationFactory.h, class LS1MDSimulation + */ + void switchOnCoupling() { _couplingEnabled = true; } + + /** + * @brief Disables coupling logic, not allowing coupling steps to run. Typically done + * when equilibrating the simulation initially. + * + * Set from within MaMiCo, check + * https://github.com/HSU-HPC/MaMiCo/blob/master/coupling/interface/MDSimulationFactory.h, + * class LS1MDSimulation + */ + void switchOffCoupling() { _couplingEnabled = false; } + + /** + * @brief Getter method for the coupling state. + * + * Not used currently, but may be used in future + */ + bool getCouplingState() const { return _couplingEnabled; } private: -#ifdef MAMICO_COUPLING - coupling::services::MacroscopicCellServiceImpl* _macroscopicCellService; -#endif + coupling::services::CouplingCellServiceImpl *_couplingCellService = nullptr; bool _couplingEnabled = false; }; + +#endif // MAMICO_COUPLING diff --git a/src/steereoCommands/sendCouplingMDCommand.cpp b/src/steereoCommands/sendCouplingMDCommand.cpp index 1ee3ca5753..60966e40c7 100644 --- a/src/steereoCommands/sendCouplingMDCommand.cpp +++ b/src/steereoCommands/sendCouplingMDCommand.cpp @@ -97,7 +97,7 @@ ReturnType SendCouplingMDCommand::executeProcessing() double rmax = moleculeContainer->getBoundingBoxMax(dim); logger->debug() << "dim is " << dim << ", dir is " << dir << std::endl; - logger->debug() << "halo is " << moleculeContainer->get_halo_L(dim) << std::endl; + logger->debug() << "halo is " << moleculeContainer->getHaloWidthForDimension(dim) << std::endl; Molecule* currentMolecule; double low_limit = rmin; // particles below this limit have to be copied or moved to the lower process diff --git a/src/utils/Math.h b/src/utils/Math.h index b1a9a1c23d..7ba96efe5e 100644 --- a/src/utils/Math.h +++ b/src/utils/Math.h @@ -10,6 +10,7 @@ #include #include +#include template struct square @@ -50,5 +51,31 @@ inline void LinearRegression(const std::vector& x, const std::vector& y, beta2 = beta1 = 0.; } +/** + * @brief Check whether two floats are within some maximum relative difference. + * Used to check for equality while taking floating-point errors into account. + * + * Taken from AutoPas - src/autopas/utils/Math.h + * + * @param a + * @param b + * @param maxRelativeDifference optional, default value 1e-9 + * @return true + * @return false + */ +bool inline isNearRel(double a, double b, double maxRelativeDifference = 1e-9) +{ + const auto greaterNumber = std::max(std::abs(a), std::abs(b)); + const auto absoluteDifference = maxRelativeDifference * greaterNumber; + const auto diff = std::abs(a - b); + return diff <= absoluteDifference; +} +/** + * @brief Find whether an int is positive or negative (zero is positive). + * + * @param n the number to be checked + * @return int -1 if the number is negative, and 1 otherwise + */ +inline short int findSign(int n) { return n < 0 ? -1 : 1; } #endif /* MATH_H_ */ diff --git a/src/utils/generator/ReplicaFiller.cpp b/src/utils/generator/ReplicaFiller.cpp index f2fc5d0fe6..eb7d054943 100644 --- a/src/utils/generator/ReplicaFiller.cpp +++ b/src/utils/generator/ReplicaFiller.cpp @@ -95,7 +95,7 @@ class ParticleContainerToBasisWrapper : public ParticleContainer { void deleteOuterParticles() override {} - double get_halo_L(int index) const override { return 0.0; } + double getHaloWidthForDimension(int index) const override { return 0.0; } double getCutoff() const override { return 0.0; }