MercuryDPM  Beta
 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-2014, 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 <math.h>
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(std::vector<double>::iterator it = cellN_.begin(); it != cellN_.end(); ++it)
41  {
42  *it = 0.0;
43  }
44 
45  //Set the minimum and maximum radius of the particles.
46  rMin_ = nextafter(problem.particleHandler.getSmallestParticle()->getInteractionRadius(), 0.0);
47  rMax_ = nextafter(nextafter(problem.particleHandler.getLargestParticle()->getInteractionRadius(), std::numeric_limits<double>::max()), std::numeric_limits<double>::max());
48 
49  //for all particles, add one to the cell in cellN_ which has is associated with
50  //the radius of that particle.
51  for(std::vector<BaseParticle*>::const_iterator it = problem.particleHandler.begin(); it != problem.particleHandler.end(); ++it)
52  {
53  cellN_[radius2Cell((*it)->getInteractionRadius())]++;
54  }
55 
56  //assign what the number of particles is.
57  unsigned int numParticles = problem.particleHandler.getNumberOfObjects();
58 
59  intCellN.push_back(1.5 * cellN_[0] - 0.5 * cellN_[1]);
60  for(std::vector<double>::iterator it = cellN_.begin(); it != cellN_.end(); ++it)
61  {
62  intCellN.push_back(*it);
63  }
64  intCellN.push_back(1.5 * cellN_[numCells_ - 1] - 0.5 * cellN_[numCells_ - 2]);
65 
66  if(verbosity > 0)
67  {
68  std::cout << "rMin=" << rMin_ << " rMax=" << rMax_ << " NParticles=" << numParticles << std::endl;
69  for(unsigned int i = 0; i < cellN_.size(); i++)
70  {
71  std::cout << "From" << std::setw(8) << cell2Min(i) << " to" << std::setw(8) << cell2Max(i) << " number=" << std::setw(5) << cellN_[i] << std::endl;
72  }
73  std::cout << "Domain size [" << problem.getXMin() << "," << problem.getXMax() << "]x[" << problem.getYMin() << "," << problem.getYMax() << "]x[" << problem.getZMin() << "," << problem.getZMax() << "]" << std::endl;
74 
75  for(unsigned int i = 0; i < intCellN.size() - 1; i++)
76  {
77  std::cout << "From" << std::setw(8) << intCell2Min(i) << " to" << std::setw(8) << intCell2Max(i) << " linear from " << std::setw(8) << intCellN[i] << " to " << std::setw(8) << intCellN[i + 1] << std::endl;
78  }
79  /*std::cout<<"["<<cellN_[0];
80  for (unsigned int i=1;i<cellN_.size();i++)
81  {
82  std::cout<<","<<cellN_[i];
83  }
84  std::cout<<"]"<<std::endl;*/
85  }
86 
87  dimension_ = problem.getSystemDimensions();
88  switch(dimension_)
89  {
90  case 1:
91  length_ = pow((problem.getXMax() - problem.getXMin()) / numParticles, 1);
92  break;
93  case 2:
94  length_ = pow((problem.getXMax() - problem.getXMin()) * (problem.getYMax() - problem.getYMin()) / numParticles, 1.0 / 2.0);
95  break;
96  case 3:
97  length_ = pow((problem.getXMax() - problem.getXMin()) * (problem.getYMax() - problem.getYMin()) * (problem.getZMax() - problem.getZMin()) / numParticles, 1.0 / 3.0);
98  break;
99  }
100 }
101 
102 void HGridOptimiser::initialisePolyFunc(double omega, std::vector<double>& coeff, unsigned int numberOfCells, int verbosity)
103 {
104  numCells_ = numberOfCells;
106  cellN_.resize(numCells_);
107  for(std::vector<double>::iterator it = cellN_.begin(); it != cellN_.end(); ++it)
108  {
109  *it = 0;
110  }
111 
112  rMin_ = nextafter(1, 0.0);
113  rMax_ = nextafter(omega, std::numeric_limits<double>::max());
114  for(unsigned int i = 0; i < cellN_.size(); i++)
115  {
116  double start = cell2Min(i);
117  double end = cell2Max(i);
118  for(unsigned int j = 0; j < coeff.size(); j++)
119  {
120  cellN_[i] += coeff[j] / (1.0 + j) * (std::pow(end, j + 1) - std::pow(start, j + 1));
121  }
122  }
123 
124  intCellN.push_back(1.5 * cellN_[0] - 0.5 * cellN_[1]);
125  for(std::vector<double>::iterator it = cellN_.begin(); it != cellN_.end(); ++it)
126  {
127  intCellN.push_back(*it);
128  }
129  intCellN.push_back(1.5 * cellN_[numCells_ - 1] - 0.5 * cellN_[numCells_ - 2]);
130 
131  dimension_ = 1;
132  length_ = 1.0;
133 
134  if(verbosity > 0)
135  {
136  std::cout << "rMin=" << rMin_ << " rMax=" << rMax_ << std::endl;
137  for(unsigned int i = 0; i < cellN_.size(); i++)
138  {
139  std::cout << "From" << std::setw(8) << cell2Min(i) << " to" << std::setw(8) << cell2Max(i) << " number=" << std::setw(5) << cellN_[i] << std::endl;
140  }
141 
142  for(unsigned int i = 0; i < intCellN.size() - 1; i++)
143  {
144  std::cout << "From" << std::setw(8) << intCell2Min(i) << " to" << std::setw(8) << intCell2Max(i) << " linear from " << std::setw(8) << intCellN[i] << " to " << std::setw(8) << intCellN[i + 1] << std::endl;
145  }
146  }
147 }
148 
156 unsigned int HGridOptimiser::radius2Cell(double r)
157 {
158  if(r < rMin_ || r >= rMax_)
159  logger(ERROR, "The radius % is not in the range [%,%]", r, rMin_, rMax_);
160  unsigned int y = static_cast<unsigned int> (floor(numCells_ * (r - rMin_) / (rMax_ - rMin_)));
161  return y;
162 }
163 
164 unsigned int HGridOptimiser::radius2IntCell(double r)
165 {
166  if(r < rMin_ || r >= rMax_)
167  throw;
168  unsigned int y = static_cast<unsigned int> (floor(numCells_ * (r - rMin_) / (rMax_ - rMin_) + 0.5));
169  return y;
170 }
171 
172 double HGridOptimiser::intCell2Min(unsigned int i)
173 {
174  if(i == 0)
175  return rMin_;
176  else
177  return rMin_ + (rMax_ - rMin_) * (i - 0.5) / numCells_;
178 }
179 
180 double HGridOptimiser::intCell2Max(unsigned int i)
181 {
182  if(i == numCells_)
183  return rMax_;
184  else
185  return rMin_ + (rMax_ - rMin_) * (i + 0.5) / numCells_;
186 }
187 
191 double HGridOptimiser::cell2Min(unsigned int i)
192 {
193  return rMin_ + (rMax_ - rMin_) * i / numCells_;
194 }
195 
199 double HGridOptimiser::cell2Max(unsigned int i)
200 {
201  return rMin_ + (rMax_ - rMin_) * (i + 1) / numCells_;
202 }
203 
204 double HGridOptimiser::pdfIntCell(double start, double end, unsigned int i, int p)
205 {
206  double a = (intCellN[i + 1] - intCellN[i]) / (intCell2Max(i) - intCell2Min(i));
207  double b = (intCellN[i] * intCell2Max(i) - intCellN[i + 1] * intCell2Min(i)) / (intCell2Max(i) - intCell2Min(i));
208  return(a / (p + 2) * (std::pow(end, p + 2) - std::pow(start, p + 2)) + b / (p + 1) * (std::pow(end, p + 1) - std::pow(start, p + 1)));
209 }
210 
214 double HGridOptimiser::pdfInt(double start, double end, int p)
215 {
216  //std::cout<<"pdfInit("<<start<<","<<end<<","<<p<<");"<<std::endl;
217  unsigned int startCell = radius2IntCell(start);
218  unsigned int endCell = radius2IntCell(end);
219 
220  double numerator = 0;
221  if(startCell == endCell)
222  {
223  numerator = pdfIntCell(start, end, startCell, p);
224  }
225  else if(endCell < startCell)
226  {
227  numerator = 0;
228  }
229  else
230  {
231  numerator = pdfIntCell(start, intCell2Max(startCell), startCell, p);
232  //std::cout<<"Teller+="<<pdfIntCell(start,intCell2Max(startCell),startCell,p)<<std::endl;
233  for(unsigned int i = startCell + 1; i < endCell; i++)
234  {
235  numerator += pdfIntCell(intCell2Min(i), intCell2Max(i), i, p);
236  //std::cout<<"Teller+="<<pdfIntCell(intCell2Min(i),intCell2Max(i),i,p)<<std::endl;
237  }
238  numerator += pdfIntCell(intCell2Min(endCell), end, endCell, p);
239  //std::cout<<"Teller+="<<pdfIntCell(intCell2Min(endCell),end,endCell,p)<<std::endl;
240  }
241  double denominator = 0;
242  for(unsigned int i = 0; i <= numCells_; i++)
243  {
244  denominator += pdfIntCell(intCell2Min(i), intCell2Max(i), i, 0);
245  //std::cout<<"Noemer+="<<pdfIntCell(intCell2Min(i),intCell2Max(i),i,0)<<std::endl;
246  }
247  return numerator / denominator;
248 }
249 
251 double HGridOptimiser::diffPdfInt(double x, int p)
252 {
253  unsigned int xCell = radius2IntCell(x);
254  double denominator = 0;
255  for(unsigned int i = 0; i <= numCells_; i++)
256  {
257  denominator += pdfIntCell(intCell2Min(i), intCell2Max(i), i, 0);
258  }
259  double a = (intCellN[xCell + 1] - intCellN[xCell]) / (intCell2Max(xCell) - intCell2Min(xCell));
260  double b = (intCellN[xCell] * intCell2Max(xCell) - intCellN[xCell + 1] * intCell2Min(xCell)) / (intCell2Max(xCell) - intCell2Min(xCell));
261  return(a * std::pow(x, p + 1) + b * std::pow(x, p)) / denominator;
262 }
263 
264 double HGridOptimiser::expectedCellsIntegralCellNumerator(double start, double end, unsigned int i, int p, double h)
265 {
266  double a = (intCellN[i + 1] - intCellN[i]) / (intCell2Max(i) - intCell2Min(i));
267  double b = (intCellN[i] * intCell2Max(i) - intCellN[i + 1] * intCell2Min(i)) / (intCell2Max(i) - intCell2Min(i));
268  double r;
269  r = start;
270  double min = std::pow(2.0 * (r + h) / h, p) * (a * p * r - a * h + a * r + b * p + 2 * b) * (r + h) / ((p + 1) * (p + 2));
271  r = end;
272  double plus = std::pow(2.0 * (r + h) / h, p) * (a * p * r - a * h + a * r + b * p + 2 * b) * (r + h) / ((p + 1) * (p + 2));
273  return plus - min;
274 }
275 
276 double HGridOptimiser::diffHExpectedCellsIntegralCellNumerator(double start, double end, unsigned int i, int p, double h)
277 {
278  double a = (intCellN[i + 1] - intCellN[i]) / (intCell2Max(i) - intCell2Min(i));
279  double b = (intCellN[i] * intCell2Max(i) - intCellN[i + 1] * intCell2Min(i)) / (intCell2Max(i) - intCell2Min(i));
280  double r;
281  r = start;
282  //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));
283  double min = (-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) * (r + h) +
284  -std::pow(2.0 * r / h + 2.0, p) * a * (r + h) +
285  std::pow(2.0 * r / h + 2.0, p) * (a * p * r - a * h + a * r + b * p + 2 * b)) / ((p + 1) * (p + 2));
286  r = end;
287  //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));
288  double plus = (-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) * (r + h) +
289  -std::pow(2.0 * r / h + 2.0, p) * a * (r + h) +
290  std::pow(2.0 * r / h + 2.0, p) * (a * p * r - a * h + a * r + b * p + 2 * b)) / ((p + 1) * (p + 2));
291  return plus - min;
292 }
293 
294 double HGridOptimiser::expectedCellsIntegralCellDenominator(double start, double end, unsigned int i)
295 {
296  double a = (intCellN[i + 1] - intCellN[i]) / (intCell2Max(i) - intCell2Min(i));
297  double b = (intCellN[i] * intCell2Max(i) - intCellN[i + 1] * intCell2Min(i)) / (intCell2Max(i) - intCell2Min(i));
298  return(a / 2.0 * (pow(end, 2) - pow(start, 2)) + b * (pow(end, 1) - pow(start, 1)));
299 }
300 
305 double HGridOptimiser::expectedCellsIntegral(double start, double end, int p, double h)
306 {
307  unsigned int startCell = radius2IntCell(start);
308  unsigned int endCell = radius2IntCell(end);
309 
310  double numerator = 0;
311  double denominator = 0;
312  if(startCell == endCell)
313  {
314  numerator = expectedCellsIntegralCellNumerator(start, end, startCell, p, h);
315  denominator = expectedCellsIntegralCellDenominator(start, end, startCell);
316  }
317  else
318  {
319  numerator = expectedCellsIntegralCellNumerator(start, intCell2Max(startCell), startCell, p, h);
320  denominator = expectedCellsIntegralCellDenominator(start, intCell2Max(startCell), startCell);
321  for(unsigned int i = startCell + 1; i < endCell; i++)
322  {
323  numerator += expectedCellsIntegralCellNumerator(intCell2Min(i), intCell2Max(i), i, p, h);
325  }
326  numerator += expectedCellsIntegralCellNumerator(intCell2Min(endCell), end, endCell, p, h);
327  denominator += expectedCellsIntegralCellDenominator(intCell2Min(endCell), end, endCell);
328 
329  }
330  // std::cout<<"New: numerator="<<numerator<<" denominator="<<denominator<<std::endl;
331  return numerator / denominator;
332 }
333 
334 double
335 HGridOptimiser::diffStartExpectedCellsIntegral(double start, double end, int p, double h)
336 {
337  unsigned int startCell = radius2IntCell(start);
338  unsigned int endCell = radius2IntCell(end);
339 
340  double numerator = 0;
341  double denominator = 0;
342 
343  double a = (intCellN[startCell + 1] - intCellN[startCell]) / (intCell2Max(startCell) - intCell2Min(startCell));
344  double b = (intCellN[startCell] * intCell2Max(startCell) - intCellN[startCell + 1] * intCell2Min(startCell)) / (intCell2Max(startCell) - intCell2Min(startCell));
345 
346  double diffTeller = -(2.0 / h * p * std::pow(2.0 * (start + h) / h, p - 1) * (a * p * start - a * h + a * start + b * p + 2 * b) * (start + h) / ((p + 1) * (p + 2)) +
347  std::pow(2.0 * (start + h) / h, p) * (a * p + a) * (start + h) / ((p + 1) * (p + 2)) +
348  std::pow(2.0 * (start + h) / h, p) * (a * p * start - a * h + a * start + b * p + 2 * b) / ((p + 1) * (p + 2)));
349  double diffNoemer = -(a * start + b);
350  if(startCell == endCell)
351  {
352  numerator = expectedCellsIntegralCellNumerator(start, end, startCell, p, h);
353  denominator = expectedCellsIntegralCellDenominator(start, end, startCell);
354  }
355  else
356  {
357  numerator = expectedCellsIntegralCellNumerator(start, intCell2Max(startCell), startCell, p, h);
358  denominator = expectedCellsIntegralCellDenominator(start, intCell2Max(startCell), startCell);
359  for(unsigned int i = startCell + 1; i < endCell; i++)
360  {
361  numerator += expectedCellsIntegralCellNumerator(intCell2Min(i), intCell2Max(i), i, p, h);
363  }
364  numerator += expectedCellsIntegralCellNumerator(intCell2Min(endCell), end, endCell, p, h);
365  denominator += expectedCellsIntegralCellDenominator(intCell2Min(endCell), end, endCell);
366 
367  }
368  //std::cout<<"new: numerator="<<numerator<<" denominator="<<denominator<<" diffTeller="<<diffTeller<<" diffNoemer="<<diffNoemer<<std::endl;
369  return(diffTeller * denominator - numerator * diffNoemer) / (denominator * denominator);
370 }
371 
372 double HGridOptimiser::diffEndExpectedCellsIntegral(double start, double end, int p, double h)
373 {
374  unsigned int startCell = radius2IntCell(start);
375  unsigned int endCell = radius2IntCell(end);
376 
377  double numerator = 0;
378  double denominator = 0;
379 
380  double a = (intCellN[endCell + 1] - intCellN[endCell]) / (intCell2Max(endCell) - intCell2Min(endCell));
381  double b = (intCellN[endCell] * intCell2Max(endCell) - intCellN[endCell + 1] * intCell2Min(endCell)) / (intCell2Max(endCell) - intCell2Min(endCell));
382 
383  double diffTeller = (2.0 / h * p * std::pow(2.0 * (end + h) / h, p - 1) * (a * p * end - a * h + a * end + b * p + 2 * b) * (end + h) / ((p + 1) * (p + 2)) +
384  std::pow(2.0 * (end + h) / h, p) * (a * p + a) * (end + h) / ((p + 1) * (p + 2)) +
385  std::pow(2.0 * (end + h) / h, p) * (a * p * end - a * h + a * end + b * p + 2 * b) / ((p + 1) * (p + 2)));
386  double diffNoemer = (a * end + b);
387  if(startCell == endCell)
388  {
389  numerator = expectedCellsIntegralCellNumerator(start, end, startCell, p, h);
390  denominator = expectedCellsIntegralCellDenominator(start, end, startCell);
391  }
392  else
393  {
394  numerator = expectedCellsIntegralCellNumerator(start, intCell2Max(startCell), startCell, p, h);
395  denominator = expectedCellsIntegralCellDenominator(start, intCell2Max(startCell), startCell);
396  for(unsigned int i = startCell + 1; i < endCell; i++)
397  {
398  numerator += expectedCellsIntegralCellNumerator(intCell2Min(i), intCell2Max(i), i, p, h);
400  }
401  numerator += expectedCellsIntegralCellNumerator(intCell2Min(endCell), end, endCell, p, h);
402  denominator += expectedCellsIntegralCellDenominator(intCell2Min(endCell), end, endCell);
403 
404  }
405  //std::cout<<"new: numerator="<<numerator<<" denominator="<<denominator<<" diffTeller="<<diffTeller<<" diffNoemer="<<diffNoemer<<std::endl;
406  return(diffTeller * denominator - numerator * diffNoemer) / (denominator * denominator);
407 }
408 
409 double HGridOptimiser::diffHExpectedCellsIntegral(double start, double end, int p, double h)
410 {
411  unsigned int startCell = radius2IntCell(start);
412  unsigned int endCell = radius2IntCell(end);
413 
414  double denominator = 0;
415  double diffTeller = 0;
416 
417  if(startCell == endCell)
418  {
419  diffTeller = diffHExpectedCellsIntegralCellNumerator(start, end, startCell, p, h);
420  denominator = expectedCellsIntegralCellDenominator(start, end, startCell);
421  }
422  else
423  {
424  diffTeller = diffHExpectedCellsIntegralCellNumerator(start, intCell2Max(startCell), startCell, p, h);
425  denominator = expectedCellsIntegralCellDenominator(start, intCell2Max(startCell), startCell);
426  for(unsigned int i = startCell + 1; i < endCell; i++)
427  {
428  diffTeller += diffHExpectedCellsIntegralCellNumerator(intCell2Min(i), intCell2Max(i), i, p, h);
430  }
431  diffTeller += diffHExpectedCellsIntegralCellNumerator(intCell2Min(endCell), end, endCell, p, h);
432  denominator += expectedCellsIntegralCellDenominator(intCell2Min(endCell), end, endCell);
433  }
434  return diffTeller / denominator;
435 }
436 
437 void HGridOptimiser::calculateDiffWork(std::vector<double>& hGridCellSizes, std::vector<double>& dfdx, HGridMethod method, int verbosity)
438 {
439  unsigned int NLevels = hGridCellSizes.size() - 1;
440  unsigned int NParams = hGridCellSizes.size();
441 
442  if(verbosity > 0)
443  {
444  std::cout << "Length scale=" << length_ << std::endl;
445  }
446  if(verbosity > 0)
447  {
448  std::cout << "Level sizes:";
449  for(unsigned int i = 0; i < NParams; i++)
450  {
451  std::cout << " " << hGridCellSizes[i];
452  }
453  std::cout << std::endl;
454  }
455 
456  std::vector<double> particlesPerLevel;
457  for(unsigned int i = 0; i < NLevels; i++)
458  {
459  particlesPerLevel.push_back(pdfInt(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], 0));
460  if(particlesPerLevel.back() < 0)
461  {
462  particlesPerLevel.back() = 0;
463  }
464  }
465  if(verbosity > 0)
466  {
467  std::cout << "Particles per level: ";
468  for(unsigned int i = 0; i < NLevels; i++)
469  {
470  std::cout << " " << particlesPerLevel[i];
471  }
472  std::cout << std::endl;
473  }
474 
475  //How changes particlesPerLeve[i] when hGridCellSize[j] is changed
476  std::vector<std::vector<double> > diffParticlesPerLevel;
477  diffParticlesPerLevel.resize(NLevels);
478  for(unsigned int i = 0; i < NLevels; i++)
479  {
480  for(unsigned int j = 0; j < NParams; j++)
481  {
482  diffParticlesPerLevel[i].push_back(0.0);
483  if(j == i)
484  {
485  diffParticlesPerLevel[i][j] += -0.5 * diffPdfInt(0.5 * hGridCellSizes[i], 0);
486  }
487  if(j == i + 1)
488  {
489  diffParticlesPerLevel[i][j] += 0.5 * diffPdfInt(0.5 * hGridCellSizes[i + 1], 0);
490  }
491  }
492  }
493  if(verbosity > 0)
494  {
495  std::cout << "Diff Particles per level: " << std::endl;
496  for(unsigned int i = 0; i < NLevels; i++)
497  {
498  for(unsigned int j = 0; j < NParams; j++)
499  {
500  std::cout << " " << std::setw(12) << diffParticlesPerLevel[i][j];
501  }
502  std::cout << std::endl;
503  }
504  }
505 
506  std::vector<double> cellsPerLevel;
507  for(unsigned int i = 0; i < NLevels; i++)
508  {
509  cellsPerLevel.push_back(pow(length_ / hGridCellSizes[i + 1], dimension_));
510  }
511  if(verbosity > 0)
512  {
513  std::cout << "Cells per level: ";
514  for(unsigned int i = 0; i < NLevels; i++)
515  {
516  std::cout << " " << cellsPerLevel[i];
517  }
518  std::cout << std::endl;
519  }
520 
521  //How changes cellsPerLevel[i] when hGridCellSize[j] is changed
522  std::vector<std::vector<double> > diffCellsPerLevel;
523  diffCellsPerLevel.resize(hGridCellSizes.size());
524  for(unsigned int i = 0; i < NLevels; i++)
525  {
526  for(unsigned int j = 0; j < NParams; j++)
527  {
528  if(j == i + 1)
529  {
530  diffCellsPerLevel[i].push_back(-pow(length_ / hGridCellSizes[i + 1], dimension_) * dimension_ / hGridCellSizes[i + 1]);
531  }
532  else
533  {
534  diffCellsPerLevel[i].push_back(0.0);
535  }
536  }
537  }
538  if(verbosity > 0)
539  {
540  std::cout << "Diff Cells per level: " << std::endl;
541  for(unsigned int i = 0; i < NLevels; i++)
542  {
543  for(unsigned int j = 0; j < NParams; j++)
544  {
545  std::cout << " " << std::setw(12) << diffCellsPerLevel[i][j];
546  }
547  std::cout << std::endl;
548  }
549  }
550 
551  std::vector<double> particlesPerCell;
552  for(unsigned int i = 0; i < NLevels; i++)
553  {
554  particlesPerCell.push_back(particlesPerLevel[i] / cellsPerLevel[i]);
555  }
556  if(verbosity > 0)
557  {
558  std::cout << "Particles per cell: ";
559  for(unsigned int i = 0; i < particlesPerCell.size(); i++)
560  {
561  std::cout << " " << particlesPerCell[i];
562  }
563  std::cout << std::endl;
564  }
565 
566  //How changes particlesPerCell[i] when hGridCellSize[j] is changed
567  std::vector<std::vector<double> > diffParticlesPerCell;
568  diffParticlesPerCell.resize(NLevels);
569  for(unsigned int i = 0; i < NLevels; i++)
570  {
571  for(unsigned int j = 0; j < NParams; j++)
572  {
573  diffParticlesPerCell[i].push_back((diffParticlesPerLevel[i][j] * cellsPerLevel[i] - particlesPerLevel[i] * diffCellsPerLevel[i][j]) / (cellsPerLevel[i] * cellsPerLevel[i]));
574  }
575  }
576  if(verbosity > 0)
577  {
578  std::cout << "Diff Particles per Cell: " << std::endl;
579  for(unsigned int i = 0; i < NLevels; i++)
580  {
581  for(unsigned int j = 0; j < NParams; j++)
582  {
583  std::cout << " " << std::setw(12) << diffParticlesPerCell[i][j];
584  }
585  std::cout << std::endl;
586  }
587  }
588 
589  double contactWork = 0, overheadWork = 0;
590  std::vector<double> diffContactWork;
591  std::vector<double> diffOverheadWork;
592  diffContactWork.resize(NParams);
593  diffOverheadWork.resize(NParams);
594  for(unsigned int j = 0; j < NParams; j++)
595  {
596  diffContactWork[j] = 0;
597  diffOverheadWork[j] = 0;
598  }
599 
600  for(unsigned int i = 0; i < NLevels; i++)
601  {
602  contactWork += 0.5 * pow(3, dimension_) * particlesPerCell[i] * particlesPerLevel[i];
603  overheadWork += 0.5 * (pow(3, dimension_) + 1) * particlesPerLevel[i];
604  for(unsigned int j = 0; j < NParams; j++)
605  {
606  diffContactWork[j] += 0.5 * pow(3, dimension_) * (diffParticlesPerCell[i][j] * particlesPerLevel[i] + particlesPerCell[i] * diffParticlesPerLevel[i][j]);
607  diffOverheadWork[j] += 0.5 * (pow(3, dimension_) + 1) * diffParticlesPerLevel[i][j];
608  }
609 
610  unsigned int startJ, endJ;
611  switch(method)
612  {
613  case BOTTOMUP:
614  {
615  startJ = i + 1;
616  endJ = NLevels;
617  break;
618  }
619  case TOPDOWN:
620  {
621  startJ = 0;
622  endJ = i;
623  break;
624  }
625  }
626  for(unsigned int j = startJ; j < endJ; j++)
627  {
628  double numberOfCellToCheck;
629  numberOfCellToCheck = expectedCellsIntegral(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], dimension_, hGridCellSizes[j + 1]);
630 
631  std::vector<double> diffNumberOfCellToCheck;
632  for(unsigned int k = 0; k < NParams; k++)
633  {
634  diffNumberOfCellToCheck.push_back(0.0);
635  if(k == i)
636  {
637  diffNumberOfCellToCheck[k] += 0.5 * diffStartExpectedCellsIntegral(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], dimension_, hGridCellSizes[j + 1]);
638  }
639  if(k == i + 1)
640  {
641  diffNumberOfCellToCheck[k] += 0.5 * diffEndExpectedCellsIntegral(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], dimension_, hGridCellSizes[j + 1]);
642  }
643  if(k == j + 1)
644  {
645  diffNumberOfCellToCheck[k] += diffHExpectedCellsIntegral(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], dimension_, hGridCellSizes[j + 1]);
646  }
647  }
648  contactWork += particlesPerLevel[i] * numberOfCellToCheck * particlesPerCell[j];
649  overheadWork += particlesPerLevel[i] * numberOfCellToCheck;
650  //std::cout<<"i= "<<i<<" j= "<<j<<" NumberOfCellToCheck= "<<numberOfCellToCheck<<" diffNumberOfCellToCheck[2]= "<<diffNumberOfCellToCheck[2]<<std::endl;
651  for(unsigned int k = 0; k < NParams; k++)
652  {
653  diffContactWork[k] += (diffParticlesPerLevel[i][k] * numberOfCellToCheck * particlesPerCell[j] + particlesPerLevel[i] * diffNumberOfCellToCheck[k] * particlesPerCell[j] + particlesPerLevel[i] * numberOfCellToCheck * diffParticlesPerCell[j][k]);
654  diffOverheadWork[k] += (diffParticlesPerLevel[i][k] * numberOfCellToCheck + particlesPerLevel[i] * diffNumberOfCellToCheck[k]);
655  }
656  }
657  }
658  if(verbosity > 0)
659  {
660  std::cout << "Contact work: " << contactWork << std::endl;
661  std::cout << "Overhead work: " << cellCheckOverContactCheckRatio_ * overheadWork << std::endl;
662  std::cout << "Sum work: " << contactWork + cellCheckOverContactCheckRatio_ * overheadWork << std::endl;
663  std::cout << "diff Contactwork: ";
664  for(unsigned int j = 0; j < NParams; j++)
665  {
666  std::cout << " " << diffContactWork[j];
667  }
668  std::cout << std::endl;
669  std::cout << "diff Overheadwork: ";
670  for(unsigned int j = 0; j < NParams; j++)
671  {
672  std::cout << " " << cellCheckOverContactCheckRatio_ * diffOverheadWork[j];
673  }
674  std::cout << std::endl;
675  std::cout << "diff sum: ";
676  for(unsigned int j = 0; j < NParams; j++)
677  {
678  std::cout << " " << diffContactWork[j] + cellCheckOverContactCheckRatio_ * diffOverheadWork[j];
679  }
680  std::cout << std::endl;
681  }
682  dfdx.clear();
683  for(unsigned int j = 0; j < NParams; j++)
684  {
685  dfdx.push_back(diffContactWork[j] + cellCheckOverContactCheckRatio_ * diffOverheadWork[j]);
686  }
687 }
688 
689 double HGridOptimiser::calculateWork(std::vector<double>& hGridCellSizes, HGridMethod method, int verbosity)
690 {
691  unsigned int NLevels = hGridCellSizes.size() - 1;
692  unsigned int NParams = hGridCellSizes.size();
693 
694  if(verbosity > 1)
695  {
696  std::cout << "Length scale=" << length_ << std::endl;
697  }
698  if(verbosity > 0)
699  {
700  std::cout << "Level sizes:" << std::endl;
701  for(unsigned int i = 0; i < NParams; i++)
702  {
703  std::cout << std::setw(10) << hGridCellSizes[i] << " ";
704  }
705  std::cout << std::endl;
706  }
707 
708  std::vector<double> particlesPerLevel;
709  for(unsigned int i = 0; i < NLevels; i++)
710  {
711  particlesPerLevel.push_back(pdfInt(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], 0));
712  if(particlesPerLevel.back() < 0)
713  {
714  particlesPerLevel.back() = 0;
715  }
716  }
717  if(verbosity > 0)
718  {
719  std::cout << "Particles per level:" << std::endl;
720  for(unsigned int i = 0; i < NLevels; i++)
721  {
722  std::cout << std::setw(10) << particlesPerLevel[i] << " ";
723  }
724  std::cout << std::endl;
725  }
726 
727  std::vector<double> cellsPerLevel;
728  for(unsigned int i = 0; i < NLevels; i++)
729  {
730  cellsPerLevel.push_back(pow(length_ / hGridCellSizes[i + 1], dimension_));
731  }
732  if(verbosity > 1)
733  {
734  std::cout << "Cells per level:" << std::endl;
735  for(unsigned int i = 0; i < NLevels; i++)
736  {
737  std::cout << std::setw(10) << cellsPerLevel[i] << " ";
738  }
739  std::cout << std::endl;
740  }
741 
742  std::vector<double> particlesPerCell;
743  for(unsigned int i = 0; i < NLevels; i++)
744  {
745  particlesPerCell.push_back(particlesPerLevel[i] / cellsPerLevel[i]);
746  }
747  if(verbosity > 1)
748  {
749  std::cout << "Particles per cell:" << std::endl;
750  for(unsigned int i = 0; i < particlesPerCell.size(); i++)
751  {
752  std::cout << std::setw(10) << particlesPerCell[i] << " ";
753  }
754  std::cout << std::endl;
755  }
756 
757  std::vector<std::vector<double> > contactWork;
758  std::vector<std::vector<double> > overheadWork;
759  contactWork.resize(NLevels);
760  overheadWork.resize(NLevels);
761  for(unsigned int i = 0; i < NLevels; i++)
762  {
763  contactWork[i].resize(NLevels);
764  overheadWork[i].resize(NLevels);
765  for(unsigned int j = 0; j < NLevels; j++)
766  {
767  contactWork[i][j] = 0;
768  overheadWork[i][j] = 0;
769  }
770  }
771 
772  for(unsigned int i = 0; i < NLevels; i++)
773  {
774  contactWork[i][i] += 0.5 * pow(3, dimension_) * particlesPerCell[i] * particlesPerLevel[i];
775  overheadWork[i][i] += 0.5 * (pow(3, dimension_) + 1) * particlesPerLevel[i];
776 
777  unsigned int startJ, endJ;
778  switch(method)
779  {
780  case BOTTOMUP:
781  {
782  startJ = i + 1;
783  endJ = NLevels;
784  break;
785  }
786  case TOPDOWN:
787  {
788  startJ = 0;
789  endJ = i;
790  break;
791  }
792  }
793  for(unsigned int j = startJ; j < endJ; j++)
794  {
795  double numberOfCellToCheck;
796  numberOfCellToCheck = expectedCellsIntegral(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], dimension_, hGridCellSizes[j + 1]);
797  contactWork[i][j] += particlesPerLevel[i] * numberOfCellToCheck * particlesPerCell[j];
798  overheadWork[i][j] += particlesPerLevel[i] * numberOfCellToCheck;
799 
800  }
801  }
802 
803  double totalContactWork = 0, totalOverheadWork = 0;
804  for(unsigned int i = 0; i < NLevels; i++)
805  {
806  for(unsigned int j = 0; j < NLevels; j++)
807  {
808  totalContactWork += contactWork[i][j];
809  totalOverheadWork += overheadWork[i][j];
810  }
811  }
812 
813  if(verbosity > 0)
814  {
815  std::cout << "Contact work: " << totalContactWork << std::endl;
816  for(unsigned int i = 0; i < NLevels; i++)
817  {
818  for(unsigned int j = 0; j < NLevels; j++)
819  {
820  std::cout << std::setw(10) << contactWork[i][j] << " ";
821  }
822  std::cout << std::endl;
823  }
824  std::cout << "Overhead work: " << totalOverheadWork << std::endl;
825  for(unsigned int i = 0; i < NLevels; i++)
826  {
827  for(unsigned int j = 0; j < NLevels; j++)
828  {
829  std::cout << std::setw(10) << overheadWork[i][j] << " ";
830  }
831  std::cout << std::endl;
832  }
833  std::cout << "Total work: " << totalContactWork + totalOverheadWork << std::endl;
834  }
835  return totalContactWork + cellCheckOverContactCheckRatio_ * totalOverheadWork;
836 }
837 
838 void HGridOptimiser::calcDfDx(std::vector<double>& hGridCellSizes, std::vector<double>& dfdx, HGridMethod method, int verbosity)
839 {
840  if(verbosity > 0)
841  {
842  std::cout << "calcDfDx" << std::endl;
843  }
844  double W = calculateWork(hGridCellSizes, method, verbosity - 1);
845  double dx = 1e-7;
846  double nW;
847  dfdx.clear();
848  for(unsigned int i = 0; i < hGridCellSizes.size(); i++)
849  {
850  hGridCellSizes[i] += dx;
851  nW = calculateWork(hGridCellSizes, method, verbosity - 1);
852  dfdx.push_back((nW - W) / dx);
853  hGridCellSizes[i] -= dx;
854  if(verbosity > 0)
855  {
856  std::cout << "i=" << i << " Value=" << hGridCellSizes[i] << " Change=" << nW - W << " dfdx=" << dfdx.back() << std::endl;
857  }
858  }
859 }
860 
861 double HGridOptimiser::checkLimit(std::vector<double>& hGridCellSizes, std::vector<double>& dfdx, int verbosity)
862 {
863  if(verbosity > 0)
864  {
865  std::cout << "checkLimits part 1 remove levels" << std::endl;
866  }
867 
868  dfdx[0] = 0;
869  dfdx.back() = 0;
870  for(unsigned int i = 1; i < hGridCellSizes.size(); i++)
871  {
872  if(std::abs(hGridCellSizes[i - 1] - hGridCellSizes[i]) < 1e-6)
873  {
874  if(dfdx[i - 1] < dfdx[i])
875  {
876  if(verbosity > 0)
877  {
878  std::cout << "Remove level" << i << std::endl;
879  }
880  if(i > 0.5 * hGridCellSizes.size())
881  {
882  hGridCellSizes.erase(hGridCellSizes.begin() + i - 1);
883  dfdx.erase(dfdx.begin() + i - 1);
884  }
885  else
886  {
887  hGridCellSizes.erase(hGridCellSizes.begin() + i);
888  dfdx.erase(dfdx.begin() + i);
889  }
890  i--;
891  }
892  }
893  }
894 
895  if(verbosity > 0)
896  {
897  std::cout << "checkLimits part 2 calculate limit" << std::endl;
898  for(unsigned int i = 0; i < hGridCellSizes.size(); i++)
899  {
900  std::cout << "i=" << i << " Value=" << hGridCellSizes[i] << " dfdx=" << dfdx[i] << std::endl;
901  }
902  }
903  double maxStepSize = -std::numeric_limits<double>::max();
904  double nmax;
905  for(unsigned int i = 1; i < hGridCellSizes.size(); i++)
906  {
907  nmax = (hGridCellSizes[i - 1] - hGridCellSizes[i]) / (dfdx[i] - dfdx[i - 1]);
908  if(verbosity > 0)
909  {
910  std::cout << "Max f=" << nmax << " because otherwise " << hGridCellSizes[i - 1] << "+x*" << dfdx[i - 1] << ">" << hGridCellSizes[i] << "+x*" << dfdx[i] << std::endl;
911  }
912  if(nmax < 0)
913  {
914  maxStepSize = std::max(nmax, maxStepSize);
915  }
916  }
917  if(verbosity > 0)
918  {
919  std::cout << "maxStepSize=" << maxStepSize << std::endl;
920  }
921  return maxStepSize;
922 }
923 
924 void HGridOptimiser::applyStep(std::vector<double>& hGridCellSizes, std::vector<double>& dfdx, double stepsize, int verbosity)
925 {
926  if(verbosity > 0)
927  {
928  std::cout << "Apply step:" << std::endl;
929  }
930  for(unsigned int i = 0; i < hGridCellSizes.size() - 1; i++)
931  {
932  hGridCellSizes[i] += stepsize * dfdx[i];
933  if(verbosity > 0)
934  {
935  std::cout << "hGridCellSizes[i]" << "+" << stepsize << "*" << dfdx[i] << "=" << hGridCellSizes[i] << std::endl;
936  }
937  }
938 }
939 
940 double HGridOptimiser::goldenSectionSearch(std::vector<double>& startHGridCellSizes, std::vector<double>& searchDirection, double min, double cur, double max, HGridMethod method, int verbosity)
941 {
942  std::vector<double> curHGridCellSizes = startHGridCellSizes;
943  std::vector<double> xHGridCellSizes = startHGridCellSizes;
944 
945  applyStep(curHGridCellSizes, searchDirection, cur, verbosity - 1);
946  double curWork = calculateWork(curHGridCellSizes, method, verbosity - 1);
947 
948  double x;
949  double resphi = 2.0 - 0.5 * (1 + std::sqrt(5));
950 
951  //Determine new probing point
952  if(max - cur > cur - min)
953  {
954  //x between cur and max
955  x = cur + resphi * (max - cur);
956  }
957  else
958  {
959  //x between min and cur
960  x = cur - resphi * (cur - min);
961  }
962 
963  if(std::abs(max - min) < 1e-7 * (std::abs(cur) + std::abs(x)))
964  {
965  return 0.5 * (min + max);
966  }
967 
968  applyStep(xHGridCellSizes, searchDirection, x, 0);
969  double xWork = calculateWork(xHGridCellSizes, method, 0);
970  if(verbosity > 0)
971  {
972  std::cout << "min=" << min << " max=" << max << " cur=" << cur << " curWork=" << curWork << " x=" << x << " xWork=" << xWork << std::endl;
973  }
974  if(xWork < curWork) //x is better
975  {
976  if(max - cur > cur - min)
977  {
978  return goldenSectionSearch(startHGridCellSizes, searchDirection, cur, x, max, method, verbosity);
979  }
980  else
981  {
982  return goldenSectionSearch(startHGridCellSizes, searchDirection, min, x, cur, method, verbosity);
983  }
984  }
985  else //cur is better
986  {
987  if(max - cur > cur - min)
988  {
989  return goldenSectionSearch(startHGridCellSizes, searchDirection, min, cur, x, method, verbosity);
990  }
991  else
992  {
993  return goldenSectionSearch(startHGridCellSizes, searchDirection, x, cur, max, method, verbosity);
994  }
995  }
996 }
997 
998 void HGridOptimiser::getOptimalDistribution(std::vector<double>& hGridCellSizes, unsigned int numberOfLevels, HGridMethod method, int verbosity)
999 {
1000  hGridCellSizes.clear();
1001  for(unsigned int i = 0; i <= numberOfLevels; i++)
1002  {
1003  hGridCellSizes.push_back(2.0 * (rMin_ + (rMax_ - rMin_) * i / numberOfLevels));
1004  }
1005 
1006  std::vector<double> dfdx;
1007  double step, max, W, oW;
1008  W = calculateWork(hGridCellSizes, method, verbosity - 3);
1009  int stepNumber = 0;
1010  do
1011  {
1012  oW = W;
1013  calculateDiffWork(hGridCellSizes, dfdx, method, verbosity - 3);
1014  max = checkLimit(hGridCellSizes, dfdx, verbosity - 3);
1015  step = goldenSectionSearch(hGridCellSizes, dfdx, max, 0.5 * max, 0, method, verbosity - 2);
1016  applyStep(hGridCellSizes, dfdx, step, verbosity - 3);
1017  W = calculateWork(hGridCellSizes, method, verbosity - 3);
1018  while(W > oW)
1019  {
1020  if(verbosity > 1)
1021  {
1022  std::cout << stepNumber << " Bad step, trying closer search range" << std::endl;
1023  }
1024  applyStep(hGridCellSizes, dfdx, -step, verbosity - 3);
1025  max /= 2;
1026  step = goldenSectionSearch(hGridCellSizes, dfdx, max, 0.5 * max, 0, method, verbosity - 2);
1027  applyStep(hGridCellSizes, dfdx, step, verbosity - 3);
1028  W = calculateWork(hGridCellSizes, method, verbosity - 3);
1029  }
1030  ++stepNumber;
1031  if(verbosity > 1)
1032  {
1033  std::cout << stepNumber << " Correct step: old work=" << oW << " new work=" << W << " difference=" << oW - W << " step=" << step / max << std::endl;
1034  }
1035  }
1036  while(oW - W > 1e-13 && stepNumber);
1037  calculateDiffWork(hGridCellSizes, dfdx, method, verbosity - 2);
1038  if(verbosity > 0)
1039  {
1040  std::cout << "Work=" << W << std::endl;
1041  std::cout << "Optimal cell sizes:";
1042  for(unsigned int i = 0; i < hGridCellSizes.size(); i++)
1043  {
1044  std::cout << " " << hGridCellSizes[i];
1045  }
1046  std::cout << std::endl;
1047  }
1048 }
1049 
1050 void HGridOptimiser::histNumberParticlesPerCell(std::vector<double>& hGridCellSizes)
1051 {
1052  unsigned int NLevels = hGridCellSizes.size() - 1;
1053 
1054  std::vector<double> particlesPerLevel;
1055  for(unsigned int i = 0; i < NLevels; i++)
1056  {
1057  particlesPerLevel.push_back(pdfInt(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], 0));
1058  if(particlesPerLevel.back() < 0)
1059  {
1060  particlesPerLevel.back() = 0;
1061  }
1062  }
1063 
1064  std::vector<double> cellsPerLevel;
1065  for(unsigned int i = 0; i < NLevels; i++)
1066  {
1067  cellsPerLevel.push_back(pow(length_ / hGridCellSizes[i + 1], dimension_));
1068  }
1069 
1070  for(unsigned int i = 0; i < NLevels; i++)
1071  {
1072  double l = particlesPerLevel[i] / cellsPerLevel[i];
1073  std::cout << "Histogram for level " << i << ":";
1074  for(unsigned int k = 0; k <= 10; k++)
1075  {
1076  std::cout << " " << std::setw(8) << std::floor(std::pow(l, k) * std::exp(-l) / mathsFunc::factorial(k)*1e4 * cellsPerLevel[i]);
1077  }
1078  std::cout << std::endl;
1079  }
1080 }
double goldenSectionSearch(std::vector< double > &startHGridCellSizes, std::vector< double > &searchDirection, double min, double cur, double max, HGridMethod method, int verbosity)
BaseParticle * getLargestParticle() const
Gets a pointer to the largest BaseParticle (by interactionRadius) in this ParticleHandler.
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 getSystemDimensions() const
Returns the dimension of the simulation. Note there is also a particle dimension. ...
Definition: DPMBase.cc:466
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
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.cc:259
Mdouble getXMin() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMin() returns XMin...
Definition: DPMBase.cc:224
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)
double intCell2Min(unsigned int i)
Mdouble getYMin() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMin() returns YMin...
Definition: DPMBase.cc:238
const std::vector< T * >::const_iterator end() const
Gets the end of the const_iterator over all BaseBoundary in this BaseHandler.
Definition: BaseHandler.h:494
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:482
Mdouble getXMax() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMax() returns XMax...
Definition: DPMBase.cc:231
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:46
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:74
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:124
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.
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:878
double checkLimit(std::vector< double > &hGridCellSizes, std::vector< double > &dfdx, int verbosity)
std::vector< double > cellN_
unsigned int getNumberOfObjects() const
Gets the number of Object in this BaseHandler.
Definition: BaseHandler.h:464
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.cc:245
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.cc:252
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...
Mdouble getInteractionRadius() const
Returns the particle's interaction radius, which might be different from radius_ (e.g., when dealing with wet particles)
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)