1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
|
/// \file
/// Leapfrog time integrator
#include "timestep.h"
#include "CoMDTypes.h"
#include "linkCells.h"
#include "parallel.h"
#include "performanceTimers.h"
static void advanceVelocity(SimFlat* s, int nBoxes, real_t dt);
static void advancePosition(SimFlat* s, int nBoxes, real_t dt);
/// Advance the simulation time to t+dt using a leap frog method
/// (equivalent to velocity verlet).
///
/// Forces must be computed before calling the integrator the first time.
///
/// - Advance velocities half time step using forces
/// - Advance positions full time step using velocities
/// - Update link cells and exchange remote particles
/// - Compute forces
/// - Update velocities half time step using forces
///
/// This leaves positions, velocities, and forces at t+dt, with the
/// forces ready to perform the half step velocity update at the top of
/// the next call.
///
/// After nSteps the kinetic energy is computed for diagnostic output.
double timestep(SimFlat* s, int nSteps, real_t dt)
{
for (int ii=0; ii<nSteps; ++ii)
{
startTimer(velocityTimer);
advanceVelocity(s, s->boxes->nLocalBoxes, 0.5*dt);
stopTimer(velocityTimer);
startTimer(positionTimer);
advancePosition(s, s->boxes->nLocalBoxes, dt);
stopTimer(positionTimer);
startTimer(redistributeTimer);
redistributeAtoms(s);
stopTimer(redistributeTimer);
startTimer(computeForceTimer);
computeForce(s);
stopTimer(computeForceTimer);
startTimer(velocityTimer);
advanceVelocity(s, s->boxes->nLocalBoxes, 0.5*dt);
stopTimer(velocityTimer);
}
kineticEnergy(s);
return s->ePotential;
}
void computeForce(SimFlat* s)
{
s->pot->force(s);
}
void advanceVelocity(SimFlat* s, int nBoxes, real_t dt)
{
for (int iBox=0; iBox<nBoxes; iBox++)
{
for (int iOff=MAXATOMS*iBox,ii=0; ii<s->boxes->nAtoms[iBox]; ii++,iOff++)
{
s->atoms->p[iOff][0] += dt*s->atoms->f[iOff][0];
s->atoms->p[iOff][1] += dt*s->atoms->f[iOff][1];
s->atoms->p[iOff][2] += dt*s->atoms->f[iOff][2];
}
}
}
void advancePosition(SimFlat* s, int nBoxes, real_t dt)
{
for (int iBox=0; iBox<nBoxes; iBox++)
{
for (int iOff=MAXATOMS*iBox,ii=0; ii<s->boxes->nAtoms[iBox]; ii++,iOff++)
{
int iSpecies = s->atoms->iSpecies[iOff];
real_t invMass = 1.0/s->species[iSpecies].mass;
s->atoms->r[iOff][0] += dt*s->atoms->p[iOff][0]*invMass;
s->atoms->r[iOff][1] += dt*s->atoms->p[iOff][1]*invMass;
s->atoms->r[iOff][2] += dt*s->atoms->p[iOff][2]*invMass;
}
}
}
/// Calculates total kinetic and potential energy across all tasks. The
/// local potential energy is a by-product of the force routine.
void kineticEnergy(SimFlat* s)
{
real_t eLocal[2];
eLocal[0] = s->ePotential;
eLocal[1] = 0;
for (int iBox=0; iBox<s->boxes->nLocalBoxes; iBox++)
{
for (int iOff=MAXATOMS*iBox,ii=0; ii<s->boxes->nAtoms[iBox]; ii++,iOff++)
{
int iSpecies = s->atoms->iSpecies[iOff];
real_t invMass = 0.5/s->species[iSpecies].mass;
eLocal[1] += ( s->atoms->p[iOff][0] * s->atoms->p[iOff][0] +
s->atoms->p[iOff][1] * s->atoms->p[iOff][1] +
s->atoms->p[iOff][2] * s->atoms->p[iOff][2] )*invMass;
}
}
real_t eSum[2];
startTimer(commReduceTimer);
addRealParallel(eLocal, eSum, 2);
stopTimer(commReduceTimer);
s->ePotential = eSum[0];
s->eKinetic = eSum[1];
}
/// \details
/// This function provides one-stop shopping for the sequence of events
/// that must occur for a proper exchange of halo atoms after the atom
/// positions have been updated by the integrator.
///
/// - updateLinkCells: Since atoms have moved, some may be in the wrong
/// link cells.
/// - haloExchange (atom version): Sends atom data to remote tasks.
/// - sort: Sort the atoms.
///
/// \see updateLinkCells
/// \see initAtomHaloExchange
/// \see sortAtomsInCell
void redistributeAtoms(SimFlat* sim)
{
updateLinkCells(sim->boxes, sim->atoms);
startTimer(atomHaloTimer);
haloExchange(sim->atomExchange, sim);
stopTimer(atomHaloTimer);
for (int ii=0; ii<sim->boxes->nTotalBoxes; ++ii)
sortAtomsInCell(sim->atoms, sim->boxes, ii);
}
|
