Merge pull request #313 from ls1mardyn/fixInitSimstep

Set initial simulation time step earlier
ls1mardyn · Jun 3, 2024 · 4e77640 · 4e77640
2 parents baff96a + 05ccdeb
commit 4e77640
Showing 1 changed file with 10 additions and 4 deletions.
diff --git a/src/Simulation.cpp b/src/Simulation.cpp
@@ -196,7 +196,8 @@ void Simulation::readXML(XMLfileUnits& xmlconfig) {
 	/* run section */
 	if(xmlconfig.changecurrentnode("run")) {
 		xmlconfig.getNodeValueReduced("currenttime", _simulationTime);
-		Log::global_log->info() << "Simulation start time: " << _simulationTime << std::endl;
+		Log::global_log->info() << "Simulation start time (might be altered by checkpoint): "
+								<< getSimulationTime() << std::endl;
 		/* steps */
 		xmlconfig.getNodeValue("equilibration/steps", _initStatistics);
 		Log::global_log->info() << "Number of equilibration steps: " << _initStatistics << std::endl;
@@ -682,11 +683,18 @@ void Simulation::readXML(XMLfileUnits& xmlconfig) {
 		bool ignoreCheckpointTime = false;
 		if(xmlconfig.getNodeValue("ignoreCheckpointTime", ignoreCheckpointTime)) {
 			if(ignoreCheckpointTime)
-				_simulationTime = 0;
+				setSimulationTime(0.0);
 		}
 	}
 
 	xmlconfig.changecurrentnode(oldpath);
+
+	// _simulationTime might be updated by readers -> also update _initSimulation
+	// In theory, the division of the physical time (_simulationTime) and the time step width,
+	// should give a whole number. However, due to numerical errors, the result of the
+	// division must be rounded
+	_simstep = _initSimulation = static_cast<unsigned long>(std::round(_simulationTime / _integrator->getTimestepLength()));
+	Log::global_log->info() << "Set initial time step to " << _initSimulation << std::endl;
 }
 
 
@@ -942,8 +950,6 @@ void Simulation::prepare_start() {
 
 	_ensemble->prepare_start();
 
-	_simstep = _initSimulation = (unsigned long) round(_simulationTime / _integrator->getTimestepLength() );
-	Log::global_log->info() << "Set initial time step to start from to " << _initSimulation << std::endl;
 	Log::global_log->info() << "System initialised with " << _domain->getglobalNumMolecules(true, _moleculeContainer, _domainDecomposition) << " molecules." << std::endl;
 
 }