Line data Source code
1 : // Copyright (C) 2015 Technische Universitaet Muenchen
2 : // This file is part of the Mamico project. For conditions of distribution
3 : // and use, please see the copyright notice in Mamico's main folder, or at
4 : // www5.in.tum.de/mamico
5 :
6 : template <class LinkedCell, unsigned int dim>
7 0 : coupling::UsherParticleInsertion<LinkedCell, dim>::UsherParticleInsertion(unsigned int insertDeleteMassEveryTimestep, double rSigmaCoeff,
8 : double meanPotentialEnergyFactor, double uOverlapCoeff, double stepRefCoeff,
9 : unsigned int iterMax, unsigned int restartMax, double tolerance,
10 : double offsetFromOuterBoundary,
11 : coupling::interface::MDSolverInterface<LinkedCell, dim>* const mdSolverInterface)
12 : : coupling::ParticleInsertion<LinkedCell, dim>(insertDeleteMassEveryTimestep),
13 : #ifdef USHER_DEBUG
14 : _energyInserted(0), _energyRemoved(0), _ZhouEnergyInserted(0), _ZhouEnergyRemoved(0), _particlesInserted(0), _particlesRemoved(0),
15 : #endif
16 0 : _mdSolverInterface(mdSolverInterface),
17 0 : _usherParams(rSigmaCoeff, meanPotentialEnergyFactor, uOverlapCoeff, stepRefCoeff, iterMax, restartMax, tolerance, offsetFromOuterBoundary) {
18 0 : }
19 :
20 : template <class LinkedCell, unsigned int dim>
21 0 : typename coupling::ParticleInsertion<LinkedCell, dim>::Action coupling::UsherParticleInsertion<LinkedCell, dim>::insertDeleteMass(
22 : coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell, const tarch::la::Vector<dim, double>& couplingCellPosition,
23 : const tarch::la::Vector<dim, double>& couplingCellSize, const tarch::la::Vector<dim, double>& meanVelocity, const double& temperature,
24 : const coupling::BoundaryForceController<LinkedCell, dim>& boundaryForceController) {
25 0 : typename coupling::ParticleInsertion<LinkedCell, dim>::Action action = coupling::ParticleInsertion<LinkedCell, dim>::NoAction;
26 0 : const double moleculeMass(_mdSolverInterface->getMoleculeMass());
27 : // if we have enough mass left in the cell, try to insert a particle and
28 : // remove one particle mass
29 0 : if (cell.getMicroscopicMass() >= moleculeMass) {
30 0 : action = insertParticle(cell, couplingCellPosition, couplingCellSize, meanVelocity, temperature, boundaryForceController);
31 0 : if (action == coupling::ParticleInsertion<LinkedCell, dim>::Insertion) {
32 0 : cell.addMicroscopicMass(-moleculeMass);
33 : }
34 : // if we need to delete a particle, we try this and - in case it works - add
35 : // a particle mass to the negative buffer
36 0 : } else if (cell.getMicroscopicMass() <= -moleculeMass) {
37 0 : action = deleteParticle(cell, boundaryForceController);
38 0 : if (action == coupling::ParticleInsertion<LinkedCell, dim>::Deletion) {
39 0 : cell.addMicroscopicMass(moleculeMass);
40 : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
41 : std::cout << "Delete particle: Success" << std::endl;
42 : #endif
43 : }
44 : }
45 :
46 0 : return action;
47 : }
48 :
49 : template <class LinkedCell, unsigned int dim>
50 0 : typename coupling::ParticleInsertion<LinkedCell, dim>::Action coupling::UsherParticleInsertion<LinkedCell, dim>::insertParticle(
51 : coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell, const tarch::la::Vector<dim, double>& couplingCellPosition,
52 : const tarch::la::Vector<dim, double>& couplingCellSize, const tarch::la::Vector<dim, double>& meanVelocity, const double& temperature,
53 : const coupling::BoundaryForceController<LinkedCell, dim>& boundaryForceController) {
54 0 : coupling::datastructures::Molecule<dim> molecule;
55 :
56 : const typename coupling::ParticleInsertion<LinkedCell, dim>::Action action =
57 0 : findParticlePosition(cell, couplingCellPosition, couplingCellSize, molecule, boundaryForceController);
58 :
59 : // if insertion was successful, initialise velocity according to temperature
60 : // in coupling cell
61 0 : if (action == coupling::ParticleInsertion<LinkedCell, dim>::Insertion) {
62 : // initialise velocity of molecule
63 0 : tarch::la::Vector<dim, double> velocity(0.0);
64 0 : _mdSolverInterface->getInitialVelocity(meanVelocity, _mdSolverInterface->getKB(), temperature, velocity);
65 0 : molecule.setVelocity(velocity);
66 :
67 : // add molecule to MD simulation and linked cell structures
68 0 : _mdSolverInterface->addMoleculeToMDSimulation(molecule);
69 : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
70 : std::cout << "Insert particle: Success " << molecule.getPosition() << std::endl;
71 : #endif
72 : }
73 : #if (COUPLING_MD_DEBUG == COUPLING_MD_YES)
74 : else {
75 : std::cout << "Insert particle: Failure" << std::endl;
76 : }
77 : #endif
78 0 : return action;
79 0 : }
80 :
81 : template <class LinkedCell, unsigned int dim>
82 : typename coupling::ParticleInsertion<LinkedCell, dim>::Action
83 0 : coupling::UsherParticleInsertion<LinkedCell, dim>::deleteParticle(coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& cell,
84 : const coupling::BoundaryForceController<LinkedCell, dim>& boundaryForceController) {
85 : // count particles with computeMassMapping
86 0 : coupling::cellmappings::ComputeMassMapping<LinkedCell, dim> computeMassMapping(_mdSolverInterface);
87 0 : cell.iterateConstCells(computeMassMapping);
88 :
89 : // take a random particle
90 0 : const unsigned int randomParticle =
91 0 : static_cast<unsigned int>(tarch::utils::RandomNumberService::getInstance().getUniformRandomNumber() * computeMassMapping.getNumberOfParticles());
92 :
93 : // delete this random particle
94 0 : coupling::cellmappings::DeleteParticleMapping<LinkedCell, dim> deleteParticleMapping(randomParticle, _mdSolverInterface);
95 0 : cell.iterateCells(deleteParticleMapping);
96 :
97 : #ifdef USHER_DEBUG
98 : _energyRemoved += deleteParticleMapping.getDeletedMolecule().getPotentialEnergy();
99 : _ZhouEnergyRemoved += boundaryForceController.getPotentialEnergy(deleteParticleMapping.getDeletedMolecule().getPosition());
100 : _particlesRemoved++;
101 : #endif
102 :
103 0 : return coupling::ParticleInsertion<LinkedCell, dim>::Deletion;
104 0 : }
105 :
106 : template <class LinkedCell, unsigned int dim>
107 0 : typename coupling::ParticleInsertion<LinkedCell, dim>::Action coupling::UsherParticleInsertion<LinkedCell, dim>::findParticlePosition(
108 : coupling::datastructures::CouplingCellWithLinkedCells<LinkedCell, dim>& thisCell, const tarch::la::Vector<dim, double>& couplingCellPosition,
109 : const tarch::la::Vector<dim, double>& couplingCellSize, coupling::datastructures::Molecule<dim>& molecule,
110 : const coupling::BoundaryForceController<LinkedCell, dim>& boundaryForceController) {
111 : // count particles with computeMassMapping
112 0 : coupling::cellmappings::ComputeMassMapping<LinkedCell, dim> computeMassMapping(_mdSolverInterface);
113 0 : thisCell.iterateConstCells(computeMassMapping);
114 : // compute number density; if it's zero, we assume to have at least one
115 : // particle in this cell (otherwise, we cannot choose an optimal step size for
116 : // USHER)
117 0 : double numberDensity = (double)computeMassMapping.getNumberOfParticles();
118 0 : if (numberDensity == 0.0) {
119 0 : numberDensity = 1.0;
120 : }
121 0 : for (unsigned int d = 0; d < dim; d++) {
122 0 : numberDensity = numberDensity / couplingCellSize[d];
123 : }
124 :
125 0 : tarch::la::Vector<dim, double> position(0.0);
126 0 : tarch::la::Vector<dim, double> positionOld(0.0);
127 0 : tarch::la::Vector<dim, double> force(0.0);
128 :
129 : // current energy
130 0 : double energy = 0.0;
131 : // absolute value of force
132 0 : double absForce = 0.0;
133 : // step size
134 0 : double stepSize = 0.0;
135 :
136 : // fluid parameters
137 0 : const double epsilon_times_4 = 4.0 * _mdSolverInterface->getMoleculeEpsilon();
138 0 : const double sigma = _mdSolverInterface->getMoleculeSigma();
139 :
140 : // -------- USHER parameters ------------
141 0 : const double rSigma = _usherParams._rSigmaCoeff * _mdSolverInterface->getMoleculeSigma();
142 : // energy level to be reached
143 0 : const double U_0 = _usherParams._meanPotentialEnergyFactor * thisCell.getPotentialEnergy();
144 : // energy implying a certain overlap with another particle (including
145 : // rescaling to LB scaling)
146 0 : const double U_overlap = _usherParams._uOverlapCoeff * _mdSolverInterface->getMoleculeEpsilon();
147 : // termination criterion (if relative energy |U-U0|/|U0| is smaller than
148 : // xiMax, search can be stopped)
149 0 : const double xiMax = _usherParams._tolerance;
150 : // maximum step size allowed
151 0 : const double stepRef = _usherParams.getStepRef(numberDensity, _mdSolverInterface->getMoleculeSigma());
152 : // number of particle movements allowed to find position with energy level U_0
153 0 : const int intIterMax = _usherParams._iterMax;
154 : // max. number of restart tries
155 0 : const int restartMax = _usherParams._restartMax;
156 :
157 : // upper right and lower left boundaries for particle insertion:
158 : // if we are close to the very outer boundary of the domain, we only allow
159 : // insertion within a distance of at least
160 : // _usherParams._offsetFromOuterBoundary. This is required since we may
161 : // encounter instabilities due to open boundary forcing which is not included
162 : // in the potential enery evaluation of the USHER scheme
163 0 : tarch::la::Vector<dim, double> upperRightBoundaries(couplingCellPosition + couplingCellSize);
164 0 : tarch::la::Vector<dim, double> lowerLeftBoundaries(couplingCellPosition);
165 0 : const tarch::la::Vector<dim, double> domainLower(_mdSolverInterface->getGlobalMDDomainOffset());
166 0 : const tarch::la::Vector<dim, double> domainUpper(_mdSolverInterface->getGlobalMDDomainSize() + domainLower);
167 0 : for (unsigned int d = 0; d < dim; d++) {
168 0 : upperRightBoundaries[d] = fmin(upperRightBoundaries[d], domainUpper[d] - _usherParams._offsetFromOuterBoundary);
169 0 : lowerLeftBoundaries[d] = fmax(lowerLeftBoundaries[d], domainLower[d] + _usherParams._offsetFromOuterBoundary);
170 : // if the offset yields, that we cannot insert any particle: return
171 0 : if (lowerLeftBoundaries[d] > upperRightBoundaries[d]) {
172 : return coupling::ParticleInsertion<LinkedCell, dim>::NoAction;
173 : }
174 : }
175 :
176 : #ifdef USHER_DEBUG
177 : std::cout << std::endl << "U_0 = " << U_0 << std::endl;
178 : #endif
179 :
180 : // try at max. restartMax times to insert this particle...
181 0 : for (int i = 0; i < restartMax; i++) {
182 :
183 : // generate random start position
184 0 : for (unsigned int d = 0; d < dim; d++) {
185 0 : position[d] = couplingCellPosition[d] + couplingCellSize[d] * tarch::utils::RandomNumberService::getInstance().getUniformRandomNumber();
186 : }
187 :
188 : // determine force and energy that act on molecule
189 0 : molecule.setPosition(position);
190 0 : _mdSolverInterface->calculateForceAndEnergy(molecule);
191 0 : energy = molecule.getPotentialEnergy();
192 0 : energy += boundaryForceController.getPotentialEnergy(position);
193 0 : force = molecule.getForce();
194 0 : force += boundaryForceController.getForce(position);
195 :
196 0 : if (energy - U_0 == 0.0) {
197 : #ifdef USHER_DEBUG
198 : std::cout << "energy-U_0 == 0.0" << std::endl;
199 : #endif
200 : return coupling::ParticleInsertion<LinkedCell, dim>::Insertion;
201 : }
202 : // determine signum
203 0 : int signAl = (int)((energy - U_0) / fabs(energy - U_0));
204 : #if (COUPLING_MD_ERROR == COUPLING_MD_YES)
205 0 : if ((signAl != 1) && (signAl != -1)) {
206 0 : std::cout << "ERROR coupling::UsherParticleInsertion::findParticlePosition(): "
207 : "wrong sign in USHER"
208 0 : << std::endl;
209 0 : exit(EXIT_FAILURE);
210 : }
211 : #endif
212 :
213 0 : absForce = std::sqrt(tarch::la::dot(force, force));
214 :
215 : // do steps towards expected energy level
216 0 : double xiloc = xiMax + 1.0;
217 0 : double xiOld = xiMax + 1.0;
218 0 : int success = 0;
219 0 : int loci = 0;
220 0 : for (; loci < intIterMax && success < 10; loci++) {
221 : // for checking, if a restart is required
222 0 : bool restartSearch = false;
223 :
224 : // if there is no force on the particle, we are in a low energy hole;
225 : // let's try to allow this
226 0 : if (absForce == 0.0 || energy - U_0 == 0.0) {
227 0 : molecule.setPosition(position);
228 0 : molecule.setForce(force);
229 0 : molecule.setPotentialEnergy(0.0);
230 : #ifdef USHER_DEBUG
231 : std::cout << "low energy hole" << std::endl;
232 : _particlesInserted++;
233 : #endif
234 0 : return coupling::ParticleInsertion<LinkedCell, dim>::Insertion;
235 : }
236 :
237 : // control step size
238 0 : if (energy > U_overlap) {
239 0 : stepSize = rSigma - sigma * pow(epsilon_times_4 / energy, (1.0 / 12.0));
240 :
241 : } else {
242 0 : stepSize = fabs(energy - U_0) / absForce;
243 0 : if (stepSize > stepRef) {
244 0 : stepSize = stepRef;
245 : }
246 : }
247 : #if (COUPLING_MD_ERROR == COUPLING_MD_YES)
248 0 : if (stepSize <= 0.0) {
249 0 : std::cout << "findParticlePosition(): ERROR "
250 : "coupling::UsherParticleInsertion::findParticlePosition():"
251 : " Stepsize is smaller than/ equal zero!"
252 0 : << std::endl;
253 0 : exit(EXIT_FAILURE);
254 : }
255 : #endif
256 :
257 : // update particle position
258 0 : positionOld = position;
259 0 : position = position + (stepSize * signAl / absForce) * force;
260 :
261 0 : molecule.setPosition(position);
262 :
263 : // restart searching if the new position of the particle is outside the
264 : // coupling cell
265 0 : for (unsigned int d = 0; d < dim; d++) {
266 0 : restartSearch = restartSearch || (position[d] >= upperRightBoundaries[d]) || (position[d] <= lowerLeftBoundaries[d]);
267 : }
268 0 : if (restartSearch) {
269 : break;
270 : }
271 :
272 : // determine force and energy that act on molecule
273 0 : _mdSolverInterface->calculateForceAndEnergy(molecule);
274 0 : energy = molecule.getPotentialEnergy();
275 0 : energy += boundaryForceController.getPotentialEnergy(position);
276 0 : force = molecule.getForce();
277 0 : force += boundaryForceController.getForce(position);
278 0 : signAl = (int)((energy - U_0) / fabs(energy - U_0));
279 :
280 0 : absForce = std::sqrt(tarch::la::dot(force, force));
281 0 : xiOld = xiloc;
282 0 : xiloc = fabs(energy - U_0) / fabs(U_0);
283 :
284 0 : if (xiloc < xiMax) {
285 0 : success++;
286 : }
287 :
288 : // restart searching if the difference energy-U_0 increases
289 0 : if (xiloc > xiOld) {
290 : break;
291 : }
292 : }
293 0 : if (success > 0) {
294 0 : molecule.setPosition(positionOld);
295 0 : _mdSolverInterface->calculateForceAndEnergy(molecule);
296 :
297 : #ifdef USHER_DEBUG
298 : _energyInserted += molecule.getPotentialEnergy();
299 : _ZhouEnergyInserted += boundaryForceController.getPotentialEnergy(positionOld);
300 : _particlesInserted++;
301 : energy = molecule.getPotentialEnergy();
302 : energy += boundaryForceController.getPotentialEnergy(positionOld);
303 : if (fabs(energy - U_0) / fabs(U_0) >= xiMax)
304 : std::cout << "USHER critical ERROR: fabs(energy-U_0)/fabs(U_0) >= xiMax" << std::endl;
305 : std::cout << "Finished with energy = " << energy << " after " << i << " restarts and " << loci << " iterations" << std::endl;
306 : #endif
307 :
308 0 : return coupling::ParticleInsertion<LinkedCell, dim>::Insertion;
309 : }
310 : }
311 : #ifdef USHER_DEBUG
312 : std::cout << "Failed" << std::endl;
313 : #endif
314 : return coupling::ParticleInsertion<LinkedCell, dim>::NoAction;
315 0 : }
|
