Page MenuHomec4science
Files F82222403

stimulation_langevin.cc
No OneTemporary
Actions

File Metadata

Created
Tue, Sep 10, 06:16

stimulation_langevin.cc
View Options

/**
* @file stimulation_langevin.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Wed Jul 09 21:59:47 2014
*
* @brief This stimulator applies a Langevin thermostat
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
#include "lm_common.hh"
#include "stimulation_langevin.hh"
#include "lib_md.hh"
#include "lib_continuum.hh"
#include "lib_dd.hh"
#include "filter_manager.hh"
#include "math_tools.hh"
#include "reference_manager.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
template <typename _Input>
StimulationLangevin<_Input>::StimulationLangevin(const std::string & name,
ComponentLMInterface & d):
Stimulation<_Input>(name,d),
friction_restitution_work("friction_restitution:"+name,5),
old_friction_per_atom("old_friction_per_atom:"+name,spatial_dimension),
old_restitution_per_atom("old_restitution_per_atom:"+name,spatial_dimension),
friction_per_atom("friction_per_atom:"+name,spatial_dimension),
restitution_per_atom("restitution_per_atom:"+name,spatial_dimension) {
damp = 100;
temperature = 100;
timestep = 1.;
seed = 0;
compute_work_flag = false;
// initialize the compute
friction_restitution_work.setDim(5);
friction_restitution_work.assign(5,0);
friction_restitution_work.name_computed.push_back("temperature");
friction_restitution_work.name_computed.push_back("frictional_work");
friction_restitution_work.name_computed.push_back("restitution_work");
friction_restitution_work.name_computed.push_back("cumulated_frictional_work");
friction_restitution_work.name_computed.push_back("cumulated_restitution_work");
cumulated_frictional_work = 0;
cumulated_restitution_work = 0;
// register for computes
FilterManager::getManager().addObject(&friction_restitution_work,false);
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
StimulationLangevin<Cont>::~StimulationLangevin(){
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
void StimulationLangevin<Cont>::stimulate(Cont & cont, UInt stage)
{
if (compute_work_flag && compute_work_peratom_flag)
this->addForce<true,true>(cont);
else if ((!compute_work_flag) && compute_work_peratom_flag)
this->addForce<false,true>(cont);
else if ((!compute_work_flag) && (!compute_work_peratom_flag))
this->addForce<false,false>(cont);
else if (compute_work_flag && (!compute_work_peratom_flag))
this->addForce<true,false>(cont);
else LM_FATAL("internal error: should not happend");
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
void StimulationLangevin<Cont>::resetOldWorkPerAtom(Cont & cont){
static const UInt Dim = Cont::Dim;
if (old_friction_per_atom.size() == 0){
old_friction_per_atom.resize(cont.nbElem()*Dim);
old_restitution_per_atom.resize(cont.nbElem()*Dim);
if (cont.hasRefManager()){
cont.getRefManager().attachVector(old_friction_per_atom,cont,Dim);
cont.getRefManager().attachVector(old_restitution_per_atom,cont,Dim);
}
}
else {
if (old_friction_per_atom.size() != cont.size()*Dim)
LM_FATAL("inconsistent previous storage of force vector (friction) "
<< old_friction_per_atom.size() << "(" << &old_friction_per_atom << ")"
<< " != "
<< cont.size()*Dim << "(" << &cont << ")");
if (old_restitution_per_atom.size() != cont.size()*Dim)
LM_FATAL("inconsistent previous storage of force vector (restitution) "
<< old_restitution_per_atom.size() << "(" << &old_restitution_per_atom << ")"
<< " != "
<< cont.size()*Dim << "(" << &cont << ")");
}
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
void StimulationLangevin<Cont>::gatherWork(Cont & cont) {
Communicator & comm = Communicator::getCommunicator();
CommGroup group = cont.getCommGroup();
if (compute_work_flag){
comm.allReduce(&frictional_work,1,group,"frictional work", OP_SUM);
comm.allReduce(&restitution_work,1,group,"restitution work", OP_SUM);
frictional_work *= 0.5*timestep*code_unit_system.fd2e;
restitution_work *= 0.5*timestep*code_unit_system.fd2e;
cumulated_frictional_work += frictional_work;
cumulated_restitution_work += restitution_work;
friction_restitution_work.empty();
friction_restitution_work.add(temperature);
friction_restitution_work.add(frictional_work);
friction_restitution_work.add(restitution_work);
friction_restitution_work.add(cumulated_frictional_work );
friction_restitution_work.add(cumulated_restitution_work);
friction_restitution_work.copyContainerInfo(cont);
}
}
/* -------------------------------------------------------------------------- */
template <typename Cont>
template <bool computeWork, bool computeWorkPerAtom>
void StimulationLangevin<Cont>::addForce(Cont & cont){
if (computeWork) resetOldWorkPerAtom(cont);
this->frictionCorrection<computeWork,computeWorkPerAtom>(cont);
this->restitutionCorrection<computeWork,computeWorkPerAtom>(cont);
if (computeWork) gatherWork(cont);
}
/* -------------------------------------------------------------------------- */
/* LMDESC LANGEVIN
This stimulator applies a Langevin thermostat to domain. The equation of motion
becomes:\\
\begin{equation}
m\ddot{x} = - \nabla \phi \underbrace{- \frac{m}{damp}{v}}_{f_{damp}} +
\underbrace{\sqrt{\frac{k_{b}\,T\,m}{dt\,damp}}R}_{f_{rand}}
\end{equation}
with $\phi$ the classical potential for interatomic forces, $T$ the temperature
of the heat bath, $damp$ the relaxation time, $v$ the velocity field, $k_b$ the boltzmann constant,
$dt$ the used timestep, $m$ the particle mass and $R$ a vector of random numbers
uniformly distributed between -1/2 and 1/2. \\
Explicitely, this stimulator append $f_{damp}$ and $f_{rand}$ to the force field during stage PRE\_STEP3.
*/
/* LMEXAMPLE STIMULATION thermostat LANGEVIN INPUT md TEMP 100 SEED 32 */
/* LMHERITANCE action_interface */
template <typename Cont>
void StimulationLangevin<Cont>::declareParams(){
Stimulation<Cont>::declareParams();
this->changeDefault("STAGE",PRE_STEP3);
/* LMKEYWORD DAMP
Set the damping parameter expressed in the units of time
*/
this->parseKeyword("DAMP",damp);
/* LMKEYWORD TEMP
Set the temperature desired to be set
*/
this->parseKeyword("TEMP",temperature);
/* LMKEYWORD SEED
Set the seed for the random generator
*/
this->parseKeyword("SEED",seed);
/* LMKEYWORD TIMESTEP
Set the time step needed for random noise contruction
*/
this->parseKeyword("TIMESTEP",timestep);
/* LMKEYWORD WORK_PERATOM
TODO
*/
this->parseTag("WORK_PERATOM",compute_work_peratom_flag,false);
/* LMKEYWORD WORK
Activates the computation of the work done by ther thermostat.
A compute named friction\_restitution:ID is created and allows the visualization
of the work done by the thermostat (ID should be the name given to this compute).
More details available
The work is computed like:\ \
$$dW^n = \int_{x^n}^{x^{n+1}} f \cdot dx$$
which for the case of the velocity verlet scheme is equivalent to
$$dW^n = \frac{1}{2}(f^n + f^{n+1}) v^{n+1/2}$$
Thus when this flag is activated a compute with
with the name "friction\_restitution:STIMULATOR\_NAME"
and 5 entries is registered. These entries are:
\begin{enumerate}
\item Requested temperature
\item frictional\_work: $dW_{damp}^n = \frac{1}{2}(f_{damp}^n + f_{damp}^{n+1}) v^{n+1/2}$
\item restitution\_work: $dW_{restitution}^n = \frac{1}{2}(f_{rand}^n + f_{rand}^{n+1}) v^{n+1/2}$
\item cumulated\_frictional\_work: $W_{damp}^n = \sum_0^n dW_{damp}^n$
\item cumulated\_restitution\_work: $W_{restitution}^n = \sum_0^n dW_{restitution}^n$
\end{enumerate}
*/
this->parseTag("WORK",compute_work_flag,false);
}
/* -------------------------------------------------------------------------- */
DECLARE_STIMULATION(StimulationLangevin,LIST_ATOM_MODEL)
DECLARE_STIMULATION(StimulationLangevin,LIST_CONTINUUM_MODEL)
DECLARE_STIMULATION(StimulationLangevin,LIST_DD_MODEL)
DECLARE_STIMULATION_REFPOINT(StimulationLangevin)
DECLARE_STIMULATION_GENERIC_MESH(StimulationLangevin)
__END_LIBMULTISCALE__

Event Timeline

c4science · Help