1. technology and computing
2. software
3. databases

How to add a new type of interaction to
ESyS-Particle
Steffen Abe
September 30, 2009
Contents
1 Unbonded Interactions
2
1.1 The Interaction and its Parameters . . . . . . . . . . . . . . . . .
2
1.1.1 Interaction group parameter class . . . . . . . . . . . . . .
2
1.1.2 Interaction class . . . . . . . . . . . . . . . . . . . . . . .
3
1.2 The Master/Worker Infrastructure . . . . . . . . . . . . . . . . .
7
1.2.1 The Master . . . . . . . . . . . . . . . . . . . . . . . . . .
7
1.2.2 The Worker . . . . . . . . . . . . . . . . . . . . . . . . . .
8
1.3 The Python Interface . . . . . . . . . . . . . . . . . . . . . . . . .
9
1.3.1 Python Wrapper for the Interaction Group Parameter Class 10
1.3.2 Exporting the Function to create the Interaction Group . 11
2 Bonded Interactions
13
3 Single Particle Interactions
14
4 Wall / Mesh Interactions
15
Introduction
To add a new type of interaction to ESyS-Particle is relatively easy. However,
due to the parallel Master/Worker architecture of the code and the Python interface the necessary bits and pieces are spread relative widely throughout the
source code. This document explains the steps needed to implement a new type
of interaction. There are 4 basic kinds of interactions in ESyS-Particle which
differ in the way some aspects are implemented - unbonded interactions between
pairs of particles, bonded interactions between pairs of particles, interactions between particles and walls or mesh surfaces and body force fields such as Gravity
which are implemented as “single particle interactions”. Unbonded interactions
are those which are generated dynamically for particles which come into contact with each other during the simulation without relying on pre-existing links
between specific pairs of particles.
Currently this tutorial only covers unbonded interactions, but it will be successively expanded to cover the other types of interactions
1
Chapter 1
Unbonded Interactions
This chapter explains how to implement a new type of unbonded interaction
between particles using elastic interactions with a Hertzian contact law (1.1) as
an example.
1.1
The Interaction and its Parameters
To implement the new interactions a header and a implementation file for the
necessary classes are needed, in our example called HertzianElasticInteraction.h and HertzianElasticInteraction.cpp. These files should be located
in the “Model” subdirectory of the source tree. In order to include these files
into the build process of ESyS-Particle both the .h and the .cpp file need to be
added to Makefile.am in the same subdirectory.
For the implementation of a new type of interaction two classes are needed,
one to represent the parameters describing the interactions - the Interaction
group parameters and the class handling the interaction itself, i.e. the Interaction class
1.1.1
Interaction group parameter class
All interaction group parameter classes are derived from class AIGParam. Therefore the file in which a new interaction group parameter class is defined needs
to include ”Model/IGParam.h”. The class needs the following member functions
and variables
• A constructor without parameters
• A constructor taking the name of the interaction group and all interaction
group parameters as arguments. The name is an STL-string.
• A “getTypeString()” function which returns the type of the interactions
as a STL-string.
2
• Member variables for the necessary interaction parameters
In general the interaction parameters would be implemented as public member variables which is not the most “clean” way of doing it in the light of the
principles of object orientation but it saves writing a lot of access functions.
In the example the interaction group parameter class needs two member
variables, one for Young’s modulus and one for the Poisson ratio the definition
(in HertzianElasticInteraction.h) looks like follows:
c l a s s C H e r t z i a n E l a s t i c I G P : public AIGParam
{
public :
double m E ;
double m nu ;
CHertzianElasticIGP ( ) ;
C H e r t z i a n E l a s t i c I G P ( const s t d : : s t r i n g &,double , double ) ;
v i r t u a l s t d : : s t r i n g g e t T y p e S t r i n g ( ) const
{ return ” H e r t z i a n E l a s t i c ” ; }
};
The implementation of the constructors (in HertzianElasticInteraction.cpp
would be as follows:
CHertzianElasticIGP : : CHertzianElasticIGP ( )
: AIGParam ( ) , m E ( 0 . 0 ) , m nu ( 0 . 0 )
{}
C H e r t z i a n E l a s t i c I G P : : C H e r t z i a n E l a s t i c I G P ( const s t d : : s t r i n g&
name , double E , double nu )
: AIGParam ( name ) , m E(E) , m nu ( nu )
{}
while the getTypeString function is already fully defined in the header file.
1.1.2
Interaction class
The interaction class itself is responsible for
• calculating the forces between interacting particles,
• applying these forces to the particles, and
• supplying data relevant to the interaction (such as potential energy, interaction count ....) to the FieldSaver functions.
Currently all interaction classes deal either with external influences on a single
particle (Damping, Gravity) or with interactions between exactly 2 particles.
The classes for interactions between 2 particles, such as the Hertzian elastic
interaction used as an example here, are derived from APairInteraction.
An interaction class needs to contain, as a minimum, the following functions:
3
• A constructor taking 2 particle pointers and a reference to the appropriate
interaction group parameter class - line 23 in the following listing.
• A calcForces() function which calculates the interaction force and applies
it to the particles - line 29 in the listing
• The necessary functions for saving of fields, i.e.
– The functions which map the field access functions to the field names,
get[Scalar,Vector]FieldFunction(...) (lines 9 & 10, related types are
defined in lines 5 & 6) and getCheckedScalarFieldFunction(...) (line 11,
related type defined in line 7). These functions take the field name
(as STL-string) as parameter and return a pointer to the function
to access the respective field. The functions need to be implemented
even if there is no field in the appropriate category (in that case they
just return a NULL pointer and an error message), otherwise the
template instantiation of the class storing the interactions will fail.
– The actual field access functions, in the example : getForce() and
getPotentialEnergy(), see lines 27 and 30 in the listing below.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
c l a s s C H e r t z i a n E l a s t i c I n t e r a c t i o n : public A P a i r I n t e r a c t i o n
{
public :
typedef double ( C H e r t z i a n E l a s t i c I n t e r a c t i o n : : ∗
S c a l a r F i e l d F u n c t i o n ) ( ) const ;
typedef Vec3 ( C H e r t z i a n E l a s t i c I n t e r a c t i o n : : ∗
V e c t o r F i e l d F u n c t i o n ) ( ) const ;
typedef p a i r <bool , double> ( C H e r t z i a n E l a s t i c I n t e r a c t i o n : : ∗
C h e c k e d S c a l a r F i e l d F u n c t i o n ) ( ) const ;
s t a t i c S c a l a r F i e l d F u n c t i o n g e t S c a l a r F i e l d F u n c t i o n ( const
s t r i n g &) ;
s t a t i c V e c t o r F i e l d F u n c t i o n g e t V e c t o r F i e l d F u n c t i o n ( const
s t r i n g &) ;
static CheckedScalarFieldFunction
g e t C h e c k e d S c a l a r F i e l d F u n c t i o n ( const s t r i n g &) ;
private :
double
double
Vec3
double
Vec3
m E ; //!< Young ’ s modulus
m nu ; //!< P o i s s o n r a t i o
m f o r c e ; // c a c h i n g f o r c e f o r E p o t
m dn ; // c a c h i n g d i s p l a c e m e n t f o r E p o t
m cpos ; // c e n t e r p o s i t i o n
public :
typedef C H e r t z i a n E l a s t i c I G P ParameterType ;
4
C H e r t z i a n E l a s t i c I n t e r a c t i o n ( C P a r t i c l e ∗ , C P a r t i c l e ∗ , const
C H e r t z i a n E l a s t i c I G P &) ;
virtual ˜ C H e r t z i a n E l a s t i c I n t e r a c t i o n ( ) {};
23
24
25
26
27
28
29
30
31
v i r t u a l Vec3 g e t P o s ( ) const { return m cpos ; } ;
double g e t P o t e n t i a l E n e r g y ( ) const ;
v i r t u a l void c a l c F o r c e s ( ) ;
Vec3 g e t F o r c e ( ) const ;
};
Calculating and Applying the Force
For the Hertzian Contact law used as an example the normal force fn between
two interacting particles can be calculated as
!
p
2E Rij
δn3/2
(1.1)
fn =
3(1 − ν 2 )
Where E is Young’s modulus (assumed to be identical for both particles), ν is
the Poisson ratio, Rij the effective radius at the contact point, i.e. 1/Rij =
1/Ri + 1/Rj and δn is the displacement of the particles relative to their equilibrium distance. The implementation of the force calculation could look like the
following code fragment:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
void C H e r t z i a n E l a s t i c I n t e r a c t i o n : : c a l c F o r c e s ( )
{
Vec3 D=m p1−>g e t P o s ( )−m p2−>g e t P o s ( ) ;
double d i s t=D∗D;
double e q d i s t=m p1−>getRad ( )+m p2−>getRad ( ) ;
i f ( d i s t <( e q d i s t ∗ e q d i s t ) ) {
double R i j = 1 . 0 / ( 1 . 0 / m p1−>getRad ( ) +1.0/ m p2−>getRad ( ) ) ;
d i s t=s q r t ( d i s t ) ;
m dn=d i s t −e q d i s t ;
Vec3 d i r=D. u n i t ( ) ;
m f o r c e=d i r ∗ (m E∗ s q r t ( R i j ) ) / ( 2 . 0 ∗ ( 1 . 0 − m nu∗m nu ) ) ∗pow (
m dn , 1 . 5 ) ;
Vec3 pos=m p2−>g e t P o s ( ) +(m p2−>getRad ( ) / e q d i s t ) ∗D;
m p2−>a p p l y F o r c e ( m f o r c e , pos ) ;
m p1−>a p p l y F o r c e ( −1.0∗ m f o r c e , pos ) ;
}
}
where lines 3–6 determine if the current distance between the particles is smaller
than the sum of their radii, i.e. if two particles are actually in contact. If they
are in contact, that is if the test in line 6 succeeds, the interaction force is
calculated (line 11 basically implements equation 1.1) and is the applied to the
particles in lines 13 and 14. One needs to be careful to make sure the sign of
the force is correct.
5
The elastic strain energy of an interaction can be calculated as
!
p
4E Rij
estrain =
δn5/2
15(1 − ν 2 )
(1.2)
or, if the force has been calculated from equation 1.1 already
estrain =
2
fn δn
5
(1.3)
Implementing equation 1.3 results in the code fragment
1
2
3
4
double C H e r t z i a n E l a s t i c I n t e r a c t i o n : : g e t P o t e n t i a l E n e r g y ( ) const
{
return 0 . 4 ∗ m f o r c e . norm ( ) ∗m dn ;
}
Field Access Functions
The interaction, as implemented in the example, has two “fields” which can be
accessed by the FieldSaver functions from a simulations script, the force and
the elastic strain energy (also called “potential energy”). Therefore two access functions have been implemented, getForce() where the return value is a
vector, i.e. a Vec3 and getPotentialEnergy() where the return value is a scalar
(double). Therefore getScalarFieldFunction needs to return a pointer to getPotentialEnergy() if the name parameter is “potential energy” and a NULL-pointer
otherwise. Likewise the getVectorFieldFunction returns a pointer to getForce()
if the name parameter is “force”. Implementation see below.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
CHertzianElasticInteraction : : ScalarFieldFunction
C H e r t z i a n E l a s t i c I n t e r a c t i o n : : g e t S c a l a r F i e l d F u n c t i o n ( const
s t r i n g& name )
{
CHertzianElasticInteraction : : ScalarFieldFunction sf ;
i f ( name==” p o t e n t i a l e n e r g y ” ) {
s f=&C H e r t z i a n E l a s t i c I n t e r a c t i o n : : g e t P o t e n t i a l E n e r g y ;
} else {
s f=NULL;
s t d : : c e r r << ”ERROR − i n v a l i d name f o r i n t e r a c t i o n s c a l a r
a c c e s s f u n c t i o n ” << s t d : : e n d l ;
}
return s f ;
}
CHertzianElasticInteraction : : VectorFieldFunction
C H e r t z i a n E l a s t i c I n t e r a c t i o n : : g e t V e c t o r F i e l d F u n c t i o n ( const
s t r i n g& name )
6
16
17
18
19
20
21
22
23
24
25
26
27
{
CHertzianElasticInteraction : : VectorFieldFunction vf ;
i f ( name==” f o r c e ” ) {
v f=&C H e r t z i a n E l a s t i c I n t e r a c t i o n : : g e t F o r c e ;
} else {
v f=NULL;
c e r r << ”ERROR − i n v a l i d name f o r i n t e r a c t i o n v e c t o r
a c c e s s f u n c t i o n ” << e n d l ;
}
return v f ;
}
1.2
The Master/Worker Infrastructure
In order to instantiate an interaction group of a given type, a number of functions
are needed in the Master and Worker. In the Master, functions are needed to
pack the interaction group parameters and send them to the Workers whereas
in the Worker, functions are needed to receive and unpack the parameters and
to instantiate the interaction group.
1.2.1
The Master
In the file LatticeMaster.h (in the “Parallel” subdirectory) it is necessary to add
an #include for the appropriate header file of the new interaction, i.e. in the example #include “Model/HertzianElasticInteraction.h”. Additionally the following
functions need to be added, declaration in LatticeMaster.h” and implementation
in LatticeMaster.cpp
• The function to add the interaction group for all particles, void addPairIG(
const CHertzianElasticIGP&).
• The function to add the interaction group only for particles with defined
tags void addTaggedPairIG(const CHertzianElasticIGP&,int,int,int,int).
To be more precise, these functions pack the parameters, send them to the
Workers and signal the Workers to add the interaction group.
The two functions (see listings below) both consist of the declaration of a
class for a message which will then be sent to the workers (line 3 in addPairIG,
line 5 in addTaggedPairIG), adding the necessary parameters to the message and
finally broadcasting the message to all workers (lines 7 and 14 respectively).
The differences are that in addPairIG the message class is of type IGPCommand
whereas in addTaggedPairIG it is of type TaggedIGPCommand. Also, in the
tagged case it is necessary to add the tag and mask parameters to the message
(lines 9–12).
7
1
2
3
4
5
6
7
8
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
void C L a t t i c e M a s t e r : : addPairIG ( const C H e r t z i a n E l a s t i c I G P &prms
)
{
IGPCommand igpCmd ( getGlobalRankAndComm ( ) ) ;
igpCmd . appendIGP ( prms ) ;
igpCmd . append ( prms . m E) ;
igpCmd . append ( prms . m nu ) ;
igpCmd . b r o a d c a s t ( ) ;
}
void C L a t t i c e M a s t e r : : addTaggedPairIG ( const C H e r t z i a n E l a s t i c I G P
&prms ,
i n t tag1 , i n t mask1 ,
i n t t a g 2 , i n t mask2 )
{
TaggedIGPCommand igpCmd ( getGlobalRankAndComm ( ) ) ;
igpCmd . appendIGP ( prms ) ;
igpCmd . append ( prms . m E) ;
igpCmd . append ( prms . m nu ) ;
igpCmd . append ( t a g 1 ) ;
igpCmd . append ( mask1 ) ;
igpCmd . append ( t a g 2 ) ;
igpCmd . append ( mask2 ) ;
igpCmd . b r o a d c a s t ( ) ;
}
1.2.2
The Worker
In the Worker, functions are needed to receive and unpack the parameters and
to instantiate the interaction group. For this purpose the following additions
need to be made to the file SubLattice.hpp:
• Add include of the header file of the new interaction type, i.e. “Model/HertzianElasticInteraction.h”
• Inside the existing function doAddPIG() add the necessary code to instantiate the interaction group, i.e.
– Add a test for the name of the new interaction type: if (type=”HertzianElastic”)
( see line 2 in the listing below)
– Unpack the interaction group parameters from the message buffer
(lines 4 and 5)
– If a tagged interaction group is to be instantiated, i.e. if “tagged”
is true (line 6) then it is also necessary to unpack the tag and mask
parameters from the message buffer (lines 7–10)
8
– Instantiate a “pair interaction storage” of the interaction group (line
11 for the tagged case, line 13 otherwise). To choose the appropriate type of “pair interaction storage” class (PIS), it is necessary to
distinguish if the interactions contain persistent data which need to
be exchanged between different worker processes should the particles
involved in the interaction move to a different domain. The options
are ParallelInteractionStorage NE : “no exchange” (interactions without persistent data) or ParallelInteractionStorage ED : “exchange, dynamic” (interactions with persistent data which are created dynamically). An added “ T” on the class name stands for the “tagged”
versions of these classes.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
....
} e l s e i f ( type==” H e r t z i a n E l a s t i c ” ) {
CHertzianElasticIGP heigp ;
h e i g p . m E=p a r a m b u f f e r . p o p d o u b l e ( ) ;
h e i g p . m nu=p a r a m b u f f e r . p o p d o u b l e ( ) ;
i f ( tagged ) {
i n t t a g 1=p a r a m b u f f e r . p o p i n t ( ) ;
i n t mask1=p a r a m b u f f e r . p o p i n t ( ) ;
i n t t a g 2=p a r a m b u f f e r . p o p i n t ( ) ;
i n t mask2=p a r a m b u f f e r . p o p i n t ( ) ;
n e w p i s=new P a r a l l e l I n t e r a c t i o n S t o r a g e N E T <T,
C H e r t z i a n E l a s t i c I n t e r a c t i o n >(m ppa , heigp , tag1 , mask1 ,
tag2 , mask2 ) ;
} else {
n e w p i s=new P a r a l l e l I n t e r a c t i o n S t o r a g e N E <T,
C H e r t z i a n E l a s t i c I n t e r a c t i o n >(m ppa , h e i g p ) ;
}
r e s=true ;
}
....
N.B. The strategy for the instantiation of interaction groups described above
only applies to “dynamic” interaction groups, i.e. those where interactions between individual particles are created dynamically from neigbour relations between particles. Bonded interaction groups, where the interactions between individual particles are determined by pre-defined connectivity relation between the
particles, are different. Those interaction groups are created via addBondedIG()
1.3
The Python Interface
In order to be able to use the new type of interaction from an ESyS-Particle
simulation script the Python interface to ESyS-Particle needs two extensions, a
Python wrapper for the interaction group parameter class and a Python wrapper
for the call to instantiate the interaction group.
9
1.3.1
Python Wrapper for the Interaction Group Parameter Class
In order to be able to create an instance of the interaction group parameters it is
necessary to expose the class implemented in section 1.1 to the python interface.
For this purpose the files InteractionParamsPy.h and InteractionParamsPy.cpp need to be extended by
• Adding include to the header file of the new interaction type (i.e. “Model/HertzianElasticInteraction.h”) to InteractionParamsPy.h
• Adding the definition of a wrapper class for the interaction group parameter class (in the example CHertzianElasticIGP in InteractionParamsPy.h.
The wrapper class needs to have as its only member function a constructor
which takes the same arguments as the wrapped class
1
2
3
4
5
c l a s s H e r t z i a n E l a s t i c P r m s P y : public C H e r t z i a n E l a s t i c I G P
{
public :
H e r t z i a n E l a s t i c P r m s P y ( const s t d : : s t r i n g &,double , double ) ;
};
• Add the implementation of the constructor in InteractionParamsPy.cpp
which in most cases will just pass the arguments on to the constructor of
the wrapped class
1
2
3
4
5
H e r t z i a n E l a s t i c P r m s P y : : H e r t z i a n E l a s t i c P r m s P y ( const s t d : :
s t r i n g &name ,
double E ,
double nu )
: C H e r t z i a n E l a s t i c I G P ( name , E , nu )
{}
• Again in InteractionParamsPy.cpp, add the wrapper class to the interface
which is exported to Python via boost. The two functions exposed are
the constructor defined above and the getName function which is inherited
from the wrapped class.
1
2
3
4
5
b o o s t : : python : : c l a s s <H e r t z i a n E l a s t i c P r m s P y , b o o s t : : python : :
b a s e s <I n t e r a c t i o n P r m s P y > >(
” HertzianElasticPrms ” ,
” Parameters f o r H e r t z i a n e l a s t i c c o n t a c t i n t e r a c t i o n s .
”,
b o o s t : : python : : i n i t <const s t d : : s t r i n g &, double ,
double>(
(
10
6
7
8
9
10
11
a r g ( ”name” ) ,
a r g ( ”E” ) ,
a r g ( ”nu” )
),
” @type name : s t r i n g \n”
”@kwarg name : Name a s s i g n e d t o t h i s group o f
i n t e r a c t i o n s . \ n”
” @type E : f l o a t \n”
”@kwarg E : Young ’ s modulus used f o r f o r c e
c a l c u l a t i o n . \ n”
” @type nu : f l o a t \n”
”@kwarg nu : p o i s s o n r a t i o used f o r f o r c e c a l c u l a t i o n
. \ n”
12
13
14
15
16
17
18
19
20
21
)
)
. def (
”getName” ,
&H e r t z i a n E l a s t i c P r m s P y : : getName ,
b o o s t : : python : : r e t u r n v a l u e p o l i c y <b o o s t : : python : :
c o p y c o n s t r e f e r e n c e >()
)
;
22
23
1.3.2
Exporting the Function to create the Interaction
Group
In order to be able to create the interaction group it is necessary to add a
function to the python interface which calls LatticeMaster::addPairIG with the
appropriate parameters. To implement this function the LsmMpiPy class needs
to be extended by:
• Adding a forward declaration of class HertzianElasticPrmsPy in LsmMpiPy.h
• Adding a declaration void createHertzianElasticIG(const HertzianElasticPrmsPy &prms) in LsmMpiPy.h
• Implementing void createHertzianElasticIG in LsmMpiPy.cpp
1
2
3
4
5
6
7
void LsmMpiPy : : c r e a t e H e r t z i a n E l a s t i c I G (
const H e r t z i a n E l a s t i c P r m s P y &prms
)
{
getNameTypeMap ( ) [ prms . getName ( ) ] = prms . g e t T y p e S t r i n g ( ) ;
g e t L a t t i c e M a s t e r ( ) . addPairIG ( prms ) ;
}
• Adding the appropriate variant of createInteractionGroup to the boost interface (void exportLsm() in LsmMpiPy.cpp)
11
1
2
3
4
5
6
7
8
9
10
. def (
” createInteractionGroup ” ,
&LsmMpiPy : : c r e a t e H e r t z i a n E l a s t i c I G ,
( a r g ( ”prms” ) ) ,
” C r e a t e s a group o f H e r t z i a n e l a s t i c i n t e r a c t i o n s . \ n”
” @type prms : L{ e s y s . lsm . I n t e r a c t i o n P r m s ”
”<e s y s . lsm . LsmPy . I n t e r a c t i o n P r m s >}\n”
”@kwarg prms : An o b j e c t d e s c r i b i n g t h e type o f
i n t e r a c t i o n and”
” any p a r a m e t e r s a s s o c i a t e d with t h a t i n t e r a c t i o n . \ n”
)
12
Chapter 2
Bonded Interactions
TBD
13
Chapter 3
Single Particle Interactions
TBD
14
Chapter 4
Wall / Mesh Interactions
TBD
15