LibMultiScale & AMEL User Guide
Guillaume ANCIAUX |
The LibMultiScale Library has been designed to be a high performance tool to study
the mutiscale methods currently employed on material simulations. The component
design of the project distinct the base code components (Continuum mechanics code
implemented through finite element, meshless, XFEM or molecular dynamics codes) from
the coupling components. Other services like dumpers ( logging the state of the
simulated system on files or memory, enabling filtering based on values averaging or
geometrical criterion...), or stimulators that implement initial conditions (temperature,
initial displacements) are provided.
At this day, the only stable coupling method is based on the Bridging Method.
The molecular dynamics component can be implemented by Lammps or Stamp while the
Continuum mechanics component use a constitutive law based on the Cauchy-Born Rule
and is implemented by a code developed in our team.
The project is actually maintained by Guillaume ANCIAUX.
Chapitre 2 Getting LibMultiScale and Compiling
Please consult the INSTALL.txt page for details on how to fully install
LibMultiScale and its client AMEL (Atomistic Multi-domain coupling with External
finite elements for LibMultiScale) that is downloadable at the page:
http://gforge.inria.fr/docman/?group_id=915
The LibMultiScale Project is actually maintained by Guillaume ANCIAUX in
a collaboration between the ScalApplix Project, INRIA Futurs, and the
CEA Ile-de-France (contract CEA/DIF N4600084494/P6H32). The license that governs the
LibMultiScale rights is a CECILL-C which can be consulted at the address :
http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
Chapitre 3 Preprocessing of the models
In order to use the LibMultiScale and its client AMEL,
one should prepare the models wanted to be coupled.
For instance, in the actual case, a mesh should be
prepared and a mono-atomic crystal.
In order to illustrate the discussion, we will describe
the construction of a concrete case.
The case studied will be a simple rectangular piece
of material with a planar and circular(composed geometries) bridged
zone. The figure 3.1
present the global aspect of such a coupled domain.
It is a square domain
with a side of L = 1000 ⋅ r0 where r0 is
the averaged inter-atomic distance. We firstly want
to build an mono-atomic crystal of argon, meaning that r0 = 1.22452143391761666.10−10.
The bridged zone
placed in the center of the domain would have an height
size h = L/10 while length is l = L/3.
Bridged zone size should not be smaller than
10 r0 in order to avoid wave reflexions. Thus it represents a bridging size of L/100.
In order to prevent reflexions, we use here a bridged size of 2 L /100.
First we prepare the set of atoms constituting the crystal
we want to modelize. Then we build the mesh using the appropriate tools.
The numerical models are associated with configuration files that will
be described in the following.
The different stages of this tutorial can be found in the examples directory of
LibMultiScale release.
In order to build the model we create a file casCustom that will initially
contain the useful values we defined. The key words are not defined systematically
in the following. For a full description, one should report to the glossary at the
end of the document. Firstofall, we need a global section where models, and other
important data will be assigned. This is done by opening a section block:
#general section
Section MultiScale RealUnits
...
endSection
The syntax is Section NAME UNITS, where the name is identified
by the MultiScale and the units declare the units used in the
section. Inside of that block, we set the declaration of reusable variables.
We naturally define the dimensions of the system by:
#definition of the global numbers
#inter-atomic distance
LET r0 = 1.2245214339176166619211e-10
#total domain size
LET L = 1000 * r0
#height of molecular domain
LET h = L / 10
#square width of molecular domain
LET l = L / 3
#the size of the bridged zone
LET bridgedSize = 2 * L / 100
#the dimensions of the coarsed zone
LET l2 = l + bridgedSize + h
LET h2 = bridgedSize + h
Any simulation will have to define its system dimension by:
DIMENSION 2
and the internal units system used by the codes. This precision is important
as the base code components for molecular mechanics have restricted unity systems.
UNITCODE RealUnits
The next section will explain how to configure the atomic section
before the continuum section and its mesh.
3.1 preparing the crystal
In order to prepare the atomic zone, the first step in the construction
is to build a squared/cubic zone full of atoms adequately disposed.
Then we will carve into that cubic zone to obtain the desired geometry.
In our concrete case we want
to modelize 2D Argon crystalline. We need to define a Bravais lattice.
There are two atoms in the reference Lattice, with coordinates
(in reduced unities):
0.0 0.0 0.5
0.5 0.5 0.5
The lattice is rectangular and follows the vectors:
1.22452143391761666e-10 0.0
0.0 2.1209333385024e-10 0.0
We will describe two ways to obtain that cubic zone. The first one will
use Stamp to replicate Bravais Lattice while the second will use a given
script to generate input for Lammps.
3.2 Generating with Stamp
In the following we may define by = 2.1209333385024.10−10.
The basic information concerning a crystal is synthetized in a file
called FAtomes. A more precise
description of this Stamp file can be found in the Stamp Documentation.
For the desired crystal, this file is:
*******************************************
***** les unites sont le MKS **************
** sauf les masses sont molaires(Kg/mol) **
*******************************************
nbAtomes 1
nom AR
masse 39.95e-3
sigma 0.11e-9
epsilon 1.6567944e-21
interPotentiel PotentielLJ
****************************
**** coordonnees du motif 1.23157195294235e-10 **
****************************
Motif
1.22452143391761666e-10 0.0 0.0
0.0 2.1209333385024e-10 0.0
0.0 0.0 7e-10
***********************************
**** coordonnees reduires [0..1] **
* X **** Y ***** Z *****At ********
***********************************
Structure
2
0.0 0.0 0.5 AR
0.5 0.5 0.5 AR
Then a Stamp configuration file must be produced. In the same way,
the appropriate structure of this file can be found in the Stamp
documentation. We give here an example:
FichierAtomes config-stamp/FAtomes2D
LongueurInput angstrom
TempsInput ps
TemperatureInput K
EnergieInput eV
EnergieOutput eV
Ensemble NVE
Schema verlet
Structure Fichier
Decoupage 474 81 1
Potentiel LJ
Rcut 2.93885144140227
Trequis 0.0
Deltatemps 0.001
RandSeed 17
ConditionX Libre
ConditionY Libre
ConditionZ Libre
*EcrireStatProc 1
StepLimit 40000
StepAvg 200
Protection 100
The Decoupage keyword control how many times the Bravais lattice is
duplicated in each directions. Accounting the bridging zone size,
a rectangular zone with dimensions l' × h' needs to be filled with atoms:
|
|
| l' × h' |
= |
(l + bridgedSize + h) × (bridgedSize + h) |
(3.1) |
| |
= |
(L ⋅ 136/300) × (L ⋅ 12 /100)
|
(3.2) |
|
The rounded-off (ceiled) rectangular zone should be:
l' × h' = 474 r0 × 140 r0
We then determine the lattice factors:
fx = 474 fy = 140 ⋅ 0.5773502691896277986 = 81
3.3 Generating with Script application
That second method will produce an input file for Lammps. In the crystalgen directory
there exist scripts helpful to produce crystal input files for Lammps. The interesting one for the
described example is gencrystal-lammps.pl. The various scripts for pre/post treatment can be found
on the project page1. We first need to provide a bravais lattice
basis vectors file:
argon2D.bravais:
1.22452143391761666 0 0
0 2.1209333385024 0
0 0 7
Then we need to provide the list of atoms present in an unity lattice in reduced coordinates:
maille-argon2D.xyz:
2
ARGN
AR 0.0 0.0 0.5
AR 0.5 0.5 0.5
Then calling the following command generate the input file for Lammps:
./gencrystal-lammps.pl argon2D.bravais maille-argon2D.xyz 474 81 1 > custom.lammps
Then a lammps configuration file must be provided. Let it be in.lj.custom
with the following containing:
# 2d Lennard-Jones Argon
units real
dimension 2
atom_style atomic
boundary s s p
read_data config-lammps/data.custom
#velocity all create 80.0 87287 loop geom
pair_style lj/cut 2.93885144140227
pair_coeff 1 1 0.238492861193117 1.1 2.93885144140227
neighbor 0.4 bin
neigh_modify every 20 check no
fix 1 all nve
restart 1000 restart.lammps
thermo 10
run 50
Remark : here the run keyword is present. For use with LibMultiScale it should not be present
as it is LibMultiScale who will performs the integration steps.
3.4 Carving the coarse model geometry
We have prepared a rectangular 2D crystal. We want to control the dimensions
using the LibMultiScale facilities. To do that, we must complete our configuration
file adding a dedicated section. The coarse geometric description will be a cube
described by its corner:
GEOMETRY CUBE 2 -l2/2 l2/2 -h2/2 h2/2 0 0
The syntax is GEOMETRY ID TYPE DIMENSION PARAMETERS. Parameters are
described in the glossary. The geometries can be defined anywhere in the configuration
files, as well as variables declarations.
#Stamp section
Section Stamp RealUnits
DOMAIN_NAME md
STAMP_FILE config-stamp/CAScfcCustom
GEOMETRY 4 CUBE 2 -l2/2 l2/2 \\
-h2/2 h2/2 0 0
DOMAIN_GEOMETRY 4
DOMAIN_BOUNDARY_GEOMETRY 4
endSection
#Lammps section
Section Lammps RealUnits
DOMAIN_NAME md
LAMMPS_FILE config-lammps/in.lj.custom
GEOMETRY 4 CUBE 2 -l2/2 l2/2 \\
-h2/2 h2/2 0 0
DOMAIN_GEOMETRY 4
DOMAIN_BOUNDARY_GEOMETRY 4
endSection
In this sections we defined the name of the domain by section DOMAIN_NAME name, the configuration
file used with STAMP/LAMMPS_FILE filename, and the geometry of the domain by DOMAIN_GEOMETRY ID.
The boundary must be provided also, but practically it is used in the coupling process, and then need
not to be coherent at this stage. Then the total geometry is assigned also to the boundary.
We now have to declare the model in the global section by adding MD_CODE STAMP 0 or
MD_CODE LAMMPS 0 (where 0 is the declared ID of the model) giving the global configuration file:
#general section
Section MultiScale RealUnits
#definition of the global numbers
#inter-atomic distance
LET r0 = 1.2245214339176166619211e-10
#total domain size
LET L = 1000 * r0
#height of molecular domain
LET h = L / 10
#square width of molecular domain
LET l = L / 3
#the size of the bridged zone
LET bridgedSize = 2 * L / 100
#the dimensions of the coarsed zone
LET l2 = l + bridgedSize + h
LET h2 = bridgedSize + h
DIMENSION 2
UNITCODE RealUnits
MD_CODE STAMP 0
DUMPER PARAVIEW 1
endSection
#Stamp section
Section Stamp RealUnits
DOMAIN_NAME md
STAMP_FILE config-stamp/CAScfcCustom
GEOMETRY 4 CUBE 2 -l2/2 l2/2 -h2/2 h2/2 0 0
DOMAIN_GEOMETRY 4
DOMAIN_BOUNDARY_GEOMETRY 4
endSection
Figure 3.2:
Coarsed crystal domain generated at stage 1.
Distance unities are angstroems.
We also added a DUMPER keyword to request the creation of a dumper component
with implementation through PARAVIEW output (the first parameter is the frequency
at which the dumper should generate outputs).
Then, launching the command ./AMEL casCustom 1 would produce
in ../analysis directory the initial dumps in VTK Xml format2. The figure 3.2
illustrate what should be obtained at this stage.
Remark: The main difference between Lammps and Stamp component is
the internal unit system used. Lammps have an internal system that use
angstroem metric system unit. Stamp uses SI metric system unit.
The difference is notified by the command UNITCODE option (option = AtomUnits/RealUnits)
In order to refine the geometric description, we have to use composed geometries.
In that manner, we can decompose the total atomic geometry in three geometries as
presented by the figure 3.3. Two circles with radius h2/2 and centers
(−l/2,0) and (l/2,0) respectively, and a rectangle with (−l/2,−h/2),
(−l/2,h/2),(l/2,−h/2) and (l/2,h/2) as corners.
This can be done by the following definitions:
#Stamp section
Section Stamp RealUnits
DOMAIN_NAME md
STAMP_FILE config-stamp/CAScfcCustom
#definition of sub geometries
#values used to define the half circles and the central square
LET ANGLE1 = pi/2
LET ANGLE2 = 3 * ANGLE1
GEOMETRY 2 BALL 2 -l/2 0 0 h2/2 ANGLE1 ANGLE2
GEOMETRY 3 BALL 2 l/2 0 0 h2/2 ANGLE2 ANGLE1
GEOMETRY 4 CUBE 2 -l/2 l/2 -h2/2 h2/2 0 0
#union of the geometries
GEOMETRY 5 UNION 2 3
GEOMETRY 6 UNION 5 4
DOMAIN_GEOMETRY 6
DOMAIN_BOUNDARY_GEOMETRY 6
endSection
Figure 3.3:
refined molecular domain.
We have prepared the crystal for the multi-scale model. We now have to focus
on the preparation of the continuum domain.
3.5 Preparing the mesh
The current way used to prepare the final mesh needs the software
Gmsh3that can build 2D/3D meshes. We can consider
a scale factor that will be passed to LibMultiScale
as a parameter. This is an option for convenient
build of the mesh. We then consider that L=2.
The first stage of building consists in placing the relevant
boundary points of the domain. These points will be added
for the external boundary:
(-1,-1) (-1,1) (1,-1) (1,1)
and for the internal boundaries:
(-0.33,-0.1) (-0.33,0.1) (0.33,-0.1) (0.33,0.1)
Then we have to introduce the points to define the
bridged zone by the points:
(-0.33,-0.14) (-0.33,0.14) (0.33,-0.14) (0.33,0.14)
Finally, in order to describe the curvatures we define
the points:
(-0.43,0) (-0.47,0) (0.43,0) (0.47,0)
The Gmsh file you should obtain is:
Point(1) = {-1,-1,0,0.1};
Point(2) = {-1,1,0,0.1};
Point(3) = {1,-1,0,0.1};
Point(4) = {1,1,0,0.1};
Point(5) = {-0.33,-0.1,0,0.1};
Point(6) = {-0.33,0.1,0,0.1};
Point(7) = {0.33,-0.1,0,0.1};
Point(8) = {0.33,0.1,0,0.1};
Point(9) = {-0.33,-0.14,0,0.1};
Point(10) = {-0.33,0.14,0,0.1};
Point(11) = {0.33,-0.14,0,0.1};
Point(12) = {0.33,0.14,0,0.1};
Point(13) = {-0.43,0,0,0.1};
Point(14) = {-0.47,0,0,0.1};
Point(15) = {0.43,0,0,0.1};
Point(16) = {0.47,0,0,0.1};
Point(17) = {-0.33,0,0,0.1};
Point(18) = {0.33,0,0,0.1};
Figure 3.4:
skeleton points for the mesh
that should be obtained
Now, we have to build the segments, and then the surfaces
in order to let Gmsh generate the elements.
The following lines constitute the geometric information of the
continuum domain (eventually, they can be built using the graphical
interface of Gmsh):
Circle(1) = {8,18,15};
Circle(2) = {15,18,7};
Circle(5) = {10,17,14};
Circle(6) = {6,17,13};
Circle(7) = {13,17,5};
Circle(8) = {9,17,14};
Circle(9) = {11,18,16};
Circle(10) = {12,18,16};
Line(11) = {6,8};
Line(12) = {5,7};
Line(13) = {9,11};
Line(14) = {10,12};
Line Loop(15) = {10,-9,-13,8,-5,14};
Line Loop(16) = {11,1,2,-12,-7,-6};
Plane Surface(17) = {15,16};
Line(18) = {4,3};
Line(19) = {1,2};
Line(20) = {2,4};
Line(21) = {3,1};
Line Loop(22) = {18,21,19,20};
Plane Surface(23) = {22,15};
Figure 3.5:
segment and surfaces for the mesh
that should be obtained
Now, one could play with the weights in order to inform Gmsh
of the surfacic constraints on elements we want to maintain.
Typically we want the mesh to be well refined near the center
of the domain, i.e. the future atomic zone. In that manner we
affect weights of 0.01 to points 5 to 12.
Now using the Gmsh mesh construction facility we may obtain
the result presented on the figure ?? by typing:
gmesh -o mymesh.unv -format unv -2 file.geo
where file.geo is the file containing all the precedent
geometric informations. Then a mesh file formatted with
UNV standard loadable by LibMesh is produced. This file
should be placed in LibMultiScale/run/config-libmesh/
in the standard file tree place.
Remark : one can refine in a different way by playing with
the weights placed onto the initial points. Creating other internal
surfaces allow refinement to be better controlled. For example,
setting the same weight to all points induce the mesh presented
by figure ??.
Figure 3.6:
Final mesh obtained
by Gmsh
We open a new section block with name LibMesh:
#LibMesh section
Section LibMesh RealUnits
...
endSection
Then we need to set the name of the domain:
DOMAIN_NAME mymesh
Then we have to set the size of the domain.
The global size has been decided to be L = 1000 r0.
Then we notify the file that the libMesh will have to load:
MESH_FILE config-libmesh/mymesh.unv L/2
The following parameter is the scaling factor. It is set in order
to scale the mesh file to the real domain size. As we choose L=2
in the mesh file, we may want to scale by half the real domain
size.
We also have to describe the geometry of the domain. This is done using
the substraction operator and the shrink operator. The shrink operator
takes any basic geometry (no combined geometries)
and reduces it following given parameters.
The geometry declaration take the form:
#full domain cube
GEOMETRY 17 CUBE 2 -L/2 L/2 -L/2 L/2 0 0
#shrinking the two circles
GEOMETRY 18 SHRINK 5 bridgedSize 0 0
#shrinking the central cube only on the y axis
GEOMETRY 19 SHRINK 4 0 bridgedSize 0
GEOMETRY 20 UNION 18 19
GEOMETRY 21 SUB 17 20
DOMAIN_GEOMETRY 21
DOMAIN_BOUNDARY_GEOMETRY 21
We have to define the constitutive relation that need to be used
as well as the mass of the monoatomic crystal.
We want to use a constitutive law defined by the Cauchy-Born rule
with a lennard-jones potential. We precise
the parameters by:
CONSTITUTIVELAW CAUCHY-BORN 2.1*r0
LENNARDJONES r0 1.6567944e-21 0.11e-9 [lattice_file]
MASS 39.95e-3
For the Cauchy-Born constitutive law to work well, one should specify a bravais
lattice section into a given file (or in the same file as the LibMesh section
if no file name was specified). This section should look like:
#crystal section
Section Lattice RealUnits
LET sqrt3 = 1.732050807568877294
BRAVAIS 0 r0 0.0 0.0
BRAVAIS 1 r0/2 r0/2*sqrt3 0.0
ATOMTYPE Ar 39.95e-3 0
NBATOM 1
ATOM Ar 0.0 0.0 0.0
endSection
which is the declaration of a mono-atomic crystal of Argon with
two atoms per lattice.
Also, the atomic mass is required in order to compute the cauchy-born rule.
And finally we declare the timestep that will be used:
TIMESTEP 1e-15
In order to dump the finite element domain, the same way that
we did for the molecular dynamics model, we have to declare it
in the global section adding:
ELAST_CODE LIBMESH 1
At the stage 3, the configuration file should look like:
File Edit Options Buffers Tools Help
#global section
Section MultiScale RealUnits
#definition of the global numbers
#inter-atomic distance
LET r0 = 1.2245214339176166619211e-10
#total domain size
LET L = 1000 * r0
#height of molecular domain
LET h = L / 10.
#square width of molecular domain
LET l = L / 3
#the size of the bridged zone
LET bridgedSize = 2 * L / 100
#the dimensions of the coarsed zone
LET l2 = l + bridgedSize + h
LET h2 = 2*bridgedSize + h
DIMENSION 2
UNITCODE RealUnits
MD_CODE STAMP 0
ELAST_CODE LIBMESH 1
DUMPER PARAVIEW 1
endSection
#Stamp section
Section Stamp RealUnits
DOMAIN_NAME md
STAMP_FILE config-stamp/CAScfcCustom
#definition of sub geometries
#values used to define the half circles and the central square
LET ANGLE1 = pi/2
LET ANGLE2 = 3 * ANGLE1
GEOMETRY 2 BALL 2 -l/2 0 0 h2/2 ANGLE1 ANGLE2
GEOMETRY 3 BALL 2 l/2 0 0 h2/2 ANGLE2 ANGLE1
GEOMETRY 4 CUBE 2 -l/2 l/2 -h2/2 h2/2 0 0
#union of the geometries
GEOMETRY 5 UNION 2 3
GEOMETRY 6 UNION 5 4
DOMAIN_GEOMETRY 6
DOMAIN_BOUNDARY_GEOMETRY 6
endSection
#LibMesh section
Section LibMesh RealUnits
CONSTITUTIVELAW CAUCHY-BORN 2.1*r0
MESH_FILE config-libmesh/mymesh.unv L/2
DOMAIN_NAME mymesh
GEOMETRY 17 CUBE 2 -L/2 L/2 -L/2 L/2 0 0
GEOMETRY 18 SHRINK 5 bridgedSize 0 0
GEOMETRY 19 SHRINK 4 0 bridgedSize 0
GEOMETRY 20 UNION 18 19
GEOMETRY 21 SUB 17 20
DOMAIN_GEOMETRY 21
DOMAIN_BOUNDARY_GEOMETRY 21
LENNARDJONES r0 1.6567944e-21 0.11e-9
TIMESTEP 1e-15
MASS 39.95e-3
endSection
#crystal section
Section Lattice RealUnits
LET sqrt3 = 1.732050807568877294
BRAVAIS 0 r0 0.0 0.0
BRAVAIS 1 r0/2 r0/2*sqrt3 0.0
ATOMTYPE Ar 39.95e-3 0
NBATOM 1
ATOM Ar 0.0 0.0 0.0
endSection
Reusing the command ./AMEL casCustom 1 will now
produce one more vtu file describing the continuum model giving the
result presented by figure ??.
3.6 Preparing the bridging zone
At this step of the construction we have defined both models.
The coupling is activated by the simple keyword
BRIDGING_CODE TYPE ID1 ID2 OPTIONS...
The bridging that will be described is based on the Bridging Method
(the only one stable at this time). We need to manage a Bridged Zone
where the DOFs of both domains will be weighted.
The detection of the DOFs involved in the bridging zone is made by a simple test: any DOFs
for which its initial position is inside the intersection of both model geometries,
is considered as a coupled DOF. The ponderation process has a detached procedure.
In order to help
the LibMultiScale to know the shape of the weighting gradients,
the user have to provide a geometry ID. This geometry
can be composed only by operator UNION. All primitive
are accepted.
We then have to describe the Bridging Zone by four sub geometries:
GEOMETRY 100 CUBE 2 -l/2 l/2 -h2/2 -h/2 0 0
GEOMETRY 101 CUBE 2 -l/2 l/2 h/2 h2/2 0 0
GEOMETRY 102 BALL 2 -l/2 0 h/2 h2/2 ANGLE1 ANGLE2
GEOMETRY 103 BALL 2 l/2 0 h/2 h2/2 ANGLE2 ANGLE1
GEOMETRY 104 UNION 100 101
GEOMETRY 105 UNION 102 103
GEOMETRY 106 UNION 104 105
BRIDGING_CODE BRIDGING 0 1 106
If we submit this configuration file, AMEL
generates additional VTU files:
-
recatomesatoms_weigth.vtu
- recelemsmesh_weigth.vtu
- grid.vtu
- fictifsatomesatoms_weigth.vtu
- fictifselemsmesh_weigth.vtu
- grid-fictifs.vtu
The two first ones are presented by figures ??
and ??.
It presents the DOFs that were detected inside of the bridging zone
based on the geometric description given. And the color field represents
the associated weight.
Figure 3.8:
Mesh part detected inside if the bridging zone at stage4
Figure 3.9:
Atoms detected inside if the bridging zone at stage4
It can be seen that it did not worked correctly. Due to the curvature
of the circles, atoms may have been inside of both domain geometric description, but not inside
of the desired elements. Then the LibMultiScale have selected additional elements
creating the “elemental picks” around the mesh. Additionally some nodes have a weight of zero
indicating that there were not inside the Bridging Zone geometry (default weight is zero).
Then, as a reaction to the problem, one should reduce with a delta parameter the atomic domain
and grow equivalently the bridging zone geometry.
The grid VTU dump files are the grids used to speed up the computation of
atom-element association.
Concerning the other files, they treat a part of the coupling we did not approach until now.
On a geometry like the one presented, boundary conditions cannot be maintained on the free
surfaces4.
The surface effect involved can be controlled by the coupling component. This is
the goal of the boundary geometries in the molecular dynamics cases. One should define a zone
where the atoms are concerned by the force field rupture near the free surface. Then the coupling component
will detect these atoms and owning elements to impose to these atoms a displacement interpolated
by the continuum domain. The two last “fictifs” files present the detection. As we declared
the boundary geometry to be the same as the domains, then all the atoms in the bridged zone
will have an imposed displacement. We must refine the boundary geometry to conform the wanted
atoms (i.e. free surface within a given cutoff radius range).
The total configuration file at stage 5 would be:
#global section
Section MultiScale RealUnits
#definition of the global numbers
#inter-atomic distance
LET r0 = 1.2245214339176166619211e-10
#total domain size
LET L = 1000 * r0
#height of molecular domain
LET h = L / 10.
#square width of molecular domain
LET l = L / 3
#the size of the bridged zone
LET bridgedSize = 2 * L / 100
#the dimensions of the coarsed zone
LET l2 = l + bridgedSize + h
LET h2 = 2*bridgedSize + h
#delta parameter
LET delta = r0/6
DIMENSION 2
UNITCODE RealUnits
MD_CODE STAMP 0
ELAST_CODE LIBMESH 1
GEOMETRY 100 CUBE 2 -l/2 l/2 -h2/2-delta -h/2 0 0
GEOMETRY 101 CUBE 2 -l/2 l/2 h/2 h2/2+delta 0 0
GEOMETRY 102 BALL 2 -l/2 0 h/2 h2/2+delta ANGLE1 ANGLE2
GEOMETRY 103 BALL 2 l/2 0 h/2 h2/2+delta ANGLE2 ANGLE1
GEOMETRY 104 UNION 100 101
GEOMETRY 105 UNION 102 103
GEOMETRY 106 UNION 104 105
BRIDGING_CODE BELYTSCHKO 0 1 106
DUMPER PARAVIEW 1
endSection
#Stamp section
Section Stamp RealUnits
DOMAIN_NAME md
STAMP_FILE config-stamp/CAScfcCustom
#definition of sub geometries
#values used to define the half circles and the central square
LET ANGLE1 = pi/2
LET ANGLE2 = 3 * ANGLE1
GEOMETRY 2 BALL 2 -l/2 0 0 h2/2-delta ANGLE1 ANGLE2
GEOMETRY 3 BALL 2 l/2 0 0 h2/2-delta ANGLE2 ANGLE1
GEOMETRY 4 CUBE 2 -l/2 l/2 -h2/2-delta h2/2-delta 0 0
#union of the geometries
GEOMETRY 5 UNION 2 3
GEOMETRY 6 UNION 5 4
#definition of the boundary
LET rcut = r0 * 2.6
GEOMETRY 7 SHRINK 5 rcut 0 0
GEOMETRY 8 SHRINK 4 0 rcut 0
GEOMETRY 9 UNION 7 8
GEOMETRY 10 GROW 5 r0 0 0
GEOMETRY 11 GROW 4 0 r0 0
GEOMETRY 12 UNION 10 11
GEOMETRY 13 SUB 12 9
DOMAIN_GEOMETRY 6
DOMAIN_BOUNDARY_GEOMETRY 13
endSection
#LibMesh section
Section LibMesh RealUnits
CONSTITUTIVELAW CAUCHY-BORN 2.1*r0
MESH_FILE config-libmesh/mymesh.unv L/2
DOMAIN_NAME mymesh
GEOMETRY 17 CUBE 2 -L/2 L/2 -L/2 L/2 0 0
GEOMETRY 18 SHRINK 5 bridgedSize 0 0
GEOMETRY 19 SHRINK 4 0 bridgedSize 0
GEOMETRY 20 UNION 18 19
GEOMETRY 21 SUB 17 20
DOMAIN_GEOMETRY 21
DOMAIN_BOUNDARY_GEOMETRY 21
LENNARDJONES r0 1.6567944e-21 0.11e-9
TIMESTEP 1e-15
MASS 39.95e-3
endSection
The user must control that the ponderation is correct, that the surface effect will affect
the appropriate atoms. All the configuration of the model was done in a sequential aspect.
Simulation can yet be launched changeing the number of step to the desired one:
./AMEL casCustom 10000
The description of the stimulation and dumpers needed to alter and monitor the simulation
is given in the glossary. Concerning the processor mapping, the user can consult the next section.
3.7 Configuring the distribution over a cluster of processors
At this day, there is only one type of distribution possible. Two models cannot coexist
on the same processor. More precisely, each processor will be the host for one model and
only one model as a UNIX processus. Then the option for the user is to decide how many processors are
attributed on a given model. This is done by the keyword PROCESSORS. A sample configuration
using 4 processors, 3 for the MD model and 1 for the FE model should be:
COM DISTRIBUTED
PROCESSORS 1 3
PROCESSORS 0 1
where the COM keyword expose the mapping solution. As we said, this is actually
the only one available.
- 1
- http://gforge.inria.fr/projects/libmultiscale/
- 2
-
readable by the opensource software ParaView
- 3
- Any mesh generator software that produce readable formats for LibMesh
should work fine with LibMultiScale.
- 4
- in simpler geometries, periodic boundary conditions could be applied, or truncated
potential
In this chapter are described all the rules and keywords for the parser of AMEL.
4.1 Section syntax
A section should be declared as follow:
Section NAME UNITS
...
endSection
where the name of the section is specified by NAME, and UNITS specify
the units that will be used in that section. This is mandatory in order
to have an automatic units conversion, ensuring common config files for
different internal units (hence for different codes). UNITS can have the
value RealUnits (S.I.) and AtomUnits (angstroems,...).
4.2 Variable declaration
LET varname = mathexpression
This is used to declare any reusable variable in the configuration files.
The variables are defined in the reading order. If a keyword, like MD_CODE
or ELAST_CODE is placed, the variables defined in the called section will
be defined for the next line of the currect section.
PRINT var
Used to dump the value of a variable. Debug usage for preparation of the
domains, geometry ... etc.
4.3 Geometry keywords
GEOMETRY ID TYPE OPTIONS ...
The description of the different types follows now.
Declaration of a BALL geometry.
GEOMETRY ID BALL 1 CX RMIN RMAX
GEOMETRY ID BALL 2 CX CY RMIN RMAX TMIN TMAX
GEOMETRY ID BALL 3 CX CY CZ RMIN RMAX TMIN TMAX PMIN PMAX
The declaration is dimension dependant. Considering the spherical coordinates,
RMIN,RMAX from the radius offset. TMIN TMAX, the angle offset for teta coordinate.
PMIN, PMAX, is the angle offset on the last coordinate.
Declaration of a CUBE geometry
GEOMETRY ID CUBE 1 XMIN XMAX HX
GEOMETRY ID CUBE 2 XMIN XMAX YMIN YMAX HX HY
GEOMETRY ID CUBE 3 XMIN XMAX YMIN YMAX ZMIN ZMAX HX HY HZ
The declaration is demention dependant. The parameters are the coordinates of
the bounding box. The last parameters (HX,HY,HZ), specify a possible hole
in the bounding box allowing various geometries like presented by figure ??.
Figure 4.1:
Here are presented the geometries declared for two different use of
the hole parameter.
Declaration of an ELLIPSOID geometry
GEOMETRY ID ELLISPOID 1 CX A1
GEOMETRY ID ELLISPOID 2 CX CY A1 A2
GEOMETRY ID ELLISPOID 3 CX CY CZ A1 A2 A3
CX,CY,CZ specify the center of the ellipsoid. the parameters A1, A2 and A3
specify the axis length of the ellipsoid.
Declaration of a composed geometry
GEOMETRY ID INTER ID1 ID2
GEOMETRY ID UNION ID1 ID2
GEOMETRY ID SUB ID1 ID2
INTER produce the intersection of geometries ID1 and ID2 while
UNION produce their union. SUB substract the geometry ID2 from
the ID1 geometry.
Shrinking/Growing a geometry
GEOMETRY ID SHRINK ID1 FX FY FY
GEOMETRY ID GROW ID1 FX FY FZ
The FX,FY,FZ parameters should positive. Only composed
geometries, CUBE and BALL may be used by this operator.
In the BALL case, only FX is relevant, modifying the RMAX value.
In cube case, the bounding box is reduced in each direction following
FX,FY,FZ.
4.4 Atomic domain keywords
Giving the domain name:
DOMAIN_NAME name
Defining the domain geometries:
DOMAIN_GEOMETRY ID
DOMAIN_BOUNDARY_GEOMETRY ID
This is the total domain geometric description and
the boundary geometry. The boundary geometry is only used
by the coupling components to correct surface effects.
Setting behavior for parallel optimisation (optional).
MIGRATION value
value can be enable or disable. It has sense only for parallel
simulations. By default it is enabled. At present time this keyword is deprecated.
4.4.1 Lammps keywords
LAMMPS_FILE config-file
4.4.2 Stamp keywords
STAMP_FILE config-file
4.4.3 MD1D keywords
This model works only in 1D.
RCUT value
Defines the cutoff range for the potential.
R0 value
Defines the interatomic distance of the simulated monoatomic crystal.
NUM_ATOMES value
Defines the number of atoms wanted. The generated crystal will be centered.
MASS value
Defines the atomic mass.
TIMESTEP value
Defines the timestep of simulation.
EPSILON value
SIGMA value
Defines the parameters for the Lennard-Jones potential.
4.4.4 DummyMD keywords
This model is used to replay a saved simulation. Thus a model have
to be reloaded and then frames are read instead of calculated.
DUMMY_TYPE type
Specify the type of simulation to reload. It can be LAMMPS, LAMMPS or MD1D.
DUMMY_FILE file
Specify the config file used to initiate the domain. It is hence a STAMP or
LAMMPS or MD1D configuration file.
4.5 Continuum keywords
Giving the domain name:
DOMAIN_NAME name
Defining the domain geometries:
DOMAIN_GEOMETRY ID
DOMAIN_BOUNDARY_GEOMETRY ID
This is the total domain geometric description and
the boundary geometry. The boundary geometry is only used
by the coupling components to correct surface effects.
ELEM_SIZE value value2D value3D
This specify the uniform size of the elements in each direction.
This is used for automatic grid meshing of the continuum domain.
CONSTITUTIVELAW type params
This specify to the constitutive law controller that is wanted. At present time,
possible values for type are CAUCHY-BORN and LAME.
For CAUCHY-BORN type, one parameter has to be specified and is a cutoff radius.
A command line specifying potential to be used is then required.
For LAME type, lame parameters are required ( Poisson modulus should be precised first ).
TIMESTEP value
Specify the explicit timestep value.
LENNARDJONES r0 epsilon sigma lattice_file
Specify the parameters for the lennard jones parameter. First one
is interatomic distance and others are classical lennard-jones parameters.
If lattice_file is specified then the Lattice section is searched in that
file. Otherwise the lattice section is searched in the current file.
ANALYTIC_EAM a p re ksi q lattice_file
Specify the parameters for the analytic eam potential. The computed potential
has the formula:
W = Σi ≠j a exp(-p(r_ij/re -1)) + ksi √Σi≠j exp(-2q(r_ij/re -1))
EAM potential_file lattice_file
Specify the parameters for the eam potential. The potential file
has the standard format described at the address http://www4.eas.asu.edu/cms/potentials/main/main.htm
MASS value
Mass used to compute initial density.
MESH_FILE file_name scalefactor
This notify the filename of the mesh to be loaded. Then the mesh will
be scaled by the given factor.
SHIFT x y z
This will translate the nodes by the vector (x,y,z), giving opportunity to
tune the mesh position with atoms.
RESTART_FILE filename
Specify a file to reload positions, velocity and forces stored at a given state.
It should be a zipped xml file with valid tags. Such files can be generated by
dumpers (see Dumper section).
RELOAD_PARTITION filename
not yet documented
BOUNDARY_CONDITION geom_id type regexp1D egexp2D regexp3D
The boundary conditions specified by this option affect the geometry given
by the geom_id parameter.
The boundary conditions can actually be of two types : DISP or PBC.
-
DISP imposes a given displacement given by the regular expression that should be
like ux= 2 uy= 0 uz= 0 .
- PBC declares a periodic boundary condition. It is functional only for auto-generated
meshes (grid meshes). Then a fusion of different nodes is operated in order
to have a periodic behavior of the simulation. the regexp parameters specify
the dimensions in which this operation is to be performed (example :
X=1 Y=0 Z=0 for a periodic boundary condition along X axis).
PRESSION value
This is used to apply a pression onto the X-Z surfaces of a continuum domain.
The value of the pression is to be precised in the used code units (no automatic conversion).
Obviously this option was created for a really specific job. Not to be used commonly.
4.6 Lattice keywords
This following keywords are to be used into a Lattice section.
BRAVAIS coordinate axisx axisy axisz
Specify a bravais vector identified by coordinate. There should be
as much bravais vectors as the dimension of the space.
ATOMTYPE type mass charge
Declare one atom type that is to be associated with mass and charge.
NBATOM value
Specify the number of atoms into the lattice. There should be
as use of ATOM keyword (see the following) as NBATOM specify.
ATOM type x y z
Declare one atom into the lattice. The coordinates are lattice relative.
Type parameter should refer to a previous declaration by ATOMTYPE.
4.7 Global keywords
This keywords should be used in the MultiScale section (which should be unique for
a given case).
UNITCODE value
This command is used to specify the internal units that will be used. The possible values are
AtomsUnits and RealUnits.
MD_CODE type ID conffile
This keyword defines a molecular dynamics domain. The type defines the code that
will be used to compute the trajectories of the atoms (MD1D, LAMMPS, STAMP or DUMMY).
The ID parameter identify by a number the domain hence declared. The conffile
parameter specify where to find the information about the configuration of this domain, which
will have to contain a section about that type of domain. If no file is specified, then the current
file is parsed for such a section.
ELAST_CODE type ID conffile
Idem for continuum mechanics code. Values for type can only be LIBMESH at present time.
BRIDGING_CODE TYPE IDMD IDCONTINUUM geomID [GRID_DIVISIONX value]
[GRID_DIVISIONY value]
[GRID_DIVISIONZ value]
This keyword declares a bridging zone coupled by a method of type TYPE. Types can be:
-
BELYTSCHKOFULL : belytschko method but with no condensation of the constraint matrix. Functional only in sequential.
- BELYTSCHKO : belytschko method but with condensation of the constraint matrix. It is the more generic approach.
IDMD and IDCONTINUUM are domain identifiers that must have been introduced by the command MD_CODE and ELAST_CODE
respectively. The geomID identify a geometry that will be used for the description of the weighting functions.
In the computation of the association between atoms and finite elements, a spatial grid is used to speed up
the research. The optional GRID_DIVISIONX parameters define the granularity of such a spatial grid along the X-axis.
Idem for the GRID_DIVISIONY and GRID_DIVISIONZ parameters.
4.8 Dumper keywords
Declaration of a dump component:
DUMPER TYPE FREQUENCY OPTIONS ...
all dumper have their frequency to be specified. The options
may change one to an other. Still, there are general options:
-
FILTERGEOM geom_id : used to request a filter on geometrical criterion. Only
dofs into the specified geometry will be dumped.
- FILTERDISP value : The criterion for filter is base on the displacement field.
If the displacement value, in absolute, exceed the given value, then the dof is
dumped.
- FILTERVEL value : same thing for velocity field.
- FILTERFORCE value : same thing for force field.
- PREFIX dir : specify the directory where the dumped file should be created.
Default is ../analysis/
The Paraview Dumper:
DUMPER PARAVIEW FREQUENCY [DISP] [VEL] [FORCE] [CONSTRAINT]
[BASE64] [TEXT] [ID] [MASS]
[PROC] [COMPRESSED] [P0]
the options are explained:
DISP : request dumping the displacements
VEL : request dumpin the velocities
PROC : request dumping the distribution
FORCE : request dumping the force/acceleration field
CONSTRAINT : request dumping the constraint field
DEFORMATION : request dumping the deformation field (only for continuum)
BASE64 : request using in the base64 mode(default)
TEXT : request dumping in the text mode
COMPRESSED : compression of dumped file is activated (only if zlib support has been compiled)
ID : internal identity number of each dof is dumped as a normal field used for debugging only
MASS : mass associated with each dof is dumped
P0 : initial position is dumped
DUMPER EPSN FREQUENCY xml_file
Unmaintained.
DUMPER DUMPER1D FREQUENCY
Produce 1D dump files analysable via gnuplot scripts. the script analysePVF.pl
should be used to produce videos or fig files.
DUMPER FIG FREQUENCY
If libMesh has been compiled with the fig support this should produce
(only in sequential) a valid xfig file (helpfull to produce documentation)
DUMPER SCOTCH FREQUENCY
This will produce (in sequential only) scotch graph format files.
Mainly used in internal to produce graph taking into account the coupling
stage.
DUMPER DUMPERECIN FREQUENCY
This will produce a single txt file with the kinetic energy for each timestep.
DUMPER XYZ FREQUENCY
This will produce an XYZ file for the atoms (only).
DUMPER LAMMPS FREQUENCY
This will produce a lammps file for the atoms (only).
This file will be loadable for a simulation using lammps module.
DUMPER RESTART FREQUENCY
This will produce restart files for the used codes. One restart file
produced by lammps may not be used with Stamp. The format is xml with
dedicated tags and optionally (depend on the compilation of LibMultiScale)
gzipped.
4.9 Stimulation keywords
This presents all the stimulation that is possible to inject.
STIMULATION type options
This declares the use of a given stimulation. There are two kinds,
initial and permanent stimulation. The following details the principal
types.
STIMULATION GAUSS lWave intensity geometryID
This produces an initial displacement with a gaussian expression onto
a given geometry. The value at the exponential term is given by lWave
while the intensity value modify the factor against the exponential
(intensity ⋅ e−(2π r/lWave)2 where r is the distance to the center of the
given geometry).
STIMULATION SIN lWave intensity geometryID
This produces an initial velocity field following a sinusoid function
intensity ⋅ sin(2π r /lWave).
STIMULATION MODULATE lWave intensity geometryID
This produces a modulation over the actual displacement by a given frequency.
The applied formula is u += intensity ⋅ cos(2π r/lWave).
STIMULATION RAZ
This is used when restarting a situation. One may want to reset displacement,
velocities, forced to zero. the initial positions are then set
STIMULATION MATIERELANCEE axe intensity geomID
This produce a uniaxial velocity field over the given geometry. The velocity
field follows a gradient oriented along a given axes (1=x, 2=y, 3=z) growing from
the center of the geometry to the external surfaces.
The intensity is the higher value imposed.
STIMULATION STABILISE
This implements the stabilization algorithm of Melh (minimization heuristic that kill velocities
which do not follow the potential slope).
4.10 Parallel keywords
In order to declare the distribution over a cluster, one should
set the following keywords.
COM DISTRIBUTED
PROCESSOR ID_DOMAIN nb_procs
In the future, other methods of distribution, accounting to the geometrical
situation will be available.
Ce document a été traduit de LATEX par HEVEA