Homec4science
Diffusion LibMultiScale (features/reshape_compute_python) Browse

libmultiscale/src/stimulation/stimulation_langevin.cc
27ca80ef7f29features/reshape_comput…

Last Change
Raw File

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 "stimulation_langevin.hh"
#include "factory_multiscale.hh"
#include "lib_continuum.hh"
#include "lib_dd.hh"
#include "lib_md.hh"
#include "lm_common.hh"
#include "math_tools.hh"
#include "reference_manager.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
StimulationLangevin::StimulationLangevin(const std::string &name)
: LMObject(name), friction_restitution_work("friction_restitution:" + name),
old_friction_per_atom("old_friction_per_atom:" + name),
old_restitution_per_atom("old_restitution_per_atom:" + name),
friction_work_per_atom("friction_work_per_atom:" + name),
restitution_work_per_atom("restitution_work_per_atom:" + name) {
damp = 100;
temperature = 100;
timestep = 1.;
seed = 0;
compute_work_flag = false;
// initialize the compute
friction_restitution_work.assign(5, 0);
friction_restitution_work.name_computed.push_back("temperature");
friction_restitution_work.name_computed.push_back("friction_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_friction_work = 0;
cumulated_restitution_work = 0;
// register for computes
FilterManager::getManager().addObject(friction_restitution_work);
}
/* -------------------------------------------------------------------------- */
StimulationLangevin::~StimulationLangevin() {}
/* -------------------------------------------------------------------------- */
template <typename Cont> void StimulationLangevin::stimulate(Cont &cont) {
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::resetOldWorkPerAtom(Cont &cont) {
constexpr UInt Dim = Cont::Dim;
if (this->compute_work_peratom_flag)
if (friction_work_per_atom.size() == 0) {
friction_work_per_atom.resize(cont.size() * Dim);
restitution_work_per_atom.resize(cont.size() * Dim);
}
if (old_friction_per_atom.size() == 0) {
old_friction_per_atom.resize(cont.size() * Dim);
old_restitution_per_atom.resize(cont.size() * 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::gatherWork(Cont &cont) {
// Communicator &comm = Communicator::getCommunicator();
CommGroup &group = cont.getCommGroup();
if (compute_work_flag) {
group.allReduce(&friction_work, 1, "frictional work", OP_SUM);
group.allReduce(&restitution_work, 1, "restitution work", OP_SUM);
friction_work *= 0.5 * timestep * code_unit_system.fd2e;
restitution_work *= 0.5 * timestep * code_unit_system.fd2e;
cumulated_friction_work += friction_work;
cumulated_restitution_work += restitution_work;
friction_restitution_work.clear();
friction_restitution_work.push_back(temperature);
friction_restitution_work.push_back(friction_work);
friction_restitution_work.push_back(restitution_work);
friction_restitution_work.push_back(cumulated_friction_work);
friction_restitution_work.push_back(cumulated_restitution_work);
friction_restitution_work.copyContainerInfo(cont);
}
}
/* -------------------------------------------------------------------------- */
template <bool computeWork, bool computeWorkPerAtom, typename Cont>
void StimulationLangevin::addForce(Cont &cont) {
if (computeWork)
resetOldWorkPerAtom(cont);
this->frictionCorrection<computeWork, computeWorkPerAtom>(cont);
this->restitutionCorrection<computeWork, computeWorkPerAtom>(cont);
if (computeWork)
gatherWork(cont);
}
/* -------------------------------------------------------------------------- */
template <bool computeWork, bool computeWorkPerAtom, typename Cont>
void __attribute__((noinline))
StimulationLangevin::frictionCorrection(Cont &cont) {
constexpr UInt Dim = Cont::Dim;
Real conv1 = code_unit_system.mv_t2f / damp;
UInt index = 0;
friction_work = 0;
for (auto &&at : cont) {
Real gamma1 = -conv1 * at.mass();
for (UInt i = 0; i < Dim; ++i, ++index) {
Real friction_term = gamma1 * at.velocity()[i];
if (computeWork)
accountFrictionWork<computeWork, computeWorkPerAtom>(
friction_term, at.velocity()[i], index);
at.force()[i] += friction_term;
}
}
}
/* -------------------------------------------------------------------------- */
template <bool computeWork, bool computeWorkPerAtom, typename Cont>
void __attribute__((noinline))
StimulationLangevin::restitutionCorrection(Cont &cont) {
constexpr UInt Dim = Cont::Dim;
Real conv2 = sqrt(24. * temperature * boltzmann * code_unit_system.kT2fd *
code_unit_system.m_tt2f_d / damp / timestep);
restitution_work = 0;
UInt index = 0;
for (auto &&at : cont) {
Real gamma2 = conv2 * sqrt(at.mass());
for (UInt i = 0; i < Dim; ++i, ++index) {
Real restitution_term = gamma2 * (MathTools::uniform() - 0.5);
if (computeWork)
accountRestitutionWork<computeWork, computeWorkPerAtom>(
restitution_term, at.velocity()[i], index);
at.force()[i] += restitution_term;
}
}
}
/* -------------------------------------------------------------------------- */
template <bool computeWork, bool computeWorkPerAtom>
inline void StimulationLangevin::accountFrictionWork(Real friction_term,
Real velocity,
UInt index) {
Real old_friction_term = old_friction_per_atom[index];
friction_work += (friction_term + old_friction_term) * velocity;
old_friction_per_atom[index] = friction_term;
if (computeWorkPerAtom)
friction_work_per_atom[index] += friction_work;
}
/* -------------------------------------------------------------------------- */
template <bool computeWork, bool computeWorkPerAtom>
inline void StimulationLangevin::accountRestitutionWork(Real restitution_term,
Real velocity,
UInt index) {
Real old_restitution_term = old_restitution_per_atom[index];
restitution_work += (restitution_term + old_restitution_term) * velocity;
old_restitution_per_atom[index] = restitution_term;
if (computeWorkPerAtom)
restitution_work_per_atom[index] += restitution_work;
}
/* -------------------------------------------------------------------------- */
/* 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 */
void StimulationLangevin::declareParams() {
StimulationInterface::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);
}
/* -------------------------------------------------------------------------- */
void StimulationLangevin::init() { MathTools::setSeed(seed); }
/* -------------------------------------------------------------------------- */
DECLARE_STIMULATION_MAKE_CALL(StimulationLangevin)
__END_LIBMULTISCALE__
c4science · Help