Skip to content

Commit

Permalink
Merge branch 'master' into reflective-boundaries
Browse files Browse the repository at this point in the history
  • Loading branch information
amartyads committed Sep 10, 2024
2 parents 6f7632a + f347b7f commit 552d79f
Show file tree
Hide file tree
Showing 4 changed files with 103 additions and 72 deletions.
30 changes: 15 additions & 15 deletions src/Domain.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,7 @@ Domain::Domain(int rank) {
this->_universalSelectiveThermostatWarning = 0;
this->_universalSelectiveThermostatError = 0;

// explosion heuristics, NOTE: turn off when using slab thermostat
// explosion heuristics, NOTE: turn off when using slab thermostat
_bDoExplosionHeuristics = true;
}

Expand Down Expand Up @@ -220,33 +220,33 @@ void Domain::calculateGlobalValues(
domainDecomp->collCommAppendUnsLong(numMolecules);
domainDecomp->collCommAppendUnsLong(rotDOF);
domainDecomp->collCommAllreduceSumAllowPrevious();
summv2 = domainDecomp->collCommGetDouble();
sumIw2 = domainDecomp->collCommGetDouble();
_globalsummv2 = domainDecomp->collCommGetDouble();
_globalsumIw2 = domainDecomp->collCommGetDouble();
numMolecules = domainDecomp->collCommGetUnsLong();
rotDOF = domainDecomp->collCommGetUnsLong();
domainDecomp->collCommFinalize();
Log::global_log->debug() << "[ thermostat ID " << thermit->first << "]\tN = " << numMolecules << "\trotDOF = " << rotDOF
<< "\tmv2 = " << summv2 << "\tIw2 = " << sumIw2 << std::endl;
<< "\tmv2 = " << _globalsummv2 << "\tIw2 = " << _globalsumIw2 << std::endl;

this->_universalThermostatN[thermit->first] = numMolecules;
this->_universalRotationalDOF[thermit->first] = rotDOF;
mardyn_assert((summv2 > 0.0) || (numMolecules == 0));
mardyn_assert((_globalsummv2 > 0.0) || (numMolecules == 0));

/* calculate the temperature of the entire system */
if(numMolecules > 0)
_globalTemperatureMap[thermit->first] =
(summv2 + sumIw2) / (double)(3*numMolecules + rotDOF);
(_globalsummv2 + _globalsumIw2) / (double)(3*numMolecules + rotDOF);
else
_globalTemperatureMap[thermit->first] = _universalTargetTemperature[thermit->first];

double Ti = Tfactor * _universalTargetTemperature[thermit->first];
if((Ti > 0.0) && (numMolecules > 0) && !_universalNVE)
{
_universalBTrans[thermit->first] = pow(3.0*numMolecules*Ti / summv2, 0.4);
if( sumIw2 == 0.0 )
_universalBTrans[thermit->first] = pow(3.0*numMolecules*Ti / _globalsummv2, 0.4);
if( _globalsumIw2 == 0.0 )
_universalBRot[thermit->first] = 1.0;
else
_universalBRot[thermit->first] = pow(rotDOF*Ti / sumIw2, 0.4);
_universalBRot[thermit->first] = pow(rotDOF*Ti / _globalsumIw2, 0.4);
}
else
{
Expand Down Expand Up @@ -519,7 +519,7 @@ void Domain::writeCheckpointHeader(std::string filename,
#ifndef NDEBUG
checkpointfilestream << "# rho\t" << this->_globalRho << "\n";
//checkpointfilestream << "# rc\t" << global_simulation->getcutoffRadius() << "\n";
checkpointfilestream << "# \n# Please address your questions and suggestions to\n# the ls1 mardyn contact point: <[email protected]>.\n# \n";
checkpointfilestream << "# \n# Please address your questions and suggestions to\n# the ls1 mardyn contact point: <[email protected]>.\n# \n";
#endif
/* by Stefan Becker: the output line "I ..." causes an error: the restart run does not start!!!
if(this->_globalUSteps > 1)
Expand Down Expand Up @@ -872,10 +872,10 @@ void Domain::setNumFluidComponents(unsigned nc){_numFluidComponent = nc;}
unsigned Domain::getNumFluidComponents(){return _numFluidComponent;}

unsigned long Domain::getNumFluidMolecules(){
unsigned long numFluidMolecules = 0;
for(unsigned i = 0; i < _numFluidComponent; i++){
Component& ci=*(global_simulation->getEnsemble()->getComponent(i));
numFluidMolecules+=ci.getNumMolecules();
}
unsigned long numFluidMolecules = 0;
for(unsigned i = 0; i < _numFluidComponent; i++){
Component& ci=*(global_simulation->getEnsemble()->getComponent(i));
numFluidMolecules+=ci.getNumMolecules();
}
return numFluidMolecules;
}
33 changes: 22 additions & 11 deletions src/Domain.h
Original file line number Diff line number Diff line change
Expand Up @@ -215,11 +215,18 @@ class Domain {
//! Before this method is called, it has to be sure that the
//! global potential has been calculated (method calculateGlobalValues)
double getAverageGlobalUpot();
double getGlobalUpot() 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();

//! @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; }

//! 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();
Expand Down Expand Up @@ -365,8 +372,8 @@ class Domain {
void evaluateRho(unsigned long localN, DomainDecompBase* comm);
void submitDU(unsigned cid, double DU, double* r);
void setLambda(double lambda) { this->_universalLambda = lambda; }
void setDensityCoefficient(float coeff) { _globalDecisiveDensity = coeff; }
void setProfiledComponentMass(double m) { _universalProfiledComponentMass = m; }
void setDensityCoefficient(float coeff) { _globalDecisiveDensity = coeff; }
void setProfiledComponentMass(double m) { _universalProfiledComponentMass = m; }

void init_cv(unsigned N, double U, double UU) {
this->_globalUSteps = N;
Expand All @@ -376,7 +383,7 @@ class Domain {
void record_cv();
double cv();

// by Stefan Becker <[email protected]>
// 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);
Expand All @@ -386,8 +393,8 @@ class Domain {
void setUpotCorr(double upotcorr){ _UpotCorr = upotcorr; }
void setVirialCorr(double virialcorr){ _VirialCorr = virialcorr; }

// explosion heuristics, NOTE: turn off when using slab thermostat
void setExplosionHeuristics(bool bVal) { _bDoExplosionHeuristics = bVal; }
// explosion heuristics, NOTE: turn off when using slab thermostat
void setExplosionHeuristics(bool bVal) { _bDoExplosionHeuristics = bVal; }

private:

Expand All @@ -408,6 +415,10 @@ class Domain {
double _globalUpot;
//! global component specific potential (fluid-fluid and fluid-solid but NOT solid-solid)
double _globalUpotCspecif;
//! global translational kinetic energy times two
double _globalsummv2;
//! global rotational kinetic energy times two
double _globalsumIw2;
//! global virial
double _globalVirial;
//! global density
Expand Down Expand Up @@ -455,9 +466,9 @@ class Domain {
double _globalSigmaUU;
//! which components should be considered?
std::map<unsigned, bool> _universalProfiledComponents;
double _universalProfiledComponentMass; // set from outside
double _universalLambda; // set from outside
float _globalDecisiveDensity; // set from outside
double _universalProfiledComponentMass; // set from outside
double _universalLambda; // set from outside
float _globalDecisiveDensity; // set from outside

int _universalSelectiveThermostatCounter;
int _universalSelectiveThermostatWarning;
Expand Down Expand Up @@ -486,8 +497,8 @@ class Domain {
//! parameter streams for each possible pair of molecule-types
Comp2Param _comp2params;

// explosion heuristics, NOTE: turn off when using slab thermostat
bool _bDoExplosionHeuristics;
// explosion heuristics, NOTE: turn off when using slab thermostat
bool _bDoExplosionHeuristics;
};


Expand Down
91 changes: 58 additions & 33 deletions src/io/ResultWriter.cpp
Original file line number Diff line number Diff line change
@@ -1,53 +1,63 @@
#include "io/ResultWriter.h"

#include <chrono>
#include <fstream>

#include "Domain.h"
#include "parallel/DomainDecompBase.h"
#include "Simulation.h"
#include "utils/Logger.h"
#include <chrono>
#include "utils/mardyn_assert.h"


void ResultWriter::readXML(XMLfileUnits& xmlconfig) {
_writeFrequency = 1;
xmlconfig.getNodeValue("writefrequency", _writeFrequency);
Log::global_log->info() << "[ResultWriter] Write frequency: " << _writeFrequency << std::endl;
if (_writeFrequency <= 0) {
Log::global_log->error() << "[ResultWriter] Write frequency must be a positive nonzero integer, but is " << _writeFrequency << std::endl;
mardyn_exit(123);
}

_outputPrefix = "mardyn";
xmlconfig.getNodeValue("outputprefix", _outputPrefix);
Log::global_log->info() << "[ResultWriter] Output prefix: " << _outputPrefix << std::endl;

size_t acc_steps = 1000;
xmlconfig.getNodeValue("accumulation_steps", acc_steps);
_U_pot_acc = new Accumulator<double>(acc_steps);
_p_acc = new Accumulator<double>(acc_steps);
_U_pot_acc = std::make_unique<Accumulator<double>>(acc_steps);
_U_kin_acc = std::make_unique<Accumulator<double>>(acc_steps);
_p_acc = std::make_unique<Accumulator<double>>(acc_steps);
Log::global_log->info() << "[ResultWriter] Accumulation steps: " << acc_steps << std::endl;

_writePrecision = 5;
xmlconfig.getNodeValue("writeprecision", _writePrecision);
Log::global_log->info() << "[ResultWriter] Write precision: " << _writePrecision << std::endl;
_writeWidth = _writePrecision + 15; // Adding a width of 15 to have enough whitespace between columns
}

void ResultWriter::init(ParticleContainer * /*particleContainer*/,
DomainDecompBase *domainDecomp, Domain * /*domain*/) {

if(domainDecomp->getRank() == 0) {
const std::string resultfile(_outputPrefix+".res");
std::ofstream resultStream;
resultStream.open(resultfile.c_str(), std::ios::out);
const auto now = std::chrono::system_clock::to_time_t(std::chrono::system_clock::now());
tm unused{};
const auto nowStr = std::put_time(localtime_r(&now, &unused), "%c");
std::string resultfile(_outputPrefix+".res");
_resultStream.open(resultfile.c_str(), std::ios::out);
_resultStream << "# ls1 MarDyn simulation started at " << nowStr << std::endl;
_resultStream << "# Averages are the accumulated values over " << _U_pot_acc->getWindowLength() << " time steps."<< std::endl;
_resultStream << std::setw(10) << "# step" << std::setw(_writePrecision+15) << "time"
<< std::setw(_writePrecision+15) << "U_pot"
<< std::setw(_writePrecision+15) << "U_pot_avg"
<< std::setw(_writePrecision+15) << "p"
<< std::setw(_writePrecision+15) << "p_avg"
<< std::setw(_writePrecision+15) << "beta_trans"
<< std::setw(_writePrecision+15) << "beta_rot"
<< std::setw(_writePrecision+15) << "c_v"
<< std::setw(_writePrecision+15) << "N"
resultStream << "# ls1 MarDyn simulation started at " << nowStr << std::endl;
resultStream << "# Averages are the accumulated values over " << _U_pot_acc->getWindowLength() << " time steps."<< std::endl;
resultStream << std::setw(10) << "# step" << std::setw(_writeWidth) << "time"
<< std::setw(_writeWidth) << "U_pot"
<< std::setw(_writeWidth) << "U_pot_avg"
<< std::setw(_writeWidth) << "U_kin"
<< std::setw(_writeWidth) << "U_kin_avg"
<< std::setw(_writeWidth) << "p"
<< std::setw(_writeWidth) << "p_avg"
<< std::setw(_writeWidth) << "beta_trans"
<< std::setw(_writeWidth) << "beta_rot"
<< std::setw(_writeWidth) << "c_v"
<< std::setw(_writeWidth) << "N"
<< std::endl;
resultStream.close();
}
}

Expand All @@ -58,31 +68,46 @@ void ResultWriter::endStep(ParticleContainer *particleContainer, DomainDecompBas

unsigned long globalNumMolecules = domain->getglobalNumMolecules(true, particleContainer, domainDecomp);
double cv = domain->cv();
double ekin = domain->getGlobalUkinTrans()+domain->getGlobalUkinRot();

_U_pot_acc->addEntry(domain->getGlobalUpot());
_U_kin_acc->addEntry(ekin);
_p_acc->addEntry(domain->getGlobalPressure());
if((domainDecomp->getRank() == 0) && (simstep % _writeFrequency == 0)){
_resultStream << std::setw(10) << simstep << std::setw(_writePrecision+15) << std::scientific << std::setprecision(_writePrecision) << _simulation.getSimulationTime()
<< std::setw(_writePrecision+15) << std::scientific << std::setprecision(_writePrecision) << domain->getGlobalUpot()
<< std::setw(_writePrecision+15) << std::scientific << std::setprecision(_writePrecision) << _U_pot_acc->getAverage()
<< std::setw(_writePrecision+15) << std::scientific << std::setprecision(_writePrecision) << domain->getGlobalPressure()
<< std::setw(_writePrecision+15) << std::scientific << std::setprecision(_writePrecision) << _p_acc->getAverage()
<< std::setw(_writePrecision+15) << std::scientific << std::setprecision(_writePrecision) << domain->getGlobalBetaTrans()
<< std::setw(_writePrecision+15) << std::scientific << std::setprecision(_writePrecision) << domain->getGlobalBetaRot()
<< std::setw(_writePrecision+15) << std::scientific << std::setprecision(_writePrecision) << cv
<< std::setw(_writePrecision+15) << std::scientific << std::setprecision(_writePrecision) << globalNumMolecules
<< std::endl;
if ((domainDecomp->getRank() == 0) && (simstep % _writeFrequency == 0)){
const std::string resultfile(_outputPrefix+".res");
std::ofstream resultStream;
resultStream.open(resultfile.c_str(), std::ios::app);
auto printOutput = [&](auto value) {
resultStream << std::setw(_writeWidth) << std::scientific << std::setprecision(_writePrecision) << value;
};
resultStream << std::setw(10) << simstep;
printOutput(_simulation.getSimulationTime());
printOutput(domain->getGlobalUpot());
printOutput(_U_pot_acc->getAverage());
printOutput(ekin);
printOutput(_U_kin_acc->getAverage());
printOutput(domain->getGlobalPressure());
printOutput(_p_acc->getAverage());
printOutput(domain->getGlobalBetaTrans());
printOutput(domain->getGlobalBetaRot());
printOutput(cv);
printOutput(globalNumMolecules);
resultStream << std::endl;
resultStream.close();
}
}

void ResultWriter::finish(ParticleContainer * /*particleContainer*/,
DomainDecompBase *domainDecomp, Domain * /*domain*/){

if(domainDecomp->getRank() == 0) {
if (domainDecomp->getRank() == 0) {
const std::string resultfile(_outputPrefix+".res");
std::ofstream resultStream;
resultStream.open(resultfile.c_str(), std::ios::app);
const auto now = std::chrono::system_clock::to_time_t(std::chrono::system_clock::now());
tm unused{};
const auto nowStr = std::put_time(localtime_r(&now, &unused), "%c");
_resultStream << "# ls1 mardyn simulation finished at " << nowStr << std::endl;
_resultStream.close();
resultStream << "# ls1 mardyn simulation finished at " << nowStr << std::endl;
resultStream.close();
}
}
21 changes: 8 additions & 13 deletions src/io/ResultWriter.h
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
#ifndef SRC_IO_RESULTWRITER_H_
#define SRC_IO_RESULTWRITER_H_

#include <fstream>
#include <memory>
#include <string>

#include "plugins/PluginBase.h"
Expand All @@ -22,12 +22,6 @@
*/
class ResultWriter : public PluginBase {
public:
ResultWriter() : _U_pot_acc(nullptr), _p_acc(nullptr) {}
~ResultWriter() {
delete _U_pot_acc;
delete _p_acc;
}

/** @brief Read in XML configuration for ResultWriter and all its included objects.
*
* The following xml object structure is handled by this method:
Expand Down Expand Up @@ -60,12 +54,13 @@ class ResultWriter : public PluginBase {
static PluginBase* createInstance() { return new ResultWriter(); }

private:
std::ofstream _resultStream;
long _writeFrequency;
int _writePrecision;
std::string _outputPrefix;
Accumulator<double> *_U_pot_acc;
Accumulator<double> *_p_acc;
long _writeFrequency{1000UL};
unsigned int _writePrecision{5};
unsigned int _writeWidth{20};
std::string _outputPrefix{"results"};
std::unique_ptr<Accumulator<double>> _U_pot_acc;
std::unique_ptr<Accumulator<double>> _U_kin_acc;
std::unique_ptr<Accumulator<double>> _p_acc;
};

#endif // SRC_IO_RESULTWRITER_H_

0 comments on commit 552d79f

Please sign in to comment.