Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F82222403
stimulation_langevin.cc
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Sep 10, 06:16
Size
9 KB
Mime Type
text/x-c++
Expires
Thu, Sep 12, 06:16 (2 d)
Engine
blob
Format
Raw Data
Handle
20665981
Attached To
rLIBMULTISCALE LibMultiScale
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
Log In to Comment