MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HGridOptimiser.cc
Go to the documentation of this file.
1 //Copyright (c) 2013-2020, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 #include <iostream>
27 #include <iomanip>
28 #include <limits>
29 #include <cmath>
30 #include "HGridOptimiser.h"
31 #include "Particles/BaseParticle.h"
32 
33 void HGridOptimiser::initialise(const MercuryBase& problem, unsigned int numberOfCells, int verbosity)
34 {
35  //assign the number of cells
36  numCells_ = numberOfCells;
38 
39  cellN_.resize(numCells_);
40  for (double& it : cellN_)
41  {
42  it = 0.0;
43  }
44 
45  //Set the minimum and maximum radius of the particles.
47  rMax_ = nextafter(nextafter(problem.particleHandler.getLargestParticle()->getMaxInteractionRadius(),
48  std::numeric_limits<double>::max()), std::numeric_limits<double>::max());
49 
50  //for all particles, add one to the cell in cellN_ which has is associated with
51  //the radius of that particle.
52  for (std::vector<BaseParticle*>::const_iterator it = problem.particleHandler.begin();
53  it != problem.particleHandler.end(); ++it)
54  {
55  cellN_[radius2Cell((*it)->getMaxInteractionRadius())]++;
56  }
57 
58  //assign what the number of particles is.
59  unsigned int numParticles = problem.particleHandler.getSize();
60 
61  intCellN.push_back(1.5 * cellN_[0] - 0.5 * cellN_[1]);
62  for (double& it : cellN_)
63  {
64  intCellN.push_back(it);
65  }
66  intCellN.push_back(1.5 * cellN_[numCells_ - 1] - 0.5 * cellN_[numCells_ - 2]);
67 
68  if (verbosity > 0)
69  {
70  std::cout << "rMin=" << rMin_ << " rMax=" << rMax_ << " NParticles=" << numParticles << std::endl;
71  for (unsigned int i = 0; i < cellN_.size(); i++)
72  {
73  std::cout << "From" << std::setw(8) << cell2Min(i) << " to" << std::setw(8) << cell2Max(i) << " number="
74  << std::setw(5) << cellN_[i] << std::endl;
75  }
76  std::cout << "Domain size [" << problem.getXMin() << "," << problem.getXMax() << "]x[" << problem.getYMin()
77  << "," << problem.getYMax() << "]x[" << problem.getZMin() << "," << problem.getZMax() << "]"
78  << std::endl;
79 
80  for (unsigned int i = 0; i < intCellN.size() - 1; i++)
81  {
82  std::cout << "From" << std::setw(8) << intCell2Min(i) << " to" << std::setw(8) << intCell2Max(i)
83  << " linear from " << std::setw(8) << intCellN[i] << " to " << std::setw(8) << intCellN[i + 1]
84  << std::endl;
85  }
86  /*std::cout<<"["<<cellN_[0];
87  for (unsigned int i=1;i<cellN_.size();i++)
88  {
89  std::cout<<","<<cellN_[i];
90  }
91  std::cout<<"]"<<std::endl;*/
92  }
93 
94  dimension_ = problem.getSystemDimensions();
95  switch (dimension_)
96  {
97  case 1:
98  length_ = pow((problem.getXMax() - problem.getXMin()) / numParticles, 1);
99  break;
100  case 2:
101  length_ = pow(
102  (problem.getXMax() - problem.getXMin()) * (problem.getYMax() - problem.getYMin()) / numParticles,
103  1.0 / 2.0);
104  break;
105  case 3:
106  length_ = pow((problem.getXMax() - problem.getXMin()) * (problem.getYMax() - problem.getYMin()) *
107  (problem.getZMax() - problem.getZMin()) / numParticles, 1.0 / 3.0);
108  break;
109  }
110 }
111 
112 void
113 HGridOptimiser::initialisePolyFunc(double omega, std::vector<double>& coeff, unsigned int numberOfCells, int verbosity)
114 {
115  numCells_ = numberOfCells;
117  cellN_.resize(numCells_);
118  for (double& it : cellN_)
119  {
120  it = 0;
121  }
122 
123  rMin_ = nextafter(1, 0.0);
124  rMax_ = nextafter(omega, std::numeric_limits<double>::max());
125  for (unsigned int i = 0; i < cellN_.size(); i++)
126  {
127  double start = cell2Min(i);
128  double end = cell2Max(i);
129  for (unsigned int j = 0; j < coeff.size(); j++)
130  {
131  cellN_[i] += coeff[j] / (1.0 + j) * (std::pow(end, j + 1) - std::pow(start, j + 1));
132  }
133  }
134 
135  intCellN.push_back(1.5 * cellN_[0] - 0.5 * cellN_[1]);
136  for (double& it : cellN_)
137  {
138  intCellN.push_back(it);
139  }
140  intCellN.push_back(1.5 * cellN_[numCells_ - 1] - 0.5 * cellN_[numCells_ - 2]);
141 
142  dimension_ = 1;
143  length_ = 1.0;
144 
145  if (verbosity > 0)
146  {
147  std::cout << "rMin=" << rMin_ << " rMax=" << rMax_ << std::endl;
148  for (unsigned int i = 0; i < cellN_.size(); i++)
149  {
150  std::cout << "From" << std::setw(8) << cell2Min(i) << " to" << std::setw(8) << cell2Max(i) << " number="
151  << std::setw(5) << cellN_[i] << std::endl;
152  }
153 
154  for (unsigned int i = 0; i < intCellN.size() - 1; i++)
155  {
156  std::cout << "From" << std::setw(8) << intCell2Min(i) << " to" << std::setw(8) << intCell2Max(i)
157  << " linear from " << std::setw(8) << intCellN[i] << " to " << std::setw(8) << intCellN[i + 1]
158  << std::endl;
159  }
160  }
161 }
162 
170 unsigned int HGridOptimiser::radius2Cell(double r)
171 {
172  if (r < rMin_ || r >= rMax_)
173  logger(ERROR, "The radius % is not in the range [%,%]", r, rMin_, rMax_);
174  unsigned int y = static_cast<unsigned int> (floor(numCells_ * (r - rMin_) / (rMax_ - rMin_)));
175  return y;
176 }
177 
178 unsigned int HGridOptimiser::radius2IntCell(double r)
179 {
180  if (r < rMin_ || r >= rMax_)
181  throw;
182  unsigned int y = static_cast<unsigned int> (floor(numCells_ * (r - rMin_) / (rMax_ - rMin_) + 0.5));
183  return y;
184 }
185 
186 double HGridOptimiser::intCell2Min(unsigned int i)
187 {
188  if (i == 0)
189  return rMin_;
190  else
191  return rMin_ + (rMax_ - rMin_) * (i - 0.5) / numCells_;
192 }
193 
194 double HGridOptimiser::intCell2Max(unsigned int i)
195 {
196  if (i == numCells_)
197  return rMax_;
198  else
199  return rMin_ + (rMax_ - rMin_) * (i + 0.5) / numCells_;
200 }
201 
205 double HGridOptimiser::cell2Min(unsigned int i)
206 {
207  return rMin_ + (rMax_ - rMin_) * i / numCells_;
208 }
209 
213 double HGridOptimiser::cell2Max(unsigned int i)
214 {
215  return rMin_ + (rMax_ - rMin_) * (i + 1) / numCells_;
216 }
217 
218 double HGridOptimiser::pdfIntCell(double start, double end, unsigned int i, int p)
219 {
220  double a = (intCellN[i + 1] - intCellN[i]) / (intCell2Max(i) - intCell2Min(i));
221  double b = (intCellN[i] * intCell2Max(i) - intCellN[i + 1] * intCell2Min(i)) / (intCell2Max(i) - intCell2Min(i));
222  return (a / (p + 2) * (std::pow(end, p + 2) - std::pow(start, p + 2)) +
223  b / (p + 1) * (std::pow(end, p + 1) - std::pow(start, p + 1)));
224 }
225 
229 double HGridOptimiser::pdfInt(double start, double end, int p)
230 {
231  //std::cout<<"pdfInit("<<start<<","<<end<<","<<p<<");"<<std::endl;
232  unsigned int startCell = radius2IntCell(start);
233  unsigned int endCell = radius2IntCell(end);
234 
235  double numerator = 0;
236  if (startCell == endCell)
237  {
238  numerator = pdfIntCell(start, end, startCell, p);
239  }
240  else if (endCell < startCell)
241  {
242  numerator = 0;
243  }
244  else
245  {
246  numerator = pdfIntCell(start, intCell2Max(startCell), startCell, p);
247  //std::cout<<"Teller+="<<pdfIntCell(start,intCell2Max(startCell),startCell,p)<<std::endl;
248  for (unsigned int i = startCell + 1; i < endCell; i++)
249  {
250  numerator += pdfIntCell(intCell2Min(i), intCell2Max(i), i, p);
251  //std::cout<<"Teller+="<<pdfIntCell(intCell2Min(i),intCell2Max(i),i,p)<<std::endl;
252  }
253  numerator += pdfIntCell(intCell2Min(endCell), end, endCell, p);
254  //std::cout<<"Teller+="<<pdfIntCell(intCell2Min(endCell),end,endCell,p)<<std::endl;
255  }
256  double denominator = 0;
257  for (unsigned int i = 0; i <= numCells_; i++)
258  {
259  denominator += pdfIntCell(intCell2Min(i), intCell2Max(i), i, 0);
260  //std::cout<<"Noemer+="<<pdfIntCell(intCell2Min(i),intCell2Max(i),i,0)<<std::endl;
261  }
262  return numerator / denominator;
263 }
264 
266 double HGridOptimiser::diffPdfInt(double x, int p)
267 {
268  unsigned int xCell = radius2IntCell(x);
269  double denominator = 0;
270  for (unsigned int i = 0; i <= numCells_; i++)
271  {
272  denominator += pdfIntCell(intCell2Min(i), intCell2Max(i), i, 0);
273  }
274  double a = (intCellN[xCell + 1] - intCellN[xCell]) / (intCell2Max(xCell) - intCell2Min(xCell));
275  double b = (intCellN[xCell] * intCell2Max(xCell) - intCellN[xCell + 1] * intCell2Min(xCell)) /
276  (intCell2Max(xCell) - intCell2Min(xCell));
277  return (a * std::pow(x, p + 1) + b * std::pow(x, p)) / denominator;
278 }
279 
280 double HGridOptimiser::expectedCellsIntegralCellNumerator(double start, double end, unsigned int i, int p, double h)
281 {
282  double a = (intCellN[i + 1] - intCellN[i]) / (intCell2Max(i) - intCell2Min(i));
283  double b = (intCellN[i] * intCell2Max(i) - intCellN[i + 1] * intCell2Min(i)) / (intCell2Max(i) - intCell2Min(i));
284  double r;
285  r = start;
286  double min = std::pow(2.0 * (r + h) / h, p) * (a * p * r - a * h + a * r + b * p + 2 * b) * (r + h) /
287  ((p + 1) * (p + 2));
288  r = end;
289  double plus = std::pow(2.0 * (r + h) / h, p) * (a * p * r - a * h + a * r + b * p + 2 * b) * (r + h) /
290  ((p + 1) * (p + 2));
291  return plus - min;
292 }
293 
294 double
295 HGridOptimiser::diffHExpectedCellsIntegralCellNumerator(double start, double end, unsigned int i, int p, double h)
296 {
297  double a = (intCellN[i + 1] - intCellN[i]) / (intCell2Max(i) - intCell2Min(i));
298  double b = (intCellN[i] * intCell2Max(i) - intCellN[i + 1] * intCell2Min(i)) / (intCell2Max(i) - intCell2Min(i));
299  double r;
300  r = start;
301  //double min =std::pow(2.0*r/h+2.0,p)*(a*p*r-a*h+a*r+b*p+2*b)*(r+h)/((p+1)*(p+2));
302  double min =
303  (-2.0 * r / h / h * p * std::pow(2.0 * r / h + 2.0, p - 1) * (a * p * r - a * h + a * r + b * p + 2 * b) *
304  (r + h) +
305  -std::pow(2.0 * r / h + 2.0, p) * a * (r + h) +
306  std::pow(2.0 * r / h + 2.0, p) * (a * p * r - a * h + a * r + b * p + 2 * b)) / ((p + 1) * (p + 2));
307  r = end;
308  //double plus=std::pow(2.0*r/h+2.0,p)*(a*p*r-a*h+a*r+b*p+2*b)*(r+h)/((p+1)*(p+2));
309  double plus =
310  (-2.0 * r / h / h * p * std::pow(2.0 * r / h + 2.0, p - 1) * (a * p * r - a * h + a * r + b * p + 2 * b) *
311  (r + h) +
312  -std::pow(2.0 * r / h + 2.0, p) * a * (r + h) +
313  std::pow(2.0 * r / h + 2.0, p) * (a * p * r - a * h + a * r + b * p + 2 * b)) / ((p + 1) * (p + 2));
314  return plus - min;
315 }
316 
317 double HGridOptimiser::expectedCellsIntegralCellDenominator(double start, double end, unsigned int i)
318 {
319  double a = (intCellN[i + 1] - intCellN[i]) / (intCell2Max(i) - intCell2Min(i));
320  double b = (intCellN[i] * intCell2Max(i) - intCellN[i + 1] * intCell2Min(i)) / (intCell2Max(i) - intCell2Min(i));
321  return (a / 2.0 * (pow(end, 2) - pow(start, 2)) + b * (pow(end, 1) - pow(start, 1)));
322 }
323 
328 double HGridOptimiser::expectedCellsIntegral(double start, double end, int p, double h)
329 {
330  unsigned int startCell = radius2IntCell(start);
331  unsigned int endCell = radius2IntCell(end);
332 
333  double numerator = 0;
334  double denominator = 0;
335  if (startCell == endCell)
336  {
337  numerator = expectedCellsIntegralCellNumerator(start, end, startCell, p, h);
338  denominator = expectedCellsIntegralCellDenominator(start, end, startCell);
339  }
340  else
341  {
342  numerator = expectedCellsIntegralCellNumerator(start, intCell2Max(startCell), startCell, p, h);
343  denominator = expectedCellsIntegralCellDenominator(start, intCell2Max(startCell), startCell);
344  for (unsigned int i = startCell + 1; i < endCell; i++)
345  {
348  }
349  numerator += expectedCellsIntegralCellNumerator(intCell2Min(endCell), end, endCell, p, h);
350  denominator += expectedCellsIntegralCellDenominator(intCell2Min(endCell), end, endCell);
351 
352  }
353  // std::cout<<"New: numerator="<<numerator<<" denominator="<<denominator<<std::endl;
354  return numerator / denominator;
355 }
356 
357 double
358 HGridOptimiser::diffStartExpectedCellsIntegral(double start, double end, int p, double h)
359 {
360  unsigned int startCell = radius2IntCell(start);
361  unsigned int endCell = radius2IntCell(end);
362 
363  double numerator = 0;
364  double denominator = 0;
365 
366  double a = (intCellN[startCell + 1] - intCellN[startCell]) / (intCell2Max(startCell) - intCell2Min(startCell));
367  double b = (intCellN[startCell] * intCell2Max(startCell) - intCellN[startCell + 1] * intCell2Min(startCell)) /
368  (intCell2Max(startCell) - intCell2Min(startCell));
369 
370  double diffTeller = -(
371  2.0 / h * p * std::pow(2.0 * (start + h) / h, p - 1) * (a * p * start - a * h + a * start + b * p + 2 * b) *
372  (start + h) / ((p + 1) * (p + 2)) +
373  std::pow(2.0 * (start + h) / h, p) * (a * p + a) * (start + h) / ((p + 1) * (p + 2)) +
374  std::pow(2.0 * (start + h) / h, p) * (a * p * start - a * h + a * start + b * p + 2 * b) /
375  ((p + 1) * (p + 2)));
376  double diffNoemer = -(a * start + b);
377  if (startCell == endCell)
378  {
379  numerator = expectedCellsIntegralCellNumerator(start, end, startCell, p, h);
380  denominator = expectedCellsIntegralCellDenominator(start, end, startCell);
381  }
382  else
383  {
384  numerator = expectedCellsIntegralCellNumerator(start, intCell2Max(startCell), startCell, p, h);
385  denominator = expectedCellsIntegralCellDenominator(start, intCell2Max(startCell), startCell);
386  for (unsigned int i = startCell + 1; i < endCell; i++)
387  {
390  }
391  numerator += expectedCellsIntegralCellNumerator(intCell2Min(endCell), end, endCell, p, h);
392  denominator += expectedCellsIntegralCellDenominator(intCell2Min(endCell), end, endCell);
393 
394  }
395  //std::cout<<"new: numerator="<<numerator<<" denominator="<<denominator<<" diffTeller="<<diffTeller<<" diffNoemer="<<diffNoemer<<std::endl;
396  return (diffTeller * denominator - numerator * diffNoemer) / (denominator * denominator);
397 }
398 
399 double HGridOptimiser::diffEndExpectedCellsIntegral(double start, double end, int p, double h)
400 {
401  unsigned int startCell = radius2IntCell(start);
402  unsigned int endCell = radius2IntCell(end);
403 
404  double numerator = 0;
405  double denominator = 0;
406 
407  double a = (intCellN[endCell + 1] - intCellN[endCell]) / (intCell2Max(endCell) - intCell2Min(endCell));
408  double b = (intCellN[endCell] * intCell2Max(endCell) - intCellN[endCell + 1] * intCell2Min(endCell)) /
409  (intCell2Max(endCell) - intCell2Min(endCell));
410 
411  double diffTeller = (
412  2.0 / h * p * std::pow(2.0 * (end + h) / h, p - 1) * (a * p * end - a * h + a * end + b * p + 2 * b) *
413  (end + h) / ((p + 1) * (p + 2)) +
414  std::pow(2.0 * (end + h) / h, p) * (a * p + a) * (end + h) / ((p + 1) * (p + 2)) +
415  std::pow(2.0 * (end + h) / h, p) * (a * p * end - a * h + a * end + b * p + 2 * b) / ((p + 1) * (p + 2)));
416  double diffNoemer = (a * end + b);
417  if (startCell == endCell)
418  {
419  numerator = expectedCellsIntegralCellNumerator(start, end, startCell, p, h);
420  denominator = expectedCellsIntegralCellDenominator(start, end, startCell);
421  }
422  else
423  {
424  numerator = expectedCellsIntegralCellNumerator(start, intCell2Max(startCell), startCell, p, h);
425  denominator = expectedCellsIntegralCellDenominator(start, intCell2Max(startCell), startCell);
426  for (unsigned int i = startCell + 1; i < endCell; i++)
427  {
430  }
431  numerator += expectedCellsIntegralCellNumerator(intCell2Min(endCell), end, endCell, p, h);
432  denominator += expectedCellsIntegralCellDenominator(intCell2Min(endCell), end, endCell);
433 
434  }
435  //std::cout<<"new: numerator="<<numerator<<" denominator="<<denominator<<" diffTeller="<<diffTeller<<" diffNoemer="<<diffNoemer<<std::endl;
436  return (diffTeller * denominator - numerator * diffNoemer) / (denominator * denominator);
437 }
438 
439 double HGridOptimiser::diffHExpectedCellsIntegral(double start, double end, int p, double h)
440 {
441  unsigned int startCell = radius2IntCell(start);
442  unsigned int endCell = radius2IntCell(end);
443 
444  double denominator = 0;
445  double diffTeller = 0;
446 
447  if (startCell == endCell)
448  {
449  diffTeller = diffHExpectedCellsIntegralCellNumerator(start, end, startCell, p, h);
450  denominator = expectedCellsIntegralCellDenominator(start, end, startCell);
451  }
452  else
453  {
454  diffTeller = diffHExpectedCellsIntegralCellNumerator(start, intCell2Max(startCell), startCell, p, h);
455  denominator = expectedCellsIntegralCellDenominator(start, intCell2Max(startCell), startCell);
456  for (unsigned int i = startCell + 1; i < endCell; i++)
457  {
460  }
461  diffTeller += diffHExpectedCellsIntegralCellNumerator(intCell2Min(endCell), end, endCell, p, h);
462  denominator += expectedCellsIntegralCellDenominator(intCell2Min(endCell), end, endCell);
463  }
464  return diffTeller / denominator;
465 }
466 
467 void
468 HGridOptimiser::calculateDiffWork(std::vector<double>& hGridCellSizes, std::vector<double>& dfdx, HGridMethod method,
469  int verbosity)
470 {
471  unsigned int NLevels = hGridCellSizes.size() - 1;
472  unsigned int NParams = hGridCellSizes.size();
473 
474  if (verbosity > 0)
475  {
476  std::cout << "Length scale=" << length_ << std::endl;
477  }
478  if (verbosity > 0)
479  {
480  std::cout << "Level sizes:";
481  for (unsigned int i = 0; i < NParams; i++)
482  {
483  std::cout << " " << hGridCellSizes[i];
484  }
485  std::cout << std::endl;
486  }
487 
488  std::vector<double> particlesPerLevel;
489  for (unsigned int i = 0; i < NLevels; i++)
490  {
491  particlesPerLevel.push_back(pdfInt(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], 0));
492  if (particlesPerLevel.back() < 0)
493  {
494  particlesPerLevel.back() = 0;
495  }
496  }
497  if (verbosity > 0)
498  {
499  std::cout << "Particles per level: ";
500  for (unsigned int i = 0; i < NLevels; i++)
501  {
502  std::cout << " " << particlesPerLevel[i];
503  }
504  std::cout << std::endl;
505  }
506 
507  //How changes particlesPerLeve[i] when hGridCellSize[j] is changed
508  std::vector<std::vector<double> > diffParticlesPerLevel;
509  diffParticlesPerLevel.resize(NLevels);
510  for (unsigned int i = 0; i < NLevels; i++)
511  {
512  for (unsigned int j = 0; j < NParams; j++)
513  {
514  diffParticlesPerLevel[i].push_back(0.0);
515  if (j == i)
516  {
517  diffParticlesPerLevel[i][j] += -0.5 * diffPdfInt(0.5 * hGridCellSizes[i], 0);
518  }
519  if (j == i + 1)
520  {
521  diffParticlesPerLevel[i][j] += 0.5 * diffPdfInt(0.5 * hGridCellSizes[i + 1], 0);
522  }
523  }
524  }
525  if (verbosity > 0)
526  {
527  std::cout << "Diff Particles per level: " << std::endl;
528  for (unsigned int i = 0; i < NLevels; i++)
529  {
530  for (unsigned int j = 0; j < NParams; j++)
531  {
532  std::cout << " " << std::setw(12) << diffParticlesPerLevel[i][j];
533  }
534  std::cout << std::endl;
535  }
536  }
537 
538  std::vector<double> cellsPerLevel;
539  for (unsigned int i = 0; i < NLevels; i++)
540  {
541  cellsPerLevel.push_back(pow(length_ / hGridCellSizes[i + 1], dimension_));
542  }
543  if (verbosity > 0)
544  {
545  std::cout << "Cells per level: ";
546  for (unsigned int i = 0; i < NLevels; i++)
547  {
548  std::cout << " " << cellsPerLevel[i];
549  }
550  std::cout << std::endl;
551  }
552 
553  //How changes cellsPerLevel[i] when hGridCellSize[j] is changed
554  std::vector<std::vector<double> > diffCellsPerLevel;
555  diffCellsPerLevel.resize(hGridCellSizes.size());
556  for (unsigned int i = 0; i < NLevels; i++)
557  {
558  for (unsigned int j = 0; j < NParams; j++)
559  {
560  if (j == i + 1)
561  {
562  diffCellsPerLevel[i].push_back(
563  -pow(length_ / hGridCellSizes[i + 1], dimension_) * dimension_ / hGridCellSizes[i + 1]);
564  }
565  else
566  {
567  diffCellsPerLevel[i].push_back(0.0);
568  }
569  }
570  }
571  if (verbosity > 0)
572  {
573  std::cout << "Diff Cells per level: " << std::endl;
574  for (unsigned int i = 0; i < NLevels; i++)
575  {
576  for (unsigned int j = 0; j < NParams; j++)
577  {
578  std::cout << " " << std::setw(12) << diffCellsPerLevel[i][j];
579  }
580  std::cout << std::endl;
581  }
582  }
583 
584  std::vector<double> particlesPerCell;
585  for (unsigned int i = 0; i < NLevels; i++)
586  {
587  particlesPerCell.push_back(particlesPerLevel[i] / cellsPerLevel[i]);
588  }
589  if (verbosity > 0)
590  {
591  std::cout << "Particles per cell: ";
592  for (double i : particlesPerCell)
593  {
594  std::cout << " " << i;
595  }
596  std::cout << std::endl;
597  }
598 
599  //How changes particlesPerCell[i] when hGridCellSize[j] is changed
600  std::vector<std::vector<double> > diffParticlesPerCell;
601  diffParticlesPerCell.resize(NLevels);
602  for (unsigned int i = 0; i < NLevels; i++)
603  {
604  for (unsigned int j = 0; j < NParams; j++)
605  {
606  diffParticlesPerCell[i].push_back(
607  (diffParticlesPerLevel[i][j] * cellsPerLevel[i] - particlesPerLevel[i] * diffCellsPerLevel[i][j]) /
608  (cellsPerLevel[i] * cellsPerLevel[i]));
609  }
610  }
611  if (verbosity > 0)
612  {
613  std::cout << "Diff Particles per Cell: " << std::endl;
614  for (unsigned int i = 0; i < NLevels; i++)
615  {
616  for (unsigned int j = 0; j < NParams; j++)
617  {
618  std::cout << " " << std::setw(12) << diffParticlesPerCell[i][j];
619  }
620  std::cout << std::endl;
621  }
622  }
623 
624  double contactWork = 0, overheadWork = 0;
625  std::vector<double> diffContactWork;
626  std::vector<double> diffOverheadWork;
627  diffContactWork.resize(NParams);
628  diffOverheadWork.resize(NParams);
629  for (unsigned int j = 0; j < NParams; j++)
630  {
631  diffContactWork[j] = 0;
632  diffOverheadWork[j] = 0;
633  }
634 
635  for (unsigned int i = 0; i < NLevels; i++)
636  {
637  contactWork += 0.5 * pow(3, dimension_) * particlesPerCell[i] * particlesPerLevel[i];
638  overheadWork += 0.5 * (pow(3, dimension_) + 1) * particlesPerLevel[i];
639  for (unsigned int j = 0; j < NParams; j++)
640  {
641  diffContactWork[j] += 0.5 * pow(3, dimension_) * (diffParticlesPerCell[i][j] * particlesPerLevel[i] +
642  particlesPerCell[i] * diffParticlesPerLevel[i][j]);
643  diffOverheadWork[j] += 0.5 * (pow(3, dimension_) + 1) * diffParticlesPerLevel[i][j];
644  }
645 
646  unsigned int startJ, endJ;
647  switch (method)
648  {
649  case BOTTOMUP:
650  {
651  startJ = i + 1;
652  endJ = NLevels;
653  break;
654  }
655  case TOPDOWN:
656  {
657  startJ = 0;
658  endJ = i;
659  break;
660  }
661  }
662  for (unsigned int j = startJ; j < endJ; j++)
663  {
664  double numberOfCellToCheck;
665  numberOfCellToCheck = expectedCellsIntegral(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1],
666  dimension_, hGridCellSizes[j + 1]);
667 
668  std::vector<double> diffNumberOfCellToCheck;
669  for (unsigned int k = 0; k < NParams; k++)
670  {
671  diffNumberOfCellToCheck.push_back(0.0);
672  if (k == i)
673  {
674  diffNumberOfCellToCheck[k] += 0.5 * diffStartExpectedCellsIntegral(0.5 * hGridCellSizes[i],
675  0.5 * hGridCellSizes[i + 1],
676  dimension_,
677  hGridCellSizes[j + 1]);
678  }
679  if (k == i + 1)
680  {
681  diffNumberOfCellToCheck[k] += 0.5 * diffEndExpectedCellsIntegral(0.5 * hGridCellSizes[i],
682  0.5 * hGridCellSizes[i + 1],
683  dimension_, hGridCellSizes[j + 1]);
684  }
685  if (k == j + 1)
686  {
687  diffNumberOfCellToCheck[k] += diffHExpectedCellsIntegral(0.5 * hGridCellSizes[i],
688  0.5 * hGridCellSizes[i + 1], dimension_,
689  hGridCellSizes[j + 1]);
690  }
691  }
692  contactWork += particlesPerLevel[i] * numberOfCellToCheck * particlesPerCell[j];
693  overheadWork += particlesPerLevel[i] * numberOfCellToCheck;
694  //std::cout<<"i= "<<i<<" j= "<<j<<" NumberOfCellToCheck= "<<numberOfCellToCheck<<" diffNumberOfCellToCheck[2]= "<<diffNumberOfCellToCheck[2]<<std::endl;
695  for (unsigned int k = 0; k < NParams; k++)
696  {
697  diffContactWork[k] += (diffParticlesPerLevel[i][k] * numberOfCellToCheck * particlesPerCell[j] +
698  particlesPerLevel[i] * diffNumberOfCellToCheck[k] * particlesPerCell[j] +
699  particlesPerLevel[i] * numberOfCellToCheck * diffParticlesPerCell[j][k]);
700  diffOverheadWork[k] += (diffParticlesPerLevel[i][k] * numberOfCellToCheck +
701  particlesPerLevel[i] * diffNumberOfCellToCheck[k]);
702  }
703  }
704  }
705  if (verbosity > 0)
706  {
707  std::cout << "Contact work: " << contactWork << std::endl;
708  std::cout << "Overhead work: " << cellCheckOverContactCheckRatio_ * overheadWork << std::endl;
709  std::cout << "Sum work: " << contactWork + cellCheckOverContactCheckRatio_ * overheadWork << std::endl;
710  std::cout << "diff Contactwork: ";
711  for (unsigned int j = 0; j < NParams; j++)
712  {
713  std::cout << " " << diffContactWork[j];
714  }
715  std::cout << std::endl;
716  std::cout << "diff Overheadwork: ";
717  for (unsigned int j = 0; j < NParams; j++)
718  {
719  std::cout << " " << cellCheckOverContactCheckRatio_ * diffOverheadWork[j];
720  }
721  std::cout << std::endl;
722  std::cout << "diff sum: ";
723  for (unsigned int j = 0; j < NParams; j++)
724  {
725  std::cout << " " << diffContactWork[j] + cellCheckOverContactCheckRatio_ * diffOverheadWork[j];
726  }
727  std::cout << std::endl;
728  }
729  dfdx.clear();
730  for (unsigned int j = 0; j < NParams; j++)
731  {
732  dfdx.push_back(diffContactWork[j] + cellCheckOverContactCheckRatio_ * diffOverheadWork[j]);
733  }
734 }
735 
736 double HGridOptimiser::calculateWork(std::vector<double>& hGridCellSizes, HGridMethod method, int verbosity)
737 {
738  unsigned int NLevels = hGridCellSizes.size() - 1;
739  unsigned int NParams = hGridCellSizes.size();
740 
741  if (verbosity > 1)
742  {
743  std::cout << "Length scale=" << length_ << std::endl;
744  }
745  if (verbosity > 0)
746  {
747  std::cout << "Level sizes:" << std::endl;
748  for (unsigned int i = 0; i < NParams; i++)
749  {
750  std::cout << std::setw(10) << hGridCellSizes[i] << " ";
751  }
752  std::cout << std::endl;
753  }
754 
755  std::vector<double> particlesPerLevel;
756  for (unsigned int i = 0; i < NLevels; i++)
757  {
758  particlesPerLevel.push_back(pdfInt(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], 0));
759  if (particlesPerLevel.back() < 0)
760  {
761  particlesPerLevel.back() = 0;
762  }
763  }
764  if (verbosity > 0)
765  {
766  std::cout << "Particles per level:" << std::endl;
767  for (unsigned int i = 0; i < NLevels; i++)
768  {
769  std::cout << std::setw(10) << particlesPerLevel[i] << " ";
770  }
771  std::cout << std::endl;
772  }
773 
774  std::vector<double> cellsPerLevel;
775  for (unsigned int i = 0; i < NLevels; i++)
776  {
777  cellsPerLevel.push_back(pow(length_ / hGridCellSizes[i + 1], dimension_));
778  }
779  if (verbosity > 1)
780  {
781  std::cout << "Cells per level:" << std::endl;
782  for (unsigned int i = 0; i < NLevels; i++)
783  {
784  std::cout << std::setw(10) << cellsPerLevel[i] << " ";
785  }
786  std::cout << std::endl;
787  }
788 
789  std::vector<double> particlesPerCell;
790  for (unsigned int i = 0; i < NLevels; i++)
791  {
792  particlesPerCell.push_back(particlesPerLevel[i] / cellsPerLevel[i]);
793  }
794  if (verbosity > 1)
795  {
796  std::cout << "Particles per cell:" << std::endl;
797  for (double i : particlesPerCell)
798  {
799  std::cout << std::setw(10) << i << " ";
800  }
801  std::cout << std::endl;
802  }
803 
804  std::vector<std::vector<double> > contactWork;
805  std::vector<std::vector<double> > overheadWork;
806  contactWork.resize(NLevels);
807  overheadWork.resize(NLevels);
808  for (unsigned int i = 0; i < NLevels; i++)
809  {
810  contactWork[i].resize(NLevels);
811  overheadWork[i].resize(NLevels);
812  for (unsigned int j = 0; j < NLevels; j++)
813  {
814  contactWork[i][j] = 0;
815  overheadWork[i][j] = 0;
816  }
817  }
818 
819  for (unsigned int i = 0; i < NLevels; i++)
820  {
821  contactWork[i][i] += 0.5 * pow(3, dimension_) * particlesPerCell[i] * particlesPerLevel[i];
822  overheadWork[i][i] += 0.5 * (pow(3, dimension_) + 1) * particlesPerLevel[i];
823 
824  unsigned int startJ, endJ;
825  switch (method)
826  {
827  case BOTTOMUP:
828  {
829  startJ = i + 1;
830  endJ = NLevels;
831  break;
832  }
833  case TOPDOWN:
834  {
835  startJ = 0;
836  endJ = i;
837  break;
838  }
839  }
840  for (unsigned int j = startJ; j < endJ; j++)
841  {
842  double numberOfCellToCheck;
843  numberOfCellToCheck = expectedCellsIntegral(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1],
844  dimension_, hGridCellSizes[j + 1]);
845  contactWork[i][j] += particlesPerLevel[i] * numberOfCellToCheck * particlesPerCell[j];
846  overheadWork[i][j] += particlesPerLevel[i] * numberOfCellToCheck;
847 
848  }
849  }
850 
851  double totalContactWork = 0, totalOverheadWork = 0;
852  for (unsigned int i = 0; i < NLevels; i++)
853  {
854  for (unsigned int j = 0; j < NLevels; j++)
855  {
856  totalContactWork += contactWork[i][j];
857  totalOverheadWork += overheadWork[i][j];
858  }
859  }
860 
861  if (verbosity > 0)
862  {
863  std::cout << "Contact work: " << totalContactWork << std::endl;
864  for (unsigned int i = 0; i < NLevels; i++)
865  {
866  for (unsigned int j = 0; j < NLevels; j++)
867  {
868  std::cout << std::setw(10) << contactWork[i][j] << " ";
869  }
870  std::cout << std::endl;
871  }
872  std::cout << "Overhead work: " << totalOverheadWork << std::endl;
873  for (unsigned int i = 0; i < NLevels; i++)
874  {
875  for (unsigned int j = 0; j < NLevels; j++)
876  {
877  std::cout << std::setw(10) << overheadWork[i][j] << " ";
878  }
879  std::cout << std::endl;
880  }
881  std::cout << "Total work: " << totalContactWork + totalOverheadWork << std::endl;
882  }
883  return totalContactWork + cellCheckOverContactCheckRatio_ * totalOverheadWork;
884 }
885 
886 void HGridOptimiser::calcDfDx(std::vector<double>& hGridCellSizes, std::vector<double>& dfdx, HGridMethod method,
887  int verbosity)
888 {
889  if (verbosity > 0)
890  {
891  std::cout << "calcDfDx" << std::endl;
892  }
893  double W = calculateWork(hGridCellSizes, method, verbosity - 1);
894  double dx = 1e-7;
895  double nW;
896  dfdx.clear();
897  for (unsigned int i = 0; i < hGridCellSizes.size(); i++)
898  {
899  hGridCellSizes[i] += dx;
900  nW = calculateWork(hGridCellSizes, method, verbosity - 1);
901  dfdx.push_back((nW - W) / dx);
902  hGridCellSizes[i] -= dx;
903  if (verbosity > 0)
904  {
905  std::cout << "i=" << i << " Value=" << hGridCellSizes[i] << " Change=" << nW - W << " dfdx=" << dfdx.back()
906  << std::endl;
907  }
908  }
909 }
910 
911 double HGridOptimiser::checkLimit(std::vector<double>& hGridCellSizes, std::vector<double>& dfdx, int verbosity)
912 {
913  if (verbosity > 0)
914  {
915  std::cout << "checkLimits part 1 remove levels" << std::endl;
916  }
917 
918  dfdx[0] = 0;
919  dfdx.back() = 0;
920  for (unsigned int i = 1; i < hGridCellSizes.size(); i++)
921  {
922  if (std::abs(hGridCellSizes[i - 1] - hGridCellSizes[i]) < 1e-6)
923  {
924  if (dfdx[i - 1] < dfdx[i])
925  {
926  if (verbosity > 0)
927  {
928  std::cout << "Remove level" << i << std::endl;
929  }
930  if (i > 0.5 * hGridCellSizes.size())
931  {
932  hGridCellSizes.erase(hGridCellSizes.begin() + i - 1);
933  dfdx.erase(dfdx.begin() + i - 1);
934  }
935  else
936  {
937  hGridCellSizes.erase(hGridCellSizes.begin() + i);
938  dfdx.erase(dfdx.begin() + i);
939  }
940  i--;
941  }
942  }
943  }
944 
945  if (verbosity > 0)
946  {
947  std::cout << "checkLimits part 2 calculate limit" << std::endl;
948  for (unsigned int i = 0; i < hGridCellSizes.size(); i++)
949  {
950  std::cout << "i=" << i << " Value=" << hGridCellSizes[i] << " dfdx=" << dfdx[i] << std::endl;
951  }
952  }
953  double maxStepSize = -std::numeric_limits<double>::max();
954  double nmax;
955  for (unsigned int i = 1; i < hGridCellSizes.size(); i++)
956  {
957  nmax = (hGridCellSizes[i - 1] - hGridCellSizes[i]) / (dfdx[i] - dfdx[i - 1]);
958  if (verbosity > 0)
959  {
960  std::cout << "Max f=" << nmax << " because otherwise " << hGridCellSizes[i - 1] << "+x*" << dfdx[i - 1]
961  << ">" << hGridCellSizes[i] << "+x*" << dfdx[i] << std::endl;
962  }
963  if (nmax < 0)
964  {
965  maxStepSize = std::max(nmax, maxStepSize);
966  }
967  }
968  if (verbosity > 0)
969  {
970  std::cout << "maxStepSize=" << maxStepSize << std::endl;
971  }
972  return maxStepSize;
973 }
974 
975 void HGridOptimiser::applyStep(std::vector<double>& hGridCellSizes, std::vector<double>& dfdx, double stepsize,
976  int verbosity)
977 {
978  if (verbosity > 0)
979  {
980  std::cout << "Apply step:" << std::endl;
981  }
982  for (unsigned int i = 0; i < hGridCellSizes.size() - 1; i++)
983  {
984  hGridCellSizes[i] += stepsize * dfdx[i];
985  if (verbosity > 0)
986  {
987  std::cout << "hGridCellSizes[i]" << "+" << stepsize << "*" << dfdx[i] << "=" << hGridCellSizes[i]
988  << std::endl;
989  }
990  }
991 }
992 
993 double
994 HGridOptimiser::goldenSectionSearch(std::vector<double>& startHGridCellSizes, std::vector<double>& searchDirection,
995  double min, double cur, double max, HGridMethod method, int verbosity)
996 {
997  std::vector<double> curHGridCellSizes = startHGridCellSizes;
998  std::vector<double> xHGridCellSizes = startHGridCellSizes;
999 
1000  applyStep(curHGridCellSizes, searchDirection, cur, verbosity - 1);
1001  double curWork = calculateWork(curHGridCellSizes, method, verbosity - 1);
1002 
1003  double x;
1004  double resphi = 2.0 - 0.5 * (1 + std::sqrt(5));
1005 
1006  //Determine new probing point
1007  if (max - cur > cur - min)
1008  {
1009  //x between cur and max
1010  x = cur + resphi * (max - cur);
1011  }
1012  else
1013  {
1014  //x between min and cur
1015  x = cur - resphi * (cur - min);
1016  }
1017 
1018  if (std::abs(max - min) < 1e-7 * (std::abs(cur) + std::abs(x)))
1019  {
1020  return 0.5 * (min + max);
1021  }
1022 
1023  applyStep(xHGridCellSizes, searchDirection, x, 0);
1024  double xWork = calculateWork(xHGridCellSizes, method, 0);
1025  if (verbosity > 0)
1026  {
1027  std::cout << "min=" << min << " max=" << max << " cur=" << cur << " curWork=" << curWork << " x=" << x
1028  << " xWork=" << xWork << std::endl;
1029  }
1030  if (xWork < curWork) //x is better
1031  {
1032  if (max - cur > cur - min)
1033  {
1034  return goldenSectionSearch(startHGridCellSizes, searchDirection, cur, x, max, method, verbosity);
1035  }
1036  else
1037  {
1038  return goldenSectionSearch(startHGridCellSizes, searchDirection, min, x, cur, method, verbosity);
1039  }
1040  }
1041  else //cur is better
1042  {
1043  if (max - cur > cur - min)
1044  {
1045  return goldenSectionSearch(startHGridCellSizes, searchDirection, min, cur, x, method, verbosity);
1046  }
1047  else
1048  {
1049  return goldenSectionSearch(startHGridCellSizes, searchDirection, x, cur, max, method, verbosity);
1050  }
1051  }
1052 }
1053 
1054 void HGridOptimiser::getOptimalDistribution(std::vector<double>& hGridCellSizes, unsigned int numberOfLevels,
1055  HGridMethod method, int verbosity)
1056 {
1057  hGridCellSizes.clear();
1058  for (unsigned int i = 0; i <= numberOfLevels; i++)
1059  {
1060  hGridCellSizes.push_back(2.0 * (rMin_ + (rMax_ - rMin_) * i / numberOfLevels));
1061  }
1062 
1063  std::vector<double> dfdx;
1064  double step, max, W, oW;
1065  W = calculateWork(hGridCellSizes, method, verbosity - 3);
1066  int stepNumber = 0;
1067  do
1068  {
1069  oW = W;
1070  calculateDiffWork(hGridCellSizes, dfdx, method, verbosity - 3);
1071  max = checkLimit(hGridCellSizes, dfdx, verbosity - 3);
1072  step = goldenSectionSearch(hGridCellSizes, dfdx, max, 0.5 * max, 0, method, verbosity - 2);
1073  applyStep(hGridCellSizes, dfdx, step, verbosity - 3);
1074  W = calculateWork(hGridCellSizes, method, verbosity - 3);
1075  while (W > oW)
1076  {
1077  if (verbosity > 1)
1078  {
1079  std::cout << stepNumber << " Bad step, trying closer search range" << std::endl;
1080  }
1081  applyStep(hGridCellSizes, dfdx, -step, verbosity - 3);
1082  max /= 2;
1083  step = goldenSectionSearch(hGridCellSizes, dfdx, max, 0.5 * max, 0, method, verbosity - 2);
1084  applyStep(hGridCellSizes, dfdx, step, verbosity - 3);
1085  W = calculateWork(hGridCellSizes, method, verbosity - 3);
1086  }
1087  ++stepNumber;
1088  if (verbosity > 1)
1089  {
1090  std::cout << stepNumber << " Correct step: old work=" << oW << " new work=" << W << " difference=" << oW - W
1091  << " step=" << step / max << std::endl;
1092  }
1093  } while (oW - W > 1e-13 && stepNumber);
1094  calculateDiffWork(hGridCellSizes, dfdx, method, verbosity - 2);
1095  if (verbosity > 0)
1096  {
1097  std::cout << "Work=" << W << std::endl;
1098  std::cout << "Optimal cell sizes:";
1099  for (double hGridCellSize : hGridCellSizes)
1100  {
1101  std::cout << " " << hGridCellSize;
1102  }
1103  std::cout << std::endl;
1104  }
1105 }
1106 
1107 void HGridOptimiser::histNumberParticlesPerCell(std::vector<double>& hGridCellSizes)
1108 {
1109  unsigned int NLevels = hGridCellSizes.size() - 1;
1110 
1111  std::vector<double> particlesPerLevel;
1112  for (unsigned int i = 0; i < NLevels; i++)
1113  {
1114  particlesPerLevel.push_back(pdfInt(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], 0));
1115  if (particlesPerLevel.back() < 0)
1116  {
1117  particlesPerLevel.back() = 0;
1118  }
1119  }
1120 
1121  std::vector<double> cellsPerLevel;
1122  for (unsigned int i = 0; i < NLevels; i++)
1123  {
1124  cellsPerLevel.push_back(pow(length_ / hGridCellSizes[i + 1], dimension_));
1125  }
1126 
1127  for (unsigned int i = 0; i < NLevels; i++)
1128  {
1129  double l = particlesPerLevel[i] / cellsPerLevel[i];
1130  std::cout << "Histogram for level " << i << ":";
1131  for (unsigned int k = 0; k <= 10; k++)
1132  {
1133  std::cout << " " << std::setw(8)
1134  << std::floor(std::pow(l, k) * std::exp(-l) / mathsFunc::factorial(k) * 1e4 * cellsPerLevel[i]);
1135  }
1136  std::cout << std::endl;
1137  }
1138 }
double goldenSectionSearch(std::vector< double > &startHGridCellSizes, std::vector< double > &searchDirection, double min, double cur, double max, HGridMethod method, int verbosity)
BaseParticle * getLargestParticle() const
Returns the pointer of the largest particle in the particle handler. When mercury is running in paral...
void getOptimalDistribution(std::vector< double > &hGridCellSizes, unsigned int numberOfLevels, HGridMethod method, int verbosity)
double diffPdfInt(double x, int power)
diff(int(f(r)*r^power*dr,r=s..e)/int(f(r)*dr,r=0..omega),e)=f(e)*e^power/int(f(r)*dr,r=0..omega)
double cell2Min(unsigned int i)
Computes the left bound of the cell with given ordinal number.
unsigned int getSize() const
Gets the size of the particleHandler (including mpi and periodic particles)
Definition: BaseHandler.h:655
unsigned int getSystemDimensions() const
Returns the system dimensionality.
Definition: DPMBase.cc:1390
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Mdouble getMaxInteractionRadius() const
Returns the particle's interaction radius, which might be different from radius_ (e.g., when dealing with wet particles)
Definition: BaseParticle.h:359
void initialisePolyFunc(double omega, std::vector< double > &coeff, unsigned int numberOfCells, int verbosity)
Mdouble getZMax() const
If the length of the problem domain in z-direction is ZMax - ZMin, then getZMax() returns ZMax...
Definition: DPMBase.h:617
Mdouble getXMin() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMin() returns XMin...
Definition: DPMBase.h:586
double calculateWork(std::vector< double > &hGridCellSizes, HGridMethod method, int verbosity)
The amount of work that has to be done to run a simulation using the HGrid, in steps.
double length_
The weighted length of the domain.
double diffStartExpectedCellsIntegral(double start, double end, int p, double h)
Mdouble exp(Mdouble Exponent)
Definition: ExtendedMath.cc:84
double intCell2Min(unsigned int i)
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
Mdouble getYMin() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMin() returns YMin...
Definition: DPMBase.h:599
const std::vector< T * >::const_iterator end() const
Gets the end of the const_iterator over all BaseBoundary in this BaseHandler.
Definition: BaseHandler.h:704
double cell2Max(unsigned int i)
Computes the right bound of the cell with given ordinal number.
double intCell2Max(unsigned int i)
const std::vector< T * >::const_iterator begin() const
Gets the begin of the const_iterator over all Object in this BaseHandler.
Definition: BaseHandler.h:690
Mdouble getXMax() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMax() returns XMax...
Definition: DPMBase.h:593
void initialise(const MercuryBase &problem, unsigned int numberOfCells, int verbosity)
HGridMethod
Enum class that indicates how particles in different levels (cross level checking) of the HGrid are c...
Definition: MercuryBase.h:43
double expectedCellsIntegral(double start, double end, int p, double h)
This function calculates: int((2*r/h+2)^d f(r) dr,r=s..e)/int(f(r) dr,r=s..e)+ Used to calculated the...
void calculateDiffWork(std::vector< double > &hGridCellSizes, std::vector< double > &dfdx, HGridMethod method, int verbosity)
This is the base class for both Mercury2D and Mercury3D. Note the actually abstract grid is defined i...
Definition: MercuryBase.h:125
void histNumberParticlesPerCell(std::vector< double > &hGridCellSizes)
double diffHExpectedCellsIntegralCellNumerator(double start, double end, unsigned int i, int p, double h)
constexpr T factorial(const T t)
factorial function
Definition: ExtendedMath.h:153
void calcDfDx(std::vector< double > &hGridCellSizes, std::vector< double > &dfdx, HGridMethod method, int verbosity)
BaseParticle * getSmallestParticle() const
Gets a pointer to the smallest BaseParticle (by interactionRadius) in this ParticleHandler of the loc...
double rMax_
Radius of the largest particle, "rounded" to the smallest double that is larger than the radius of ea...
double expectedCellsIntegralCellDenominator(double start, double end, unsigned int i)
ParticleHandler particleHandler
An object of the class ParticleHandler, contains the pointers to all the particles created...
Definition: DPMBase.h:1329
double checkLimit(std::vector< double > &hGridCellSizes, std::vector< double > &dfdx, int verbosity)
std::vector< double > cellN_
void applyStep(std::vector< double > &hGridCellSizes, std::vector< double > &dfdx, double stepsize, int verbosity)
Mdouble getYMax() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMax() returns XMax...
Definition: DPMBase.h:605
double rMin_
Radius of the smallest particle, "rounded" to the largest double that is smaller than the radius of e...
std::vector< double > intCellN
unsigned int dimension_
The dimension of the system, usually 3, sometimes 2 or 1.
Mdouble getZMin() const
If the length of the problem domain in z-direction is ZMax - ZMin, then getZMin() returns ZMin...
Definition: DPMBase.h:611
double cellCheckOverContactCheckRatio_
The ratio of the time required for a single geometric contact detection over the time required to ret...
double diffHExpectedCellsIntegral(double start, double end, int p, double h)
unsigned int radius2Cell(double r)
Assigns a BaseParticle of given radius to a certain cell.
double pdfIntCell(double start, double end, unsigned int i, int p)
double pdfInt(double start, double end, int power)
This function calculates: int(f(r)*r^power*dr,r=start..end)/int(f(r)*dr,r=0..omega) with r=a*r+b...
double expectedCellsIntegralCellNumerator(double start, double end, unsigned int i, int p, double h)
unsigned int numCells_
Number of cells, usually called levels in the HGrid.
double diffEndExpectedCellsIntegral(double start, double end, int p, double h)
unsigned int radius2IntCell(double r)