MercuryDPM  Beta
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HGridOptimiser Class Reference

#include <HGridOptimiser.h>

Public Member Functions

void initialise (const MercuryBase &problem, unsigned int numberOfCells, int verbosity)
 
void initialisePolyFunc (double omega, std::vector< double > &coeff, unsigned int numberOfCells, int verbosity)
 
unsigned int radius2Cell (double r)
 Assigns a BaseParticle of given radius to a certain cell. More...
 
unsigned int radius2IntCell (double r)
 
double intCell2Min (unsigned int i)
 
double intCell2Max (unsigned int i)
 
double cell2Min (unsigned int i)
 Computes the left bound of the cell with given ordinal number. More...
 
double cell2Max (unsigned int i)
 Computes the right bound of the cell with given ordinal number. More...
 
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. More...
 
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) More...
 
double expectedCellsIntegralCellNumerator (double start, double end, unsigned int i, int p, double h)
 
double diffHExpectedCellsIntegralCellNumerator (double start, double end, unsigned int i, int p, double h)
 
double expectedCellsIntegralCellDenominator (double start, double end, unsigned int i)
 
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 expected number of cells to check at the level with maximum size h for particle radius between start and end. More...
 
double diffStartExpectedCellsIntegral (double start, double end, int p, double h)
 
double diffEndExpectedCellsIntegral (double start, double end, int p, double h)
 
double diffHExpectedCellsIntegral (double start, double end, int p, double h)
 
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. More...
 
void calculateDiffWork (std::vector< double > &hGridCellSizes, std::vector< double > &dfdx, HGridMethod method, int verbosity)
 
void calcDfDx (std::vector< double > &hGridCellSizes, std::vector< double > &dfdx, HGridMethod method, int verbosity)
 
double checkLimit (std::vector< double > &hGridCellSizes, std::vector< double > &dfdx, int verbosity)
 
void applyStep (std::vector< double > &hGridCellSizes, std::vector< double > &dfdx, double stepsize, int verbosity)
 
double goldenSectionSearch (std::vector< double > &startHGridCellSizes, std::vector< double > &searchDirection, double min, double cur, double max, HGridMethod method, int verbosity)
 
void getOptimalDistribution (std::vector< double > &hGridCellSizes, unsigned int numberOfLevels, HGridMethod method, int verbosity)
 
void histNumberParticlesPerCell (std::vector< double > &hGridCellSizes)
 

Private Attributes

unsigned int numCells_
 Number of cells, usually called levels in the HGrid. More...
 
double rMin_
 Radius of the smallest particle, "rounded" to the largest double that is smaller than the radius of each particle. More...
 
double rMax_
 Radius of the largest particle, "rounded" to the smallest double that is larger than the radius of each particle. More...
 
double length_
 The weighted length of the domain. More...
 
double cellCheckOverContactCheckRatio_
 The ratio of the time required for a single geometric contact detection over the time required to retrieve information from a cell. This seems to be only used for checking the effort required by the HGrid, not to compute the cell sizes. More...
 
unsigned int dimension_
 The dimension of the system, usually 3, sometimes 2 or 1. More...
 
std::vector< double > cellN_
 
std::vector< double > intCellN
 

Detailed Description

Todo:
Find out what this class does and document it.

Definition at line 36 of file HGridOptimiser.h.

Member Function Documentation

void HGridOptimiser::applyStep ( std::vector< double > &  hGridCellSizes,
std::vector< double > &  dfdx,
double  stepsize,
int  verbosity 
)

Definition at line 924 of file HGridOptimiser.cc.

Referenced by getOptimalDistribution(), and goldenSectionSearch().

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 }
void HGridOptimiser::calcDfDx ( std::vector< double > &  hGridCellSizes,
std::vector< double > &  dfdx,
HGridMethod  method,
int  verbosity 
)

Definition at line 838 of file HGridOptimiser.cc.

References calculateWork().

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 }
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.
void HGridOptimiser::calculateDiffWork ( std::vector< double > &  hGridCellSizes,
std::vector< double > &  dfdx,
HGridMethod  method,
int  verbosity 
)

Definition at line 437 of file HGridOptimiser.cc.

References BOTTOMUP, cellCheckOverContactCheckRatio_, diffEndExpectedCellsIntegral(), diffHExpectedCellsIntegral(), diffPdfInt(), diffStartExpectedCellsIntegral(), dimension_, expectedCellsIntegral(), length_, pdfInt(), and TOPDOWN.

Referenced by getOptimalDistribution().

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 }
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 length_
The weighted length of the domain.
double diffStartExpectedCellsIntegral(double start, double end, int p, double h)
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...
unsigned int dimension_
The dimension of the system, usually 3, sometimes 2 or 1.
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)
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 diffEndExpectedCellsIntegral(double start, double end, int p, double h)
double HGridOptimiser::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.

Definition at line 689 of file HGridOptimiser.cc.

References BOTTOMUP, cellCheckOverContactCheckRatio_, dimension_, expectedCellsIntegral(), length_, pdfInt(), and TOPDOWN.

Referenced by calcDfDx(), getOptimalDistribution(), and goldenSectionSearch().

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 }
double length_
The weighted length of the domain.
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...
unsigned int dimension_
The dimension of the system, usually 3, sometimes 2 or 1.
double cellCheckOverContactCheckRatio_
The ratio of the time required for a single geometric contact detection over the time required to ret...
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 HGridOptimiser::cell2Max ( unsigned int  i)

Computes the right bound of the cell with given ordinal number.

Computes the right bound of the cell with given ordinal number.

Definition at line 199 of file HGridOptimiser.cc.

References numCells_, rMax_, and rMin_.

Referenced by initialise(), and initialisePolyFunc().

200 {
201  return rMin_ + (rMax_ - rMin_) * (i + 1) / numCells_;
202 }
double rMax_
Radius of the largest particle, "rounded" to the smallest double that is larger than the radius of ea...
double rMin_
Radius of the smallest particle, "rounded" to the largest double that is smaller than the radius of e...
unsigned int numCells_
Number of cells, usually called levels in the HGrid.
double HGridOptimiser::cell2Min ( unsigned int  i)

Computes the left bound of the cell with given ordinal number.

Computes the left bound of the cell with given ordinal number.

Definition at line 191 of file HGridOptimiser.cc.

References numCells_, rMax_, and rMin_.

Referenced by initialise(), and initialisePolyFunc().

192 {
193  return rMin_ + (rMax_ - rMin_) * i / numCells_;
194 }
double rMax_
Radius of the largest particle, "rounded" to the smallest double that is larger than the radius of ea...
double rMin_
Radius of the smallest particle, "rounded" to the largest double that is smaller than the radius of e...
unsigned int numCells_
Number of cells, usually called levels in the HGrid.
double HGridOptimiser::checkLimit ( std::vector< double > &  hGridCellSizes,
std::vector< double > &  dfdx,
int  verbosity 
)

Definition at line 861 of file HGridOptimiser.cc.

Referenced by getOptimalDistribution().

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 }
double HGridOptimiser::diffEndExpectedCellsIntegral ( double  start,
double  end,
int  p,
double  h 
)

Definition at line 372 of file HGridOptimiser.cc.

References expectedCellsIntegralCellDenominator(), expectedCellsIntegralCellNumerator(), intCell2Max(), intCell2Min(), intCellN, and radius2IntCell().

Referenced by calculateDiffWork().

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 }
double intCell2Min(unsigned int i)
double intCell2Max(unsigned int i)
double expectedCellsIntegralCellDenominator(double start, double end, unsigned int i)
std::vector< double > intCellN
double expectedCellsIntegralCellNumerator(double start, double end, unsigned int i, int p, double h)
unsigned int radius2IntCell(double r)
double HGridOptimiser::diffHExpectedCellsIntegral ( double  start,
double  end,
int  p,
double  h 
)

Definition at line 409 of file HGridOptimiser.cc.

References diffHExpectedCellsIntegralCellNumerator(), expectedCellsIntegralCellDenominator(), intCell2Max(), intCell2Min(), and radius2IntCell().

Referenced by calculateDiffWork().

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 }
double intCell2Min(unsigned int i)
double intCell2Max(unsigned int i)
double diffHExpectedCellsIntegralCellNumerator(double start, double end, unsigned int i, int p, double h)
double expectedCellsIntegralCellDenominator(double start, double end, unsigned int i)
unsigned int radius2IntCell(double r)
double HGridOptimiser::diffHExpectedCellsIntegralCellNumerator ( double  start,
double  end,
unsigned int  i,
int  p,
double  h 
)

Definition at line 276 of file HGridOptimiser.cc.

References intCell2Max(), intCell2Min(), and intCellN.

Referenced by diffHExpectedCellsIntegral().

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 }
double intCell2Min(unsigned int i)
double intCell2Max(unsigned int i)
std::vector< double > intCellN
double HGridOptimiser::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)

Definition at line 251 of file HGridOptimiser.cc.

References intCell2Max(), intCell2Min(), intCellN, numCells_, pdfIntCell(), and radius2IntCell().

Referenced by calculateDiffWork().

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 }
double intCell2Min(unsigned int i)
double intCell2Max(unsigned int i)
std::vector< double > intCellN
double pdfIntCell(double start, double end, unsigned int i, int p)
unsigned int numCells_
Number of cells, usually called levels in the HGrid.
unsigned int radius2IntCell(double r)
double HGridOptimiser::diffStartExpectedCellsIntegral ( double  start,
double  end,
int  p,
double  h 
)

Definition at line 335 of file HGridOptimiser.cc.

References expectedCellsIntegralCellDenominator(), expectedCellsIntegralCellNumerator(), intCell2Max(), intCell2Min(), intCellN, and radius2IntCell().

Referenced by calculateDiffWork().

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 }
double intCell2Min(unsigned int i)
double intCell2Max(unsigned int i)
double expectedCellsIntegralCellDenominator(double start, double end, unsigned int i)
std::vector< double > intCellN
double expectedCellsIntegralCellNumerator(double start, double end, unsigned int i, int p, double h)
unsigned int radius2IntCell(double r)
double HGridOptimiser::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 expected number of cells to check at the level with maximum size h for particle radius between start and end.

Definition at line 305 of file HGridOptimiser.cc.

References expectedCellsIntegralCellDenominator(), expectedCellsIntegralCellNumerator(), intCell2Max(), intCell2Min(), and radius2IntCell().

Referenced by calculateDiffWork(), and calculateWork().

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 }
double intCell2Min(unsigned int i)
double intCell2Max(unsigned int i)
double expectedCellsIntegralCellDenominator(double start, double end, unsigned int i)
double expectedCellsIntegralCellNumerator(double start, double end, unsigned int i, int p, double h)
unsigned int radius2IntCell(double r)
double HGridOptimiser::expectedCellsIntegralCellDenominator ( double  start,
double  end,
unsigned int  i 
)

Definition at line 294 of file HGridOptimiser.cc.

References intCell2Max(), intCell2Min(), and intCellN.

Referenced by diffEndExpectedCellsIntegral(), diffHExpectedCellsIntegral(), diffStartExpectedCellsIntegral(), and expectedCellsIntegral().

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 }
double intCell2Min(unsigned int i)
double intCell2Max(unsigned int i)
std::vector< double > intCellN
double HGridOptimiser::expectedCellsIntegralCellNumerator ( double  start,
double  end,
unsigned int  i,
int  p,
double  h 
)

Definition at line 264 of file HGridOptimiser.cc.

References intCell2Max(), intCell2Min(), and intCellN.

Referenced by diffEndExpectedCellsIntegral(), diffStartExpectedCellsIntegral(), and expectedCellsIntegral().

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 }
double intCell2Min(unsigned int i)
double intCell2Max(unsigned int i)
std::vector< double > intCellN
void HGridOptimiser::getOptimalDistribution ( std::vector< double > &  hGridCellSizes,
unsigned int  numberOfLevels,
HGridMethod  method,
int  verbosity 
)

Definition at line 998 of file HGridOptimiser.cc.

References applyStep(), calculateDiffWork(), calculateWork(), checkLimit(), goldenSectionSearch(), rMax_, and rMin_.

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 }
double goldenSectionSearch(std::vector< double > &startHGridCellSizes, std::vector< double > &searchDirection, double min, double cur, double max, HGridMethod method, int verbosity)
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.
void calculateDiffWork(std::vector< double > &hGridCellSizes, std::vector< double > &dfdx, HGridMethod method, int verbosity)
double rMax_
Radius of the largest particle, "rounded" to the smallest double that is larger than the radius of ea...
double checkLimit(std::vector< double > &hGridCellSizes, std::vector< double > &dfdx, int verbosity)
void applyStep(std::vector< double > &hGridCellSizes, std::vector< double > &dfdx, double stepsize, int verbosity)
double rMin_
Radius of the smallest particle, "rounded" to the largest double that is smaller than the radius of e...
double HGridOptimiser::goldenSectionSearch ( std::vector< double > &  startHGridCellSizes,
std::vector< double > &  searchDirection,
double  min,
double  cur,
double  max,
HGridMethod  method,
int  verbosity 
)

Definition at line 940 of file HGridOptimiser.cc.

References applyStep(), and calculateWork().

Referenced by getOptimalDistribution().

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 }
double goldenSectionSearch(std::vector< double > &startHGridCellSizes, std::vector< double > &searchDirection, double min, double cur, double max, HGridMethod method, int verbosity)
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.
void applyStep(std::vector< double > &hGridCellSizes, std::vector< double > &dfdx, double stepsize, int verbosity)
void HGridOptimiser::histNumberParticlesPerCell ( std::vector< double > &  hGridCellSizes)

Definition at line 1050 of file HGridOptimiser.cc.

References dimension_, mathsFunc::factorial(), length_, and pdfInt().

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 length_
The weighted length of the domain.
constexpr T factorial(const T t)
factorial function
Definition: ExtendedMath.h:124
unsigned int dimension_
The dimension of the system, usually 3, sometimes 2 or 1.
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...
void HGridOptimiser::initialise ( const MercuryBase problem,
unsigned int  numberOfCells,
int  verbosity 
)

Definition at line 33 of file HGridOptimiser.cc.

References BaseHandler< T >::begin(), cell2Max(), cell2Min(), cellCheckOverContactCheckRatio_, cellN_, dimension_, BaseHandler< T >::end(), BaseParticle::getInteractionRadius(), ParticleHandler::getLargestParticle(), BaseHandler< T >::getNumberOfObjects(), ParticleHandler::getSmallestParticle(), DPMBase::getSystemDimensions(), DPMBase::getXMax(), DPMBase::getXMin(), DPMBase::getYMax(), DPMBase::getYMin(), DPMBase::getZMax(), DPMBase::getZMin(), intCell2Max(), intCell2Min(), intCellN, length_, numCells_, DPMBase::particleHandler, radius2Cell(), rMax_, and rMin_.

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 }
BaseParticle * getLargestParticle() const
Gets a pointer to the largest BaseParticle (by interactionRadius) in this ParticleHandler.
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
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 length_
The weighted length of the domain.
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
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...
ParticleHandler particleHandler
An object of the class ParticleHandler, contains the pointers to all the particles created...
Definition: DPMBase.h:878
std::vector< double > cellN_
unsigned int getNumberOfObjects() const
Gets the number of Object in this BaseHandler.
Definition: BaseHandler.h:464
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...
unsigned int radius2Cell(double r)
Assigns a BaseParticle of given radius to a certain cell.
Mdouble getInteractionRadius() const
Returns the particle's interaction radius, which might be different from radius_ (e.g., when dealing with wet particles)
unsigned int numCells_
Number of cells, usually called levels in the HGrid.
void HGridOptimiser::initialisePolyFunc ( double  omega,
std::vector< double > &  coeff,
unsigned int  numberOfCells,
int  verbosity 
)

Definition at line 102 of file HGridOptimiser.cc.

References cell2Max(), cell2Min(), cellCheckOverContactCheckRatio_, cellN_, dimension_, intCell2Max(), intCell2Min(), intCellN, length_, numCells_, rMax_, and rMin_.

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 }
double cell2Min(unsigned int i)
Computes the left bound of the cell with given ordinal number.
double length_
The weighted length of the domain.
double intCell2Min(unsigned int i)
double cell2Max(unsigned int i)
Computes the right bound of the cell with given ordinal number.
double intCell2Max(unsigned int i)
double rMax_
Radius of the largest particle, "rounded" to the smallest double that is larger than the radius of ea...
std::vector< double > cellN_
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.
double cellCheckOverContactCheckRatio_
The ratio of the time required for a single geometric contact detection over the time required to ret...
unsigned int numCells_
Number of cells, usually called levels in the HGrid.
double HGridOptimiser::intCell2Max ( unsigned int  i)

Definition at line 180 of file HGridOptimiser.cc.

References numCells_, rMax_, and rMin_.

Referenced by diffEndExpectedCellsIntegral(), diffHExpectedCellsIntegral(), diffHExpectedCellsIntegralCellNumerator(), diffPdfInt(), diffStartExpectedCellsIntegral(), expectedCellsIntegral(), expectedCellsIntegralCellDenominator(), expectedCellsIntegralCellNumerator(), initialise(), initialisePolyFunc(), pdfInt(), and pdfIntCell().

181 {
182  if(i == numCells_)
183  return rMax_;
184  else
185  return rMin_ + (rMax_ - rMin_) * (i + 0.5) / numCells_;
186 }
double rMax_
Radius of the largest particle, "rounded" to the smallest double that is larger than the radius of ea...
double rMin_
Radius of the smallest particle, "rounded" to the largest double that is smaller than the radius of e...
unsigned int numCells_
Number of cells, usually called levels in the HGrid.
double HGridOptimiser::intCell2Min ( unsigned int  i)

Definition at line 172 of file HGridOptimiser.cc.

References numCells_, rMax_, and rMin_.

Referenced by diffEndExpectedCellsIntegral(), diffHExpectedCellsIntegral(), diffHExpectedCellsIntegralCellNumerator(), diffPdfInt(), diffStartExpectedCellsIntegral(), expectedCellsIntegral(), expectedCellsIntegralCellDenominator(), expectedCellsIntegralCellNumerator(), initialise(), initialisePolyFunc(), pdfInt(), and pdfIntCell().

173 {
174  if(i == 0)
175  return rMin_;
176  else
177  return rMin_ + (rMax_ - rMin_) * (i - 0.5) / numCells_;
178 }
double rMax_
Radius of the largest particle, "rounded" to the smallest double that is larger than the radius of ea...
double rMin_
Radius of the smallest particle, "rounded" to the largest double that is smaller than the radius of e...
unsigned int numCells_
Number of cells, usually called levels in the HGrid.
double HGridOptimiser::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.

Definition at line 214 of file HGridOptimiser.cc.

References intCell2Max(), intCell2Min(), numCells_, pdfIntCell(), and radius2IntCell().

Referenced by calculateDiffWork(), calculateWork(), and histNumberParticlesPerCell().

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 }
double intCell2Min(unsigned int i)
double intCell2Max(unsigned int i)
double pdfIntCell(double start, double end, unsigned int i, int p)
unsigned int numCells_
Number of cells, usually called levels in the HGrid.
unsigned int radius2IntCell(double r)
double HGridOptimiser::pdfIntCell ( double  start,
double  end,
unsigned int  i,
int  p 
)

Definition at line 204 of file HGridOptimiser.cc.

References intCell2Max(), intCell2Min(), and intCellN.

Referenced by diffPdfInt(), and pdfInt().

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 }
double intCell2Min(unsigned int i)
double intCell2Max(unsigned int i)
std::vector< double > intCellN
unsigned int HGridOptimiser::radius2Cell ( double  r)

Assigns a BaseParticle of given radius to a certain cell.

Assigns a cell to a BaseParticle with the given radius. Note that the index of the cells are linear in the radius of a particle. For example, if numCells_ = 10 and we are looking in the radius range [0,10], than the cell number of the particle with radius r is the number r rounded down to an integer.

Definition at line 156 of file HGridOptimiser.cc.

References ERROR, logger, numCells_, rMax_, and rMin_.

Referenced by initialise().

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 }
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double rMax_
Radius of the largest particle, "rounded" to the smallest double that is larger than the radius of ea...
double rMin_
Radius of the smallest particle, "rounded" to the largest double that is smaller than the radius of e...
unsigned int numCells_
Number of cells, usually called levels in the HGrid.
unsigned int HGridOptimiser::radius2IntCell ( double  r)

Definition at line 164 of file HGridOptimiser.cc.

References numCells_, rMax_, and rMin_.

Referenced by diffEndExpectedCellsIntegral(), diffHExpectedCellsIntegral(), diffPdfInt(), diffStartExpectedCellsIntegral(), expectedCellsIntegral(), and pdfInt().

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 }
double rMax_
Radius of the largest particle, "rounded" to the smallest double that is larger than the radius of ea...
double rMin_
Radius of the smallest particle, "rounded" to the largest double that is smaller than the radius of e...
unsigned int numCells_
Number of cells, usually called levels in the HGrid.

Member Data Documentation

double HGridOptimiser::cellCheckOverContactCheckRatio_
private

The ratio of the time required for a single geometric contact detection over the time required to retrieve information from a cell. This seems to be only used for checking the effort required by the HGrid, not to compute the cell sizes.

Definition at line 198 of file HGridOptimiser.h.

Referenced by calculateDiffWork(), calculateWork(), initialise(), and initialisePolyFunc().

std::vector<double> HGridOptimiser::cellN_
private

Definition at line 207 of file HGridOptimiser.h.

Referenced by initialise(), and initialisePolyFunc().

unsigned int HGridOptimiser::dimension_
private

The dimension of the system, usually 3, sometimes 2 or 1.

Definition at line 202 of file HGridOptimiser.h.

Referenced by calculateDiffWork(), calculateWork(), histNumberParticlesPerCell(), initialise(), and initialisePolyFunc().

double HGridOptimiser::length_
private

The weighted length of the domain.

The weighted length is computed by multiplying the lengths of all directions with each other, and then taking the appropriate type of root so that the unit is the same as that of a length (square root for 2D, cube root for 3D).

Definition at line 191 of file HGridOptimiser.h.

Referenced by calculateDiffWork(), calculateWork(), histNumberParticlesPerCell(), initialise(), and initialisePolyFunc().

unsigned int HGridOptimiser::numCells_
private

Number of cells, usually called levels in the HGrid.

Definition at line 173 of file HGridOptimiser.h.

Referenced by cell2Max(), cell2Min(), diffPdfInt(), initialise(), initialisePolyFunc(), intCell2Max(), intCell2Min(), pdfInt(), radius2Cell(), and radius2IntCell().

double HGridOptimiser::rMax_
private

Radius of the largest particle, "rounded" to the smallest double that is larger than the radius of each particle.

Definition at line 183 of file HGridOptimiser.h.

Referenced by cell2Max(), cell2Min(), getOptimalDistribution(), initialise(), initialisePolyFunc(), intCell2Max(), intCell2Min(), radius2Cell(), and radius2IntCell().

double HGridOptimiser::rMin_
private

Radius of the smallest particle, "rounded" to the largest double that is smaller than the radius of each particle.

Definition at line 178 of file HGridOptimiser.h.

Referenced by cell2Max(), cell2Min(), getOptimalDistribution(), initialise(), initialisePolyFunc(), intCell2Max(), intCell2Min(), radius2Cell(), and radius2IntCell().


The documentation for this class was generated from the following files: