Skip to content

Commit

Permalink
Andersen Thermostat integrated in TemperatureControl
Browse files Browse the repository at this point in the history 
git-svn-id: https://scm.projects.hlrs.de/svn/ls1/MarDyn/trunk@6368 a63bd714-7e14-4b5e-94e7-302c8c8ff188
  • Loading branch information
@kruegener
kruegener committed Jul 11, 2018
1 parent 21ea4dc commit 474a4b4
Show file tree
Hide file tree
Showing 4 changed files with 243 additions and 94 deletions.
26 changes: 10 additions & 16 deletions src/Simulation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -401,22 +401,16 @@ void Simulation::readXML(XMLfileUnits& xmlconfig) {
global_log->info() << "Adding velocity scaling thermostat for component '" << componentName << "' (ID: " << componentId << "), T = " << temperature << endl;
}
}
else if(thermostattype == "TemperatureControl")
{
if(NULL == _temperatureControl)
{
_temperatureControl = new TemperatureControl();
_temperatureControl->readXML(xmlconfig);
}
else
{
global_log->error() << "Instance of TemperatureControl allready exist! Programm exit ..." << endl;
Simulation::exit(-1);
}
}
else if(thermostattype == "Andersen"){

}
else if(thermostattype == "TemperatureControl") {
if (NULL == _temperatureControl) {
_temperatureControl = new TemperatureControl();
_temperatureControl->readXML(xmlconfig);
} else {
global_log->error() << "Instance of TemperatureControl allready exist! Programm exit ..."
<< endl;
Simulation::exit(-1);
}
}
else
{
global_log->warning() << "Unknown thermostat " << thermostattype << endl;
Expand Down
282 changes: 205 additions & 77 deletions src/thermostats/TemperatureControl.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,24 +28,24 @@ unsigned short ControlRegionT::_nStaticID = 0;
// class ControlRegionT

ControlRegionT::ControlRegionT()
:
_nID(0),
_dLowerCorner{0.,0.,0.},
_dUpperCorner{0.,0.,0.},
_nNumSlabs(0),
_dSlabWidth(0.0),
_thermVars(),
_dTargetTemperature(0.0),
_dTemperatureExponent(0.0),
_nTargetComponentID(0),
_nNumThermostatedTransDirections(0),
_nRegionID(0),
_accumulator(nullptr),
_strFilenamePrefixBetaLog("beta_log"),
_nWriteFreqBeta(1000),
_numSampledConfigs(0),
_dBetaTransSumGlobal(0.0),
_dBetaRotSumGlobal(0.0)
:
_nID(0),
_dLowerCorner{0.,0.,0.},
_dUpperCorner{0.,0.,0.},
_nNumSlabs(0),
_dSlabWidth(0.0),
_thermVars(),
_dTargetTemperature(0.0),
_dTemperatureExponent(0.0),
_nTargetComponentID(0),
_nNumThermostatedTransDirections(0),
_nRegionID(0),
_accumulator(nullptr),
_strFilenamePrefixBetaLog("beta_log"),
_nWriteFreqBeta(1000),
_numSampledConfigs(0),
_dBetaTransSumGlobal(0.0),
_dBetaRotSumGlobal(0.0)
{
// ID
_nID = ++_nStaticID;
Expand Down Expand Up @@ -134,28 +134,74 @@ void ControlRegionT::readXML(XMLfileUnits& xmlconfig)
xmlconfig.getNodeValue("target/temperature", _dTargetTemperature);
xmlconfig.getNodeValue("target/component", _nTargetComponentID);

// settings
xmlconfig.getNodeValue("settings/numslabs", _nNumSlabs);
xmlconfig.getNodeValue("settings/exponent", _dTemperatureExponent);
xmlconfig.getNodeValue("settings/directions", strDirections);
// calc slab width
_dSlabWidth = this->GetWidth(1) / ( (double)(_nNumSlabs) );
// create accumulator instance
_accumulator = this->CreateAccumulatorInstance(strDirections);

// write control for beta_trans and beta_rot log file
_nWriteFreqBeta = 1000;
_strFilenamePrefixBetaLog = "beta_log";
xmlconfig.getNodeValue("writefreq", _nWriteFreqBeta);
xmlconfig.getNodeValue("fileprefix", _strFilenamePrefixBetaLog);
if(_nWriteFreqBeta==0) {
global_log->warning() << "Temperature Control: write Frequency was specified to be zero. This is NOT allowed. Reset it to 1000." << std::endl;
_nWriteFreqBeta = 1000;
// ControlMethod "Default/Andersen/Mixed"
std::string methods = "";
xmlconfig.getNodeValue("method", methods);
if(methods != "") {
if(methods == "Default") {
_localMethod = Default;

// settings
xmlconfig.getNodeValue("settings/numslabs", _nNumSlabs);
xmlconfig.getNodeValue("settings/exponent", _dTemperatureExponent);
xmlconfig.getNodeValue("settings/directions", strDirections);
// calc slab width
_dSlabWidth = this->GetWidth(1) / ( (double)(_nNumSlabs) );
// create accumulator instance
_accumulator = this->CreateAccumulatorInstance(strDirections);

// write control for beta_trans and beta_rot log file
_nWriteFreqBeta = 1000;
_strFilenamePrefixBetaLog = "beta_log";
xmlconfig.getNodeValue("writefreq", _nWriteFreqBeta);
xmlconfig.getNodeValue("fileprefix", _strFilenamePrefixBetaLog);
if(_nWriteFreqBeta==0) {
global_log->warning() << "Temperature Control: write Frequency was specified to be zero. This is NOT allowed. Reset it to 1000." << std::endl;
_nWriteFreqBeta = 1000;
}
this->InitBetaLogfile();

// init data structures
this->Init();
}
else if(methods == "Andersen") {
_localMethod = Andersen;
xmlconfig.getNodeValue("settings/nu", _nuAndersen);
xmlconfig.getNodeValue("settings/timestep", _timestep);
}
else {
global_log -> error() << "[TemperatureControl] REGION: Invalid 'method' param: " << methods << std::endl;
Simulation::exit(-1);
}
global_log -> info() << "[TemperatureControl] REGION 'method' param: " << methods << std::endl;
}
this->InitBetaLogfile();
else {
_localMethod = Default;
global_log -> info() << "[TemperatureControl] REGION: no method specified, selecting Default" << std::endl;

// settings
xmlconfig.getNodeValue("settings/numslabs", _nNumSlabs);
xmlconfig.getNodeValue("settings/exponent", _dTemperatureExponent);
xmlconfig.getNodeValue("settings/directions", strDirections);
// calc slab width
_dSlabWidth = this->GetWidth(1) / ( (double)(_nNumSlabs) );
// create accumulator instance
_accumulator = this->CreateAccumulatorInstance(strDirections);

// write control for beta_trans and beta_rot log file
_nWriteFreqBeta = 1000;
_strFilenamePrefixBetaLog = "beta_log";
xmlconfig.getNodeValue("writefreq", _nWriteFreqBeta);
xmlconfig.getNodeValue("fileprefix", _strFilenamePrefixBetaLog);
if(_nWriteFreqBeta==0) {
global_log->warning() << "Temperature Control: write Frequency was specified to be zero. This is NOT allowed. Reset it to 1000." << std::endl;
_nWriteFreqBeta = 1000;
}
this->InitBetaLogfile();

// init data structures
this->Init();
// init data structures
this->Init();
}
}

void ControlRegionT::Init()
Expand Down Expand Up @@ -291,36 +337,53 @@ void ControlRegionT::ControlTemperature(Molecule* mol)
return;
}

unsigned int nPosIndex;
unsigned int nIndexMax = _nNumSlabs - 1;
if(_localMethod == Default) {
unsigned int nPosIndex;
unsigned int nIndexMax = _nNumSlabs - 1;

// calc position index
double* dLowerCorner = this->GetLowerCorner();
double dPosRelative = mol->r(1) - dLowerCorner[1];
// calc position index
double *dLowerCorner = this->GetLowerCorner();
double dPosRelative = mol->r(1) - dLowerCorner[1];

nPosIndex = (unsigned int) floor(dPosRelative / _dSlabWidth);
nPosIndex = (unsigned int) floor(dPosRelative / _dSlabWidth);

// ignore outer (halo) molecules
if(nPosIndex > nIndexMax) // negative values will be ignored to: cast to unsigned int --> high value
return;
// ignore outer (halo) molecules
if (nPosIndex > nIndexMax) // negative values will be ignored to: cast to unsigned int --> high value
return;

GlobalThermostatVariables & globalTV = _thermVars[nPosIndex]._global; // do not forget &
if(globalTV._numMolecules < 1)
return;
GlobalThermostatVariables &globalTV = _thermVars[nPosIndex]._global; // do not forget &
if (globalTV._numMolecules < 1)
return;

// scale velocity
double vcorr = 2. - 1. / globalTV._betaTrans;
double Dcorr = 2. - 1. / globalTV._betaRot;
// scale velocity
double vcorr = 2. - 1. / globalTV._betaTrans;
double Dcorr = 2. - 1. / globalTV._betaRot;

/*
mol->setv(0, mol->v(0) * vcorr);
// mol->setv(1, mol->v(1) * vcorr);
mol->setv(2, mol->v(2) * vcorr);
*/

_accumulator->ScaleVelocityComponents(mol, vcorr);
_accumulator->ScaleVelocityComponents(mol, vcorr);

mol->scale_D(Dcorr);
mol->scale_D(Dcorr);
}
else if(_localMethod == Andersen) {
double stdDevTrans, stdDevRot;
if(_rand.rnd() < _nuAndersen*_timestep) {
stdDevTrans = sqrt(_dTargetTemperature/mol->mass());
for(unsigned short d = 0; d < 3; d++) {
stdDevRot = sqrt(_dTargetTemperature*mol->getI(d));
mol->setv(d, _rand.gaussDeviate(stdDevTrans));
mol->setD(d, _rand.gaussDeviate(stdDevRot));
}
}
}
else {
global_log -> error() << "[TemperatureControl] Invalid localMethod param: " << _localMethod << std::endl;
Simulation::exit(-1);
}
}

void ControlRegionT::ResetLocalValues()
Expand All @@ -334,23 +397,25 @@ void ControlRegionT::ResetLocalValues()

void ControlRegionT::InitBetaLogfile()
{
DomainDecompBase* domainDecomp = &(global_simulation->domainDecomposition() );
if(_localMethod == Default) {
DomainDecompBase *domainDecomp = &(global_simulation->domainDecomposition());

#ifdef ENABLE_MPI
int rank = domainDecomp->getRank();
// int numprocs = domainDecomp->getNumProcs();
if (rank!= 0)
return;
int rank = domainDecomp->getRank();
// int numprocs = domainDecomp->getNumProcs();
if (rank!= 0)
return;
#endif

std::stringstream filenamestream;
filenamestream << _strFilenamePrefixBetaLog << "_reg" << this->GetID() << ".dat";
std::stringstream outputstream;
outputstream.write(reinterpret_cast<const char*>(&_nWriteFreqBeta), 8);
std::stringstream filenamestream;
filenamestream << _strFilenamePrefixBetaLog << "_reg" << this->GetID() << ".dat";
std::stringstream outputstream;
outputstream.write(reinterpret_cast<const char *>(&_nWriteFreqBeta), 8);

ofstream fileout(filenamestream.str().c_str(), std::ios::out | std::ios::binary);
fileout << outputstream.str();
fileout.close();
ofstream fileout(filenamestream.str().c_str(), std::ios::out | std::ios::binary);
fileout << outputstream.str();
fileout.close();
}
}

void ControlRegionT::WriteBetaLogfile(unsigned long simstep)
Expand Down Expand Up @@ -406,6 +471,29 @@ void TemperatureControl::readXML(XMLfileUnits& xmlconfig)
global_log->info() << "Control with frequency: " << _nControlFreq << endl;
global_log->info() << "Stop control at simstep: " << _nStop << endl;

// ControlMethod "Default/Andersen/Mixed"
std::string methods = "";
xmlconfig.getNodeValue("control/methods", methods);
if(methods != "") {
if(methods == "Default") {
_method = Default;
}
else if(methods == "Andersen") {
_method = Andersen;
}
else if(methods == "Mixed") {
_method = Mixed;
}
else {
global_log -> error() << "[TemperatureControl] Invalid 'methods' param: " << methods << std::endl;
Simulation::exit(-1);
}
}
else {
_method = Default;
global_log -> info() << "[TemperatureControl] no method specified, selecting Default" << std::endl;
}
global_log -> info() << "[TemperatureControl] 'methods' param: " << methods << std::endl;
// turn on/off explosion heuristics
// domain->SetExplosionHeuristics(bUseExplosionHeuristics);

Expand Down Expand Up @@ -482,8 +570,47 @@ void TemperatureControl::WriteBetaLogfiles(unsigned long simstep)
reg->WriteBetaLogfile(simstep);
}

/**
* @brief Decide which ControlMethod to use
*
* @param domainDecomposition
* @param particleContainer
* @param simstep
*/
void TemperatureControl::DoLoopsOverMolecules(DomainDecompBase* domainDecomposition, ParticleContainer* particleContainer, unsigned long simstep)
{
if(_method == Default){
this->DoLoopsDefault(domainDecomposition, particleContainer, simstep);
}
else if(_method == Andersen){
this->DoLoopsAndersen(domainDecomposition, particleContainer, simstep);
}
else if(_method == Mixed) {
this->DoLoopsDefault(domainDecomposition, particleContainer, simstep);
this->DoLoopsAndersen(domainDecomposition, particleContainer, simstep);
}
else{
global_log -> error() << "Unknown TemperatureControlMethod" << _method << std::endl;
Simulation::exit(-1);
}
ParticleIterator tM;
for( tM = particleContainer->iterator();
tM.hasNext();
tM.next())
{
// control temperature
this->ControlTemperature(&(*tM), simstep);
}
}

/**
* @brief Use default control method and Loop over Molecules
*
* @param domainDecomposition
* @param particleContainer
* @param simstep
*/
void TemperatureControl::DoLoopsDefault(DomainDecompBase* domainDecomposition, ParticleContainer* particleContainer, unsigned long simstep) {
// respect start/stop
if(this->GetStart() <= simstep && this->GetStop() > simstep)
{
Expand All @@ -509,18 +636,19 @@ void TemperatureControl::DoLoopsOverMolecules(DomainDecompBase* domainDecomposit

// write beta_trans, beta_rot log-files
this->WriteBetaLogfiles(simstep);

for( tM = particleContainer->iterator();
tM.hasNext();
tM.next())
{
// control temperature
this->ControlTemperature(&(*tM), simstep);

// cout << "id = " << tM->id() << ", (vx,vy,vz) = " << tM->v(0) << ", " << tM->v(1) << ", " << tM->v(2) << endl;
}
}
}

/**
* @brief Use Andersen Thermostat Control and Loop over Molecules
*
* @param domainDecomposition
* @param particleContainer
* @param simstep
*/
void TemperatureControl::DoLoopsAndersen(DomainDecompBase* domainDecomposition, ParticleContainer* particleContainer, unsigned long simstep) {
global_log->info() << "Andersen Thermostat\n";
}



0 comments on commit 474a4b4

Please sign in to comment.