Fixes to mamico plugin after mamico refactor

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()

examples/all-options.xml

@@ -197,6 +197,16 @@
           <overlappingCollectives>False</overlappingCollectives> <!-- true if overlapping collectives should be used, false otherwise. REQUIRES MPI>=3! -->
           <overlappingStartAtStep>5</overlappingStartAtStep> <!-- Start overlapping at given step (default: 5), only relevant if overlappingCollectives==True -->
           <overlappingP2P>False</overlappingP2P> <!-- Defines whether to use overlapping p2p communication or not. Default: False -->
+          <!-- Select the boundary type for each dimension. Available options are reflecting/reflective, outflow and periodic (default).
+            particles interacting with a periodic boundary for a certain axis are copied over to the corresponding subdomain, as usual
+            particles interacting with a reflecting boundary for a certain axis have their velocities reversed for that axis
+            particles interacting with an outflow boundary with any axis get deleted, irrespective of the other boundary conditions
+            -->
+          <boundaries> 
+              <x>reflective</x>
+              <y>outflow</y>
+              <z>periodic</z>
+          </boundaries> 
       </parallelisation>
       <parallelisation type="StaticIrregDomainDecomposition">
           <!-- Provide weights that determine the relative widths of each subdomain

examples/simple-boundary-test/ls1config.xml

@@ -0,0 +1,76 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<mardyn version="20100525">
+	<refunits type="SI">
+		<length unit="nm">0.34</length>
+		<mass unit="u">39.948</mass>
+		<energy unit="K">120</energy>
+	</refunits>
+	<simulation type="MD">
+		<integrator type="Leapfrog">
+			<timestep unit="reduced">0.005</timestep>
+		</integrator>
+		<run>
+			<currenttime>0</currenttime>
+			<equilibration><steps>0</steps></equilibration>
+			<production>
+				<steps>800</steps>
+			</production>
+		</run>
+		<ensemble type="NVT">
+			<temperature unit="reduced">1.1</temperature>
+			<domain type="box">
+				<lx>10</lx>
+				<ly>10</ly>
+				<lz>10</lz>
+			</domain>
+			<components>
+				<moleculetype id="1" name="Arbitrary">
+					<site type="LJ126" id="1" >
+						<coords> <x>0.0</x> <y>0.0</y> <z>0.0</z> </coords>
+						<mass>1.0</mass>
+						<sigma>1.0</sigma>
+						<epsilon>1.0</epsilon>
+						<shifted>true</shifted>
+					</site>
+				</moleculetype>
+			</components>
+			<phasespacepoint>
+				<file type="ASCII">simple_checkpoint.inp</file>
+			</phasespacepoint>
+		</ensemble>
+		<algorithm>
+			<parallelisation type="DomainDecomposition">
+			<!--<updateFrequency> 10000 </updateFrequency>
+			<loadBalancer type="none">
+
+                </loadBalancer>-->
+				<!--<MPIGridDims> <x>2</x><y>2</y> <z>2</z> </MPIGridDims>-->
+				<boundaries> <x>reflective</x><y>outflow</y><z>reflective</z> </boundaries>
+			</parallelisation>
+			<datastructure type="AutoPas">
+			</datastructure>
+			<cutoffs type="CenterOfMass">
+				<radiusLJ unit="reduced">2.2</radiusLJ>
+			</cutoffs>
+			<electrostatic type="ReactionField">
+				<epsilon>1.0e+10</epsilon>
+			</electrostatic>
+		</algorithm>
+		<output>
+			<outputplugin name="VTKMoleculeWriter">
+				<outputprefix>vtkOutput</outputprefix>
+				<writefrequency>1</writefrequency>
+			</ouputplugin>
+			<outputplugin name="Adios2Writer">
+		  		<outputfile>particles.bp</outputfile>
+		  		<adios2enginetype>BP4</adios2enginetype>
+		  		<writefrequency>1</writefrequency>
+			</outputplugin>
+			<!--<outputplugin name="CheckpointWriter">
+        		<type>ASCII</type>
+        		<writefrequency>1</writefrequency>
+        		<outputprefix>cp_binary</outputprefix>
+			</outputplugin>-->
+		</output>
+	</simulation>
+</mardyn>

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

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}

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 <[email protected]>
 // 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));

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 <[email protected]>
 	/* 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; }

src/Simulation.cpp

@@ -86,6 +86,8 @@
 #include <coupling/interface/impl/ls1/LS1StaticCommData.h>
 #endif
 
+#include "parallel/boundaries/BoundaryUtils.h"
+
 Simulation* global_simulation;
 
 Simulation::Simulation()
@@ -424,6 +426,34 @@ void Simulation::readXML(XMLfileUnits& xmlconfig) {
 			}
 			_lastTraversalTimeHistory.setCapacity(timerForLoadAveragingLength);
 
+			if(xmlconfig.changecurrentnode("boundaries")) {
+				std::string xBoundaryFromConfig, yBoundaryFromConfig, zBoundaryFromConfig;
+				xmlconfig.getNodeValue("x", xBoundaryFromConfig);
+				xmlconfig.getNodeValue("y", yBoundaryFromConfig);
+				xmlconfig.getNodeValue("z", zBoundaryFromConfig);
+				BoundaryUtils::BoundaryType xBoundary = BoundaryUtils::convertStringToBoundary(xBoundaryFromConfig);
+				BoundaryUtils::BoundaryType yBoundary = BoundaryUtils::convertStringToBoundary(yBoundaryFromConfig);
+				BoundaryUtils::BoundaryType zBoundary = BoundaryUtils::convertStringToBoundary(zBoundaryFromConfig);
+				_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()) {
+					MARDYN_EXIT("Invalid boundary type! Please check the config file");
+				}
+				if(_overlappingP2P && _domainDecomposition->hasGlobalNonPeriodicBoundary()) {
+					MARDYN_EXIT("Non-periodic boundaries not supported with overlappingP2P enabled! Exiting...");
+				}
+				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->setBoundsAndGlobalBoundaries(_domain, _ensemble);
+				xmlconfig.changecurrentnode("..");
+			}
+
 			xmlconfig.changecurrentnode("..");
 		}
 		else {
@@ -1014,13 +1044,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 +1108,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 +1437,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");