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 : /** index of linked cell (within linked cell domain, NOT within coupling
130 : * cell!) where molecule is placed in */
131 0 : tarch::la::Vector<dim, unsigned int> linkedCellIndex(0);
132 :
133 : // current energy
134 0 : double energy = 0.0;
135 : // absolute value of force
136 0 : double absForce = 0.0;
137 : // step size
138 0 : double stepSize = 0.0;
139 :
140 : // fluid parameters
141 0 : const double epsilon_times_4 = 4.0 * _mdSolverInterface->getMoleculeEpsilon();
142 0 : const double sigma = _mdSolverInterface->getMoleculeSigma();
143 :
144 : // -------- USHER parameters ------------
145 0 : const double rSigma = _usherParams._rSigmaCoeff * _mdSolverInterface->getMoleculeSigma();
146 : // energy level to be reached
147 0 : const double U_0 = _usherParams._meanPotentialEnergyFactor * thisCell.getPotentialEnergy();
148 : // energy implying a certain overlap with another particle (including
149 : // rescaling to LB scaling)
150 0 : const double U_overlap = _usherParams._uOverlapCoeff * _mdSolverInterface->getMoleculeEpsilon();
151 : // termination criterion (if relative energy |U-U0|/|U0| is smaller than
152 : // xiMax, search can be stopped)
153 0 : const double xiMax = _usherParams._tolerance;
154 : // maximum step size allowed
155 0 : const double stepRef = _usherParams.getStepRef(numberDensity, _mdSolverInterface->getMoleculeSigma());
156 : // number of particle movements allowed to find position with energy level U_0
157 0 : const int intIterMax = _usherParams._iterMax;
158 : // max. number of restart tries
159 0 : const int restartMax = _usherParams._restartMax;
160 :
161 : // upper right and lower left boundaries for particle insertion:
162 : // if we are close to the very outer boundary of the domain, we only allow
163 : // insertion within a distance of at least
164 : // _usherParams._offsetFromOuterBoundary. This is required since we may
165 : // encounter instabilities due to open boundary forcing which is not included
166 : // in the potential enery evaluation of the USHER scheme
167 0 : tarch::la::Vector<dim, double> upperRightBoundaries(couplingCellPosition + couplingCellSize);
168 0 : tarch::la::Vector<dim, double> lowerLeftBoundaries(couplingCellPosition);
169 0 : const tarch::la::Vector<dim, double> domainLower(_mdSolverInterface->getGlobalMDDomainOffset());
170 0 : const tarch::la::Vector<dim, double> domainUpper(_mdSolverInterface->getGlobalMDDomainSize() + domainLower);
171 0 : for (unsigned int d = 0; d < dim; d++) {
172 0 : upperRightBoundaries[d] = fmin(upperRightBoundaries[d], domainUpper[d] - _usherParams._offsetFromOuterBoundary);
173 0 : lowerLeftBoundaries[d] = fmax(lowerLeftBoundaries[d], domainLower[d] + _usherParams._offsetFromOuterBoundary);
174 : // if the offset yields, that we cannot insert any particle: return
175 0 : if (lowerLeftBoundaries[d] > upperRightBoundaries[d]) {
176 : return coupling::ParticleInsertion<LinkedCell, dim>::NoAction;
177 : }
178 : }
179 :
180 : #ifdef USHER_DEBUG
181 : std::cout << std::endl << "U_0 = " << U_0 << std::endl;
182 : #endif
183 :
184 : // try at max. restartMax times to insert this particle...
185 0 : for (int i = 0; i < restartMax; i++) {
186 :
187 : // generate random start position
188 0 : for (unsigned int d = 0; d < dim; d++) {
189 0 : position[d] = couplingCellPosition[d] + couplingCellSize[d] * tarch::utils::RandomNumberService::getInstance().getUniformRandomNumber();
190 : }
191 :
192 : // determine force and energy that act on molecule
193 0 : molecule.setPosition(position);
194 0 : _mdSolverInterface->calculateForceAndEnergy(molecule);
195 0 : energy = molecule.getPotentialEnergy();
196 0 : energy += boundaryForceController.getPotentialEnergy(position);
197 0 : force = molecule.getForce();
198 0 : force += boundaryForceController.getForce(position);
199 :
200 0 : if (energy - U_0 == 0.0) {
201 : #ifdef USHER_DEBUG
202 : std::cout << "energy-U_0 == 0.0" << std::endl;
203 : #endif
204 : return coupling::ParticleInsertion<LinkedCell, dim>::Insertion;
205 : }
206 : // determine signum
207 0 : int signAl = (int)((energy - U_0) / fabs(energy - U_0));
208 : #if (COUPLING_MD_ERROR == COUPLING_MD_YES)
209 0 : if ((signAl != 1) && (signAl != -1)) {
210 0 : std::cout << "ERROR coupling::UsherParticleInsertion::findParticlePosition(): "
211 : "wrong sign in USHER"
212 0 : << std::endl;
213 0 : exit(EXIT_FAILURE);
214 : }
215 : #endif
216 :
217 0 : absForce = std::sqrt(tarch::la::dot(force, force));
218 :
219 : // do steps towards expected energy level
220 0 : double xiloc = xiMax + 1.0;
221 0 : double xiOld = xiMax + 1.0;
222 0 : int success = 0;
223 0 : int loci = 0;
224 0 : for (; loci < intIterMax && success < 10; loci++) {
225 : // for checking, if a restart is required
226 0 : bool restartSearch = false;
227 :
228 : // if there is no force on the particle, we are in a low energy hole;
229 : // let's try to allow this
230 0 : if (absForce == 0.0 || energy - U_0 == 0.0) {
231 0 : molecule.setPosition(position);
232 0 : molecule.setForce(force);
233 0 : molecule.setPotentialEnergy(0.0);
234 : #ifdef USHER_DEBUG
235 : std::cout << "low energy hole" << std::endl;
236 : _particlesInserted++;
237 : #endif
238 0 : return coupling::ParticleInsertion<LinkedCell, dim>::Insertion;
239 : }
240 :
241 : // control step size
242 0 : if (energy > U_overlap) {
243 0 : stepSize = rSigma - sigma * pow(epsilon_times_4 / energy, (1.0 / 12.0));
244 :
245 : } else {
246 0 : stepSize = fabs(energy - U_0) / absForce;
247 0 : if (stepSize > stepRef) {
248 0 : stepSize = stepRef;
249 : }
250 : }
251 : #if (COUPLING_MD_ERROR == COUPLING_MD_YES)
252 0 : if (stepSize <= 0.0) {
253 0 : std::cout << "findParticlePosition(): ERROR "
254 : "coupling::UsherParticleInsertion::findParticlePosition():"
255 : " Stepsize is smaller than/ equal zero!"
256 0 : << std::endl;
257 0 : exit(EXIT_FAILURE);
258 : }
259 : #endif
260 :
261 : // update particle position
262 0 : positionOld = position;
263 0 : position = position + (stepSize * signAl / absForce) * force;
264 :
265 0 : molecule.setPosition(position);
266 :
267 : // restart searching if the new position of the particle is outside the
268 : // coupling cell
269 0 : for (unsigned int d = 0; d < dim; d++) {
270 0 : restartSearch = restartSearch || (position[d] >= upperRightBoundaries[d]) || (position[d] <= lowerLeftBoundaries[d]);
271 : }
272 0 : if (restartSearch) {
273 : break;
274 : }
275 :
276 : // determine force and energy that act on molecule
277 0 : _mdSolverInterface->calculateForceAndEnergy(molecule);
278 0 : energy = molecule.getPotentialEnergy();
279 0 : energy += boundaryForceController.getPotentialEnergy(position);
280 0 : force = molecule.getForce();
281 0 : force += boundaryForceController.getForce(position);
282 0 : signAl = (int)((energy - U_0) / fabs(energy - U_0));
283 :
284 0 : absForce = std::sqrt(tarch::la::dot(force, force));
285 0 : xiOld = xiloc;
286 0 : xiloc = fabs(energy - U_0) / fabs(U_0);
287 :
288 0 : if (xiloc < xiMax) {
289 0 : success++;
290 : }
291 :
292 : // restart searching if the difference energy-U_0 increases
293 0 : if (xiloc > xiOld) {
294 : break;
295 : }
296 : }
297 0 : if (success > 0) {
298 0 : molecule.setPosition(positionOld);
299 0 : _mdSolverInterface->calculateForceAndEnergy(molecule);
300 :
301 : #ifdef USHER_DEBUG
302 : _energyInserted += molecule.getPotentialEnergy();
303 : _ZhouEnergyInserted += boundaryForceController.getPotentialEnergy(positionOld);
304 : _particlesInserted++;
305 : energy = molecule.getPotentialEnergy();
306 : energy += boundaryForceController.getPotentialEnergy(positionOld);
307 : if (fabs(energy - U_0) / fabs(U_0) >= xiMax)
308 : std::cout << "USHER critical ERROR: fabs(energy-U_0)/fabs(U_0) >= xiMax" << std::endl;
309 : std::cout << "Finished with energy = " << energy << " after " << i << " restarts and " << loci << " iterations" << std::endl;
310 : #endif
311 :
312 0 : return coupling::ParticleInsertion<LinkedCell, dim>::Insertion;
313 : }
314 : }
315 : #ifdef USHER_DEBUG
316 : std::cout << "Failed" << std::endl;
317 : #endif
318 : return coupling::ParticleInsertion<LinkedCell, dim>::NoAction;
319 0 : }
|
