MercuryDPM  Trunk
 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 975 of file HGridOptimiser.cc.

References constants::i.

Referenced by getOptimalDistribution(), and goldenSectionSearch().

977 {
978  if (verbosity > 0)
979  {
980  std::cout << "Apply step:" << std::endl;
981  }
982  for (unsigned int i = 0; i < hGridCellSizes.size() - 1; i++)
983  {
984  hGridCellSizes[i] += stepsize * dfdx[i];
985  if (verbosity > 0)
986  {
987  std::cout << "hGridCellSizes[i]" << "+" << stepsize << "*" << dfdx[i] << "=" << hGridCellSizes[i]
988  << std::endl;
989  }
990  }
991 }
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
void HGridOptimiser::calcDfDx ( std::vector< double > &  hGridCellSizes,
std::vector< double > &  dfdx,
HGridMethod  method,
int  verbosity 
)

Definition at line 886 of file HGridOptimiser.cc.

References calculateWork(), and constants::i.

888 {
889  if (verbosity > 0)
890  {
891  std::cout << "calcDfDx" << std::endl;
892  }
893  double W = calculateWork(hGridCellSizes, method, verbosity - 1);
894  double dx = 1e-7;
895  double nW;
896  dfdx.clear();
897  for (unsigned int i = 0; i < hGridCellSizes.size(); i++)
898  {
899  hGridCellSizes[i] += dx;
900  nW = calculateWork(hGridCellSizes, method, verbosity - 1);
901  dfdx.push_back((nW - W) / dx);
902  hGridCellSizes[i] -= dx;
903  if (verbosity > 0)
904  {
905  std::cout << "i=" << i << " Value=" << hGridCellSizes[i] << " Change=" << nW - W << " dfdx=" << dfdx.back()
906  << std::endl;
907  }
908  }
909 }
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.
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
void HGridOptimiser::calculateDiffWork ( std::vector< double > &  hGridCellSizes,
std::vector< double > &  dfdx,
HGridMethod  method,
int  verbosity 
)

Definition at line 468 of file HGridOptimiser.cc.

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

Referenced by getOptimalDistribution().

470 {
471  unsigned int NLevels = hGridCellSizes.size() - 1;
472  unsigned int NParams = hGridCellSizes.size();
473 
474  if (verbosity > 0)
475  {
476  std::cout << "Length scale=" << length_ << std::endl;
477  }
478  if (verbosity > 0)
479  {
480  std::cout << "Level sizes:";
481  for (unsigned int i = 0; i < NParams; i++)
482  {
483  std::cout << " " << hGridCellSizes[i];
484  }
485  std::cout << std::endl;
486  }
487 
488  std::vector<double> particlesPerLevel;
489  for (unsigned int i = 0; i < NLevels; i++)
490  {
491  particlesPerLevel.push_back(pdfInt(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], 0));
492  if (particlesPerLevel.back() < 0)
493  {
494  particlesPerLevel.back() = 0;
495  }
496  }
497  if (verbosity > 0)
498  {
499  std::cout << "Particles per level: ";
500  for (unsigned int i = 0; i < NLevels; i++)
501  {
502  std::cout << " " << particlesPerLevel[i];
503  }
504  std::cout << std::endl;
505  }
506 
507  //How changes particlesPerLeve[i] when hGridCellSize[j] is changed
508  std::vector<std::vector<double> > diffParticlesPerLevel;
509  diffParticlesPerLevel.resize(NLevels);
510  for (unsigned int i = 0; i < NLevels; i++)
511  {
512  for (unsigned int j = 0; j < NParams; j++)
513  {
514  diffParticlesPerLevel[i].push_back(0.0);
515  if (j == i)
516  {
517  diffParticlesPerLevel[i][j] += -0.5 * diffPdfInt(0.5 * hGridCellSizes[i], 0);
518  }
519  if (j == i + 1)
520  {
521  diffParticlesPerLevel[i][j] += 0.5 * diffPdfInt(0.5 * hGridCellSizes[i + 1], 0);
522  }
523  }
524  }
525  if (verbosity > 0)
526  {
527  std::cout << "Diff Particles per level: " << std::endl;
528  for (unsigned int i = 0; i < NLevels; i++)
529  {
530  for (unsigned int j = 0; j < NParams; j++)
531  {
532  std::cout << " " << std::setw(12) << diffParticlesPerLevel[i][j];
533  }
534  std::cout << std::endl;
535  }
536  }
537 
538  std::vector<double> cellsPerLevel;
539  for (unsigned int i = 0; i < NLevels; i++)
540  {
541  cellsPerLevel.push_back(pow(length_ / hGridCellSizes[i + 1], dimension_));
542  }
543  if (verbosity > 0)
544  {
545  std::cout << "Cells per level: ";
546  for (unsigned int i = 0; i < NLevels; i++)
547  {
548  std::cout << " " << cellsPerLevel[i];
549  }
550  std::cout << std::endl;
551  }
552 
553  //How changes cellsPerLevel[i] when hGridCellSize[j] is changed
554  std::vector<std::vector<double> > diffCellsPerLevel;
555  diffCellsPerLevel.resize(hGridCellSizes.size());
556  for (unsigned int i = 0; i < NLevels; i++)
557  {
558  for (unsigned int j = 0; j < NParams; j++)
559  {
560  if (j == i + 1)
561  {
562  diffCellsPerLevel[i].push_back(
563  -pow(length_ / hGridCellSizes[i + 1], dimension_) * dimension_ / hGridCellSizes[i + 1]);
564  }
565  else
566  {
567  diffCellsPerLevel[i].push_back(0.0);
568  }
569  }
570  }
571  if (verbosity > 0)
572  {
573  std::cout << "Diff Cells per level: " << std::endl;
574  for (unsigned int i = 0; i < NLevels; i++)
575  {
576  for (unsigned int j = 0; j < NParams; j++)
577  {
578  std::cout << " " << std::setw(12) << diffCellsPerLevel[i][j];
579  }
580  std::cout << std::endl;
581  }
582  }
583 
584  std::vector<double> particlesPerCell;
585  for (unsigned int i = 0; i < NLevels; i++)
586  {
587  particlesPerCell.push_back(particlesPerLevel[i] / cellsPerLevel[i]);
588  }
589  if (verbosity > 0)
590  {
591  std::cout << "Particles per cell: ";
592  for (double i : particlesPerCell)
593  {
594  std::cout << " " << i;
595  }
596  std::cout << std::endl;
597  }
598 
599  //How changes particlesPerCell[i] when hGridCellSize[j] is changed
600  std::vector<std::vector<double> > diffParticlesPerCell;
601  diffParticlesPerCell.resize(NLevels);
602  for (unsigned int i = 0; i < NLevels; i++)
603  {
604  for (unsigned int j = 0; j < NParams; j++)
605  {
606  diffParticlesPerCell[i].push_back(
607  (diffParticlesPerLevel[i][j] * cellsPerLevel[i] - particlesPerLevel[i] * diffCellsPerLevel[i][j]) /
608  (cellsPerLevel[i] * cellsPerLevel[i]));
609  }
610  }
611  if (verbosity > 0)
612  {
613  std::cout << "Diff Particles per Cell: " << std::endl;
614  for (unsigned int i = 0; i < NLevels; i++)
615  {
616  for (unsigned int j = 0; j < NParams; j++)
617  {
618  std::cout << " " << std::setw(12) << diffParticlesPerCell[i][j];
619  }
620  std::cout << std::endl;
621  }
622  }
623 
624  double contactWork = 0, overheadWork = 0;
625  std::vector<double> diffContactWork;
626  std::vector<double> diffOverheadWork;
627  diffContactWork.resize(NParams);
628  diffOverheadWork.resize(NParams);
629  for (unsigned int j = 0; j < NParams; j++)
630  {
631  diffContactWork[j] = 0;
632  diffOverheadWork[j] = 0;
633  }
634 
635  for (unsigned int i = 0; i < NLevels; i++)
636  {
637  contactWork += 0.5 * pow(3, dimension_) * particlesPerCell[i] * particlesPerLevel[i];
638  overheadWork += 0.5 * (pow(3, dimension_) + 1) * particlesPerLevel[i];
639  for (unsigned int j = 0; j < NParams; j++)
640  {
641  diffContactWork[j] += 0.5 * pow(3, dimension_) * (diffParticlesPerCell[i][j] * particlesPerLevel[i] +
642  particlesPerCell[i] * diffParticlesPerLevel[i][j]);
643  diffOverheadWork[j] += 0.5 * (pow(3, dimension_) + 1) * diffParticlesPerLevel[i][j];
644  }
645 
646  unsigned int startJ, endJ;
647  switch (method)
648  {
649  case BOTTOMUP:
650  {
651  startJ = i + 1;
652  endJ = NLevels;
653  break;
654  }
655  case TOPDOWN:
656  {
657  startJ = 0;
658  endJ = i;
659  break;
660  }
661  }
662  for (unsigned int j = startJ; j < endJ; j++)
663  {
664  double numberOfCellToCheck;
665  numberOfCellToCheck = expectedCellsIntegral(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1],
666  dimension_, hGridCellSizes[j + 1]);
667 
668  std::vector<double> diffNumberOfCellToCheck;
669  for (unsigned int k = 0; k < NParams; k++)
670  {
671  diffNumberOfCellToCheck.push_back(0.0);
672  if (k == i)
673  {
674  diffNumberOfCellToCheck[k] += 0.5 * diffStartExpectedCellsIntegral(0.5 * hGridCellSizes[i],
675  0.5 * hGridCellSizes[i + 1],
676  dimension_,
677  hGridCellSizes[j + 1]);
678  }
679  if (k == i + 1)
680  {
681  diffNumberOfCellToCheck[k] += 0.5 * diffEndExpectedCellsIntegral(0.5 * hGridCellSizes[i],
682  0.5 * hGridCellSizes[i + 1],
683  dimension_, hGridCellSizes[j + 1]);
684  }
685  if (k == j + 1)
686  {
687  diffNumberOfCellToCheck[k] += diffHExpectedCellsIntegral(0.5 * hGridCellSizes[i],
688  0.5 * hGridCellSizes[i + 1], dimension_,
689  hGridCellSizes[j + 1]);
690  }
691  }
692  contactWork += particlesPerLevel[i] * numberOfCellToCheck * particlesPerCell[j];
693  overheadWork += particlesPerLevel[i] * numberOfCellToCheck;
694  //std::cout<<"i= "<<i<<" j= "<<j<<" NumberOfCellToCheck= "<<numberOfCellToCheck<<" diffNumberOfCellToCheck[2]= "<<diffNumberOfCellToCheck[2]<<std::endl;
695  for (unsigned int k = 0; k < NParams; k++)
696  {
697  diffContactWork[k] += (diffParticlesPerLevel[i][k] * numberOfCellToCheck * particlesPerCell[j] +
698  particlesPerLevel[i] * diffNumberOfCellToCheck[k] * particlesPerCell[j] +
699  particlesPerLevel[i] * numberOfCellToCheck * diffParticlesPerCell[j][k]);
700  diffOverheadWork[k] += (diffParticlesPerLevel[i][k] * numberOfCellToCheck +
701  particlesPerLevel[i] * diffNumberOfCellToCheck[k]);
702  }
703  }
704  }
705  if (verbosity > 0)
706  {
707  std::cout << "Contact work: " << contactWork << std::endl;
708  std::cout << "Overhead work: " << cellCheckOverContactCheckRatio_ * overheadWork << std::endl;
709  std::cout << "Sum work: " << contactWork + cellCheckOverContactCheckRatio_ * overheadWork << std::endl;
710  std::cout << "diff Contactwork: ";
711  for (unsigned int j = 0; j < NParams; j++)
712  {
713  std::cout << " " << diffContactWork[j];
714  }
715  std::cout << std::endl;
716  std::cout << "diff Overheadwork: ";
717  for (unsigned int j = 0; j < NParams; j++)
718  {
719  std::cout << " " << cellCheckOverContactCheckRatio_ * diffOverheadWork[j];
720  }
721  std::cout << std::endl;
722  std::cout << "diff sum: ";
723  for (unsigned int j = 0; j < NParams; j++)
724  {
725  std::cout << " " << diffContactWork[j] + cellCheckOverContactCheckRatio_ * diffOverheadWork[j];
726  }
727  std::cout << std::endl;
728  }
729  dfdx.clear();
730  for (unsigned int j = 0; j < NParams; j++)
731  {
732  dfdx.push_back(diffContactWork[j] + cellCheckOverContactCheckRatio_ * diffOverheadWork[j]);
733  }
734 }
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)
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
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 736 of file HGridOptimiser.cc.

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

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

737 {
738  unsigned int NLevels = hGridCellSizes.size() - 1;
739  unsigned int NParams = hGridCellSizes.size();
740 
741  if (verbosity > 1)
742  {
743  std::cout << "Length scale=" << length_ << std::endl;
744  }
745  if (verbosity > 0)
746  {
747  std::cout << "Level sizes:" << std::endl;
748  for (unsigned int i = 0; i < NParams; i++)
749  {
750  std::cout << std::setw(10) << hGridCellSizes[i] << " ";
751  }
752  std::cout << std::endl;
753  }
754 
755  std::vector<double> particlesPerLevel;
756  for (unsigned int i = 0; i < NLevels; i++)
757  {
758  particlesPerLevel.push_back(pdfInt(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], 0));
759  if (particlesPerLevel.back() < 0)
760  {
761  particlesPerLevel.back() = 0;
762  }
763  }
764  if (verbosity > 0)
765  {
766  std::cout << "Particles per level:" << std::endl;
767  for (unsigned int i = 0; i < NLevels; i++)
768  {
769  std::cout << std::setw(10) << particlesPerLevel[i] << " ";
770  }
771  std::cout << std::endl;
772  }
773 
774  std::vector<double> cellsPerLevel;
775  for (unsigned int i = 0; i < NLevels; i++)
776  {
777  cellsPerLevel.push_back(pow(length_ / hGridCellSizes[i + 1], dimension_));
778  }
779  if (verbosity > 1)
780  {
781  std::cout << "Cells per level:" << std::endl;
782  for (unsigned int i = 0; i < NLevels; i++)
783  {
784  std::cout << std::setw(10) << cellsPerLevel[i] << " ";
785  }
786  std::cout << std::endl;
787  }
788 
789  std::vector<double> particlesPerCell;
790  for (unsigned int i = 0; i < NLevels; i++)
791  {
792  particlesPerCell.push_back(particlesPerLevel[i] / cellsPerLevel[i]);
793  }
794  if (verbosity > 1)
795  {
796  std::cout << "Particles per cell:" << std::endl;
797  for (double i : particlesPerCell)
798  {
799  std::cout << std::setw(10) << i << " ";
800  }
801  std::cout << std::endl;
802  }
803 
804  std::vector<std::vector<double> > contactWork;
805  std::vector<std::vector<double> > overheadWork;
806  contactWork.resize(NLevels);
807  overheadWork.resize(NLevels);
808  for (unsigned int i = 0; i < NLevels; i++)
809  {
810  contactWork[i].resize(NLevels);
811  overheadWork[i].resize(NLevels);
812  for (unsigned int j = 0; j < NLevels; j++)
813  {
814  contactWork[i][j] = 0;
815  overheadWork[i][j] = 0;
816  }
817  }
818 
819  for (unsigned int i = 0; i < NLevels; i++)
820  {
821  contactWork[i][i] += 0.5 * pow(3, dimension_) * particlesPerCell[i] * particlesPerLevel[i];
822  overheadWork[i][i] += 0.5 * (pow(3, dimension_) + 1) * particlesPerLevel[i];
823 
824  unsigned int startJ, endJ;
825  switch (method)
826  {
827  case BOTTOMUP:
828  {
829  startJ = i + 1;
830  endJ = NLevels;
831  break;
832  }
833  case TOPDOWN:
834  {
835  startJ = 0;
836  endJ = i;
837  break;
838  }
839  }
840  for (unsigned int j = startJ; j < endJ; j++)
841  {
842  double numberOfCellToCheck;
843  numberOfCellToCheck = expectedCellsIntegral(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1],
844  dimension_, hGridCellSizes[j + 1]);
845  contactWork[i][j] += particlesPerLevel[i] * numberOfCellToCheck * particlesPerCell[j];
846  overheadWork[i][j] += particlesPerLevel[i] * numberOfCellToCheck;
847 
848  }
849  }
850 
851  double totalContactWork = 0, totalOverheadWork = 0;
852  for (unsigned int i = 0; i < NLevels; i++)
853  {
854  for (unsigned int j = 0; j < NLevels; j++)
855  {
856  totalContactWork += contactWork[i][j];
857  totalOverheadWork += overheadWork[i][j];
858  }
859  }
860 
861  if (verbosity > 0)
862  {
863  std::cout << "Contact work: " << totalContactWork << std::endl;
864  for (unsigned int i = 0; i < NLevels; i++)
865  {
866  for (unsigned int j = 0; j < NLevels; j++)
867  {
868  std::cout << std::setw(10) << contactWork[i][j] << " ";
869  }
870  std::cout << std::endl;
871  }
872  std::cout << "Overhead work: " << totalOverheadWork << std::endl;
873  for (unsigned int i = 0; i < NLevels; i++)
874  {
875  for (unsigned int j = 0; j < NLevels; j++)
876  {
877  std::cout << std::setw(10) << overheadWork[i][j] << " ";
878  }
879  std::cout << std::endl;
880  }
881  std::cout << "Total work: " << totalContactWork + totalOverheadWork << std::endl;
882  }
883  return totalContactWork + cellCheckOverContactCheckRatio_ * totalOverheadWork;
884 }
double length_
The weighted length of the domain.
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
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 213 of file HGridOptimiser.cc.

References numCells_, rMax_, and rMin_.

Referenced by initialise(), and initialisePolyFunc().

214 {
215  return rMin_ + (rMax_ - rMin_) * (i + 1) / numCells_;
216 }
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 205 of file HGridOptimiser.cc.

References numCells_, rMax_, and rMin_.

Referenced by initialise(), and initialisePolyFunc().

206 {
207  return rMin_ + (rMax_ - rMin_) * i / numCells_;
208 }
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 911 of file HGridOptimiser.cc.

References constants::i.

Referenced by getOptimalDistribution().

912 {
913  if (verbosity > 0)
914  {
915  std::cout << "checkLimits part 1 remove levels" << std::endl;
916  }
917 
918  dfdx[0] = 0;
919  dfdx.back() = 0;
920  for (unsigned int i = 1; i < hGridCellSizes.size(); i++)
921  {
922  if (std::abs(hGridCellSizes[i - 1] - hGridCellSizes[i]) < 1e-6)
923  {
924  if (dfdx[i - 1] < dfdx[i])
925  {
926  if (verbosity > 0)
927  {
928  std::cout << "Remove level" << i << std::endl;
929  }
930  if (i > 0.5 * hGridCellSizes.size())
931  {
932  hGridCellSizes.erase(hGridCellSizes.begin() + i - 1);
933  dfdx.erase(dfdx.begin() + i - 1);
934  }
935  else
936  {
937  hGridCellSizes.erase(hGridCellSizes.begin() + i);
938  dfdx.erase(dfdx.begin() + i);
939  }
940  i--;
941  }
942  }
943  }
944 
945  if (verbosity > 0)
946  {
947  std::cout << "checkLimits part 2 calculate limit" << std::endl;
948  for (unsigned int i = 0; i < hGridCellSizes.size(); i++)
949  {
950  std::cout << "i=" << i << " Value=" << hGridCellSizes[i] << " dfdx=" << dfdx[i] << std::endl;
951  }
952  }
953  double maxStepSize = -std::numeric_limits<double>::max();
954  double nmax;
955  for (unsigned int i = 1; i < hGridCellSizes.size(); i++)
956  {
957  nmax = (hGridCellSizes[i - 1] - hGridCellSizes[i]) / (dfdx[i] - dfdx[i - 1]);
958  if (verbosity > 0)
959  {
960  std::cout << "Max f=" << nmax << " because otherwise " << hGridCellSizes[i - 1] << "+x*" << dfdx[i - 1]
961  << ">" << hGridCellSizes[i] << "+x*" << dfdx[i] << std::endl;
962  }
963  if (nmax < 0)
964  {
965  maxStepSize = std::max(nmax, maxStepSize);
966  }
967  }
968  if (verbosity > 0)
969  {
970  std::cout << "maxStepSize=" << maxStepSize << std::endl;
971  }
972  return maxStepSize;
973 }
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
double HGridOptimiser::diffEndExpectedCellsIntegral ( double  start,
double  end,
int  p,
double  h 
)

Definition at line 399 of file HGridOptimiser.cc.

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

Referenced by calculateDiffWork().

400 {
401  unsigned int startCell = radius2IntCell(start);
402  unsigned int endCell = radius2IntCell(end);
403 
404  double numerator = 0;
405  double denominator = 0;
406 
407  double a = (intCellN[endCell + 1] - intCellN[endCell]) / (intCell2Max(endCell) - intCell2Min(endCell));
408  double b = (intCellN[endCell] * intCell2Max(endCell) - intCellN[endCell + 1] * intCell2Min(endCell)) /
409  (intCell2Max(endCell) - intCell2Min(endCell));
410 
411  double diffTeller = (
412  2.0 / h * p * std::pow(2.0 * (end + h) / h, p - 1) * (a * p * end - a * h + a * end + b * p + 2 * b) *
413  (end + h) / ((p + 1) * (p + 2)) +
414  std::pow(2.0 * (end + h) / h, p) * (a * p + a) * (end + h) / ((p + 1) * (p + 2)) +
415  std::pow(2.0 * (end + h) / h, p) * (a * p * end - a * h + a * end + b * p + 2 * b) / ((p + 1) * (p + 2)));
416  double diffNoemer = (a * end + b);
417  if (startCell == endCell)
418  {
419  numerator = expectedCellsIntegralCellNumerator(start, end, startCell, p, h);
420  denominator = expectedCellsIntegralCellDenominator(start, end, startCell);
421  }
422  else
423  {
424  numerator = expectedCellsIntegralCellNumerator(start, intCell2Max(startCell), startCell, p, h);
425  denominator = expectedCellsIntegralCellDenominator(start, intCell2Max(startCell), startCell);
426  for (unsigned int i = startCell + 1; i < endCell; i++)
427  {
430  }
431  numerator += expectedCellsIntegralCellNumerator(intCell2Min(endCell), end, endCell, p, h);
432  denominator += expectedCellsIntegralCellDenominator(intCell2Min(endCell), end, endCell);
433 
434  }
435  //std::cout<<"new: numerator="<<numerator<<" denominator="<<denominator<<" diffTeller="<<diffTeller<<" diffNoemer="<<diffNoemer<<std::endl;
436  return (diffTeller * denominator - numerator * diffNoemer) / (denominator * denominator);
437 }
double intCell2Min(unsigned int i)
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
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 439 of file HGridOptimiser.cc.

References diffHExpectedCellsIntegralCellNumerator(), expectedCellsIntegralCellDenominator(), constants::i, intCell2Max(), intCell2Min(), and radius2IntCell().

Referenced by calculateDiffWork().

440 {
441  unsigned int startCell = radius2IntCell(start);
442  unsigned int endCell = radius2IntCell(end);
443 
444  double denominator = 0;
445  double diffTeller = 0;
446 
447  if (startCell == endCell)
448  {
449  diffTeller = diffHExpectedCellsIntegralCellNumerator(start, end, startCell, p, h);
450  denominator = expectedCellsIntegralCellDenominator(start, end, startCell);
451  }
452  else
453  {
454  diffTeller = diffHExpectedCellsIntegralCellNumerator(start, intCell2Max(startCell), startCell, p, h);
455  denominator = expectedCellsIntegralCellDenominator(start, intCell2Max(startCell), startCell);
456  for (unsigned int i = startCell + 1; i < endCell; i++)
457  {
460  }
461  diffTeller += diffHExpectedCellsIntegralCellNumerator(intCell2Min(endCell), end, endCell, p, h);
462  denominator += expectedCellsIntegralCellDenominator(intCell2Min(endCell), end, endCell);
463  }
464  return diffTeller / denominator;
465 }
double intCell2Min(unsigned int i)
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
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 295 of file HGridOptimiser.cc.

References constants::i, intCell2Max(), intCell2Min(), and intCellN.

Referenced by diffHExpectedCellsIntegral().

296 {
297  double a = (intCellN[i + 1] - intCellN[i]) / (intCell2Max(i) - intCell2Min(i));
298  double b = (intCellN[i] * intCell2Max(i) - intCellN[i + 1] * intCell2Min(i)) / (intCell2Max(i) - intCell2Min(i));
299  double r;
300  r = start;
301  //double min =std::pow(2.0*r/h+2.0,p)*(a*p*r-a*h+a*r+b*p+2*b)*(r+h)/((p+1)*(p+2));
302  double min =
303  (-2.0 * r / h / h * p * std::pow(2.0 * r / h + 2.0, p - 1) * (a * p * r - a * h + a * r + b * p + 2 * b) *
304  (r + h) +
305  -std::pow(2.0 * r / h + 2.0, p) * a * (r + h) +
306  std::pow(2.0 * r / h + 2.0, p) * (a * p * r - a * h + a * r + b * p + 2 * b)) / ((p + 1) * (p + 2));
307  r = end;
308  //double plus=std::pow(2.0*r/h+2.0,p)*(a*p*r-a*h+a*r+b*p+2*b)*(r+h)/((p+1)*(p+2));
309  double plus =
310  (-2.0 * r / h / h * p * std::pow(2.0 * r / h + 2.0, p - 1) * (a * p * r - a * h + a * r + b * p + 2 * b) *
311  (r + h) +
312  -std::pow(2.0 * r / h + 2.0, p) * a * (r + h) +
313  std::pow(2.0 * r / h + 2.0, p) * (a * p * r - a * h + a * r + b * p + 2 * b)) / ((p + 1) * (p + 2));
314  return plus - min;
315 }
double intCell2Min(unsigned int i)
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
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 266 of file HGridOptimiser.cc.

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

Referenced by calculateDiffWork().

267 {
268  unsigned int xCell = radius2IntCell(x);
269  double denominator = 0;
270  for (unsigned int i = 0; i <= numCells_; i++)
271  {
272  denominator += pdfIntCell(intCell2Min(i), intCell2Max(i), i, 0);
273  }
274  double a = (intCellN[xCell + 1] - intCellN[xCell]) / (intCell2Max(xCell) - intCell2Min(xCell));
275  double b = (intCellN[xCell] * intCell2Max(xCell) - intCellN[xCell + 1] * intCell2Min(xCell)) /
276  (intCell2Max(xCell) - intCell2Min(xCell));
277  return (a * std::pow(x, p + 1) + b * std::pow(x, p)) / denominator;
278 }
double intCell2Min(unsigned int i)
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
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 358 of file HGridOptimiser.cc.

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

Referenced by calculateDiffWork().

359 {
360  unsigned int startCell = radius2IntCell(start);
361  unsigned int endCell = radius2IntCell(end);
362 
363  double numerator = 0;
364  double denominator = 0;
365 
366  double a = (intCellN[startCell + 1] - intCellN[startCell]) / (intCell2Max(startCell) - intCell2Min(startCell));
367  double b = (intCellN[startCell] * intCell2Max(startCell) - intCellN[startCell + 1] * intCell2Min(startCell)) /
368  (intCell2Max(startCell) - intCell2Min(startCell));
369 
370  double diffTeller = -(
371  2.0 / h * p * std::pow(2.0 * (start + h) / h, p - 1) * (a * p * start - a * h + a * start + b * p + 2 * b) *
372  (start + h) / ((p + 1) * (p + 2)) +
373  std::pow(2.0 * (start + h) / h, p) * (a * p + a) * (start + h) / ((p + 1) * (p + 2)) +
374  std::pow(2.0 * (start + h) / h, p) * (a * p * start - a * h + a * start + b * p + 2 * b) /
375  ((p + 1) * (p + 2)));
376  double diffNoemer = -(a * start + b);
377  if (startCell == endCell)
378  {
379  numerator = expectedCellsIntegralCellNumerator(start, end, startCell, p, h);
380  denominator = expectedCellsIntegralCellDenominator(start, end, startCell);
381  }
382  else
383  {
384  numerator = expectedCellsIntegralCellNumerator(start, intCell2Max(startCell), startCell, p, h);
385  denominator = expectedCellsIntegralCellDenominator(start, intCell2Max(startCell), startCell);
386  for (unsigned int i = startCell + 1; i < endCell; i++)
387  {
390  }
391  numerator += expectedCellsIntegralCellNumerator(intCell2Min(endCell), end, endCell, p, h);
392  denominator += expectedCellsIntegralCellDenominator(intCell2Min(endCell), end, endCell);
393 
394  }
395  //std::cout<<"new: numerator="<<numerator<<" denominator="<<denominator<<" diffTeller="<<diffTeller<<" diffNoemer="<<diffNoemer<<std::endl;
396  return (diffTeller * denominator - numerator * diffNoemer) / (denominator * denominator);
397 }
double intCell2Min(unsigned int i)
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
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 328 of file HGridOptimiser.cc.

References expectedCellsIntegralCellDenominator(), expectedCellsIntegralCellNumerator(), constants::i, intCell2Max(), intCell2Min(), and radius2IntCell().

Referenced by calculateDiffWork(), and calculateWork().

329 {
330  unsigned int startCell = radius2IntCell(start);
331  unsigned int endCell = radius2IntCell(end);
332 
333  double numerator = 0;
334  double denominator = 0;
335  if (startCell == endCell)
336  {
337  numerator = expectedCellsIntegralCellNumerator(start, end, startCell, p, h);
338  denominator = expectedCellsIntegralCellDenominator(start, end, startCell);
339  }
340  else
341  {
342  numerator = expectedCellsIntegralCellNumerator(start, intCell2Max(startCell), startCell, p, h);
343  denominator = expectedCellsIntegralCellDenominator(start, intCell2Max(startCell), startCell);
344  for (unsigned int i = startCell + 1; i < endCell; i++)
345  {
348  }
349  numerator += expectedCellsIntegralCellNumerator(intCell2Min(endCell), end, endCell, p, h);
350  denominator += expectedCellsIntegralCellDenominator(intCell2Min(endCell), end, endCell);
351 
352  }
353  // std::cout<<"New: numerator="<<numerator<<" denominator="<<denominator<<std::endl;
354  return numerator / denominator;
355 }
double intCell2Min(unsigned int i)
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
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 317 of file HGridOptimiser.cc.

References constants::i, intCell2Max(), intCell2Min(), and intCellN.

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

318 {
319  double a = (intCellN[i + 1] - intCellN[i]) / (intCell2Max(i) - intCell2Min(i));
320  double b = (intCellN[i] * intCell2Max(i) - intCellN[i + 1] * intCell2Min(i)) / (intCell2Max(i) - intCell2Min(i));
321  return (a / 2.0 * (pow(end, 2) - pow(start, 2)) + b * (pow(end, 1) - pow(start, 1)));
322 }
double intCell2Min(unsigned int i)
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
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 280 of file HGridOptimiser.cc.

References constants::i, intCell2Max(), intCell2Min(), and intCellN.

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

281 {
282  double a = (intCellN[i + 1] - intCellN[i]) / (intCell2Max(i) - intCell2Min(i));
283  double b = (intCellN[i] * intCell2Max(i) - intCellN[i + 1] * intCell2Min(i)) / (intCell2Max(i) - intCell2Min(i));
284  double r;
285  r = start;
286  double min = std::pow(2.0 * (r + h) / h, p) * (a * p * r - a * h + a * r + b * p + 2 * b) * (r + h) /
287  ((p + 1) * (p + 2));
288  r = end;
289  double plus = std::pow(2.0 * (r + h) / h, p) * (a * p * r - a * h + a * r + b * p + 2 * b) * (r + h) /
290  ((p + 1) * (p + 2));
291  return plus - min;
292 }
double intCell2Min(unsigned int i)
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
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 1054 of file HGridOptimiser.cc.

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

1056 {
1057  hGridCellSizes.clear();
1058  for (unsigned int i = 0; i <= numberOfLevels; i++)
1059  {
1060  hGridCellSizes.push_back(2.0 * (rMin_ + (rMax_ - rMin_) * i / numberOfLevels));
1061  }
1062 
1063  std::vector<double> dfdx;
1064  double step, max, W, oW;
1065  W = calculateWork(hGridCellSizes, method, verbosity - 3);
1066  int stepNumber = 0;
1067  do
1068  {
1069  oW = W;
1070  calculateDiffWork(hGridCellSizes, dfdx, method, verbosity - 3);
1071  max = checkLimit(hGridCellSizes, dfdx, verbosity - 3);
1072  step = goldenSectionSearch(hGridCellSizes, dfdx, max, 0.5 * max, 0, method, verbosity - 2);
1073  applyStep(hGridCellSizes, dfdx, step, verbosity - 3);
1074  W = calculateWork(hGridCellSizes, method, verbosity - 3);
1075  while (W > oW)
1076  {
1077  if (verbosity > 1)
1078  {
1079  std::cout << stepNumber << " Bad step, trying closer search range" << std::endl;
1080  }
1081  applyStep(hGridCellSizes, dfdx, -step, verbosity - 3);
1082  max /= 2;
1083  step = goldenSectionSearch(hGridCellSizes, dfdx, max, 0.5 * max, 0, method, verbosity - 2);
1084  applyStep(hGridCellSizes, dfdx, step, verbosity - 3);
1085  W = calculateWork(hGridCellSizes, method, verbosity - 3);
1086  }
1087  ++stepNumber;
1088  if (verbosity > 1)
1089  {
1090  std::cout << stepNumber << " Correct step: old work=" << oW << " new work=" << W << " difference=" << oW - W
1091  << " step=" << step / max << std::endl;
1092  }
1093  } while (oW - W > 1e-13 && stepNumber);
1094  calculateDiffWork(hGridCellSizes, dfdx, method, verbosity - 2);
1095  if (verbosity > 0)
1096  {
1097  std::cout << "Work=" << W << std::endl;
1098  std::cout << "Optimal cell sizes:";
1099  for (double hGridCellSize : hGridCellSizes)
1100  {
1101  std::cout << " " << hGridCellSize;
1102  }
1103  std::cout << std::endl;
1104  }
1105 }
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.
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
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 994 of file HGridOptimiser.cc.

References applyStep(), and calculateWork().

Referenced by getOptimalDistribution().

996 {
997  std::vector<double> curHGridCellSizes = startHGridCellSizes;
998  std::vector<double> xHGridCellSizes = startHGridCellSizes;
999 
1000  applyStep(curHGridCellSizes, searchDirection, cur, verbosity - 1);
1001  double curWork = calculateWork(curHGridCellSizes, method, verbosity - 1);
1002 
1003  double x;
1004  double resphi = 2.0 - 0.5 * (1 + std::sqrt(5));
1005 
1006  //Determine new probing point
1007  if (max - cur > cur - min)
1008  {
1009  //x between cur and max
1010  x = cur + resphi * (max - cur);
1011  }
1012  else
1013  {
1014  //x between min and cur
1015  x = cur - resphi * (cur - min);
1016  }
1017 
1018  if (std::abs(max - min) < 1e-7 * (std::abs(cur) + std::abs(x)))
1019  {
1020  return 0.5 * (min + max);
1021  }
1022 
1023  applyStep(xHGridCellSizes, searchDirection, x, 0);
1024  double xWork = calculateWork(xHGridCellSizes, method, 0);
1025  if (verbosity > 0)
1026  {
1027  std::cout << "min=" << min << " max=" << max << " cur=" << cur << " curWork=" << curWork << " x=" << x
1028  << " xWork=" << xWork << std::endl;
1029  }
1030  if (xWork < curWork) //x is better
1031  {
1032  if (max - cur > cur - min)
1033  {
1034  return goldenSectionSearch(startHGridCellSizes, searchDirection, cur, x, max, method, verbosity);
1035  }
1036  else
1037  {
1038  return goldenSectionSearch(startHGridCellSizes, searchDirection, min, x, cur, method, verbosity);
1039  }
1040  }
1041  else //cur is better
1042  {
1043  if (max - cur > cur - min)
1044  {
1045  return goldenSectionSearch(startHGridCellSizes, searchDirection, min, cur, x, method, verbosity);
1046  }
1047  else
1048  {
1049  return goldenSectionSearch(startHGridCellSizes, searchDirection, x, cur, max, method, verbosity);
1050  }
1051  }
1052 }
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 1107 of file HGridOptimiser.cc.

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

1108 {
1109  unsigned int NLevels = hGridCellSizes.size() - 1;
1110 
1111  std::vector<double> particlesPerLevel;
1112  for (unsigned int i = 0; i < NLevels; i++)
1113  {
1114  particlesPerLevel.push_back(pdfInt(0.5 * hGridCellSizes[i], 0.5 * hGridCellSizes[i + 1], 0));
1115  if (particlesPerLevel.back() < 0)
1116  {
1117  particlesPerLevel.back() = 0;
1118  }
1119  }
1120 
1121  std::vector<double> cellsPerLevel;
1122  for (unsigned int i = 0; i < NLevels; i++)
1123  {
1124  cellsPerLevel.push_back(pow(length_ / hGridCellSizes[i + 1], dimension_));
1125  }
1126 
1127  for (unsigned int i = 0; i < NLevels; i++)
1128  {
1129  double l = particlesPerLevel[i] / cellsPerLevel[i];
1130  std::cout << "Histogram for level " << i << ":";
1131  for (unsigned int k = 0; k <= 10; k++)
1132  {
1133  std::cout << " " << std::setw(8)
1134  << std::floor(std::pow(l, k) * std::exp(-l) / mathsFunc::factorial(k) * 1e4 * cellsPerLevel[i]);
1135  }
1136  std::cout << std::endl;
1137  }
1138 }
double length_
The weighted length of the domain.
Mdouble exp(Mdouble Exponent)
Definition: ExtendedMath.cc:84
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
constexpr T factorial(const T t)
factorial function
Definition: ExtendedMath.h:153
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(), ParticleHandler::getLargestParticle(), BaseParticle::getMaxInteractionRadius(), BaseHandler< T >::getSize(), ParticleHandler::getSmallestParticle(), DPMBase::getSystemDimensions(), DPMBase::getXMax(), DPMBase::getXMin(), DPMBase::getYMax(), DPMBase::getYMin(), DPMBase::getZMax(), DPMBase::getZMin(), constants::i, 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 (double& it : cellN_)
41  {
42  it = 0.0;
43  }
44 
45  //Set the minimum and maximum radius of the particles.
47  rMax_ = nextafter(nextafter(problem.particleHandler.getLargestParticle()->getMaxInteractionRadius(),
48  std::numeric_limits<double>::max()), std::numeric_limits<double>::max());
49 
50  //for all particles, add one to the cell in cellN_ which has is associated with
51  //the radius of that particle.
52  for (std::vector<BaseParticle*>::const_iterator it = problem.particleHandler.begin();
53  it != problem.particleHandler.end(); ++it)
54  {
55  cellN_[radius2Cell((*it)->getMaxInteractionRadius())]++;
56  }
57 
58  //assign what the number of particles is.
59  unsigned int numParticles = problem.particleHandler.getSize();
60 
61  intCellN.push_back(1.5 * cellN_[0] - 0.5 * cellN_[1]);
62  for (double& it : cellN_)
63  {
64  intCellN.push_back(it);
65  }
66  intCellN.push_back(1.5 * cellN_[numCells_ - 1] - 0.5 * cellN_[numCells_ - 2]);
67 
68  if (verbosity > 0)
69  {
70  std::cout << "rMin=" << rMin_ << " rMax=" << rMax_ << " NParticles=" << numParticles << std::endl;
71  for (unsigned int i = 0; i < cellN_.size(); i++)
72  {
73  std::cout << "From" << std::setw(8) << cell2Min(i) << " to" << std::setw(8) << cell2Max(i) << " number="
74  << std::setw(5) << cellN_[i] << std::endl;
75  }
76  std::cout << "Domain size [" << problem.getXMin() << "," << problem.getXMax() << "]x[" << problem.getYMin()
77  << "," << problem.getYMax() << "]x[" << problem.getZMin() << "," << problem.getZMax() << "]"
78  << std::endl;
79 
80  for (unsigned int i = 0; i < intCellN.size() - 1; i++)
81  {
82  std::cout << "From" << std::setw(8) << intCell2Min(i) << " to" << std::setw(8) << intCell2Max(i)
83  << " linear from " << std::setw(8) << intCellN[i] << " to " << std::setw(8) << intCellN[i + 1]
84  << std::endl;
85  }
86  /*std::cout<<"["<<cellN_[0];
87  for (unsigned int i=1;i<cellN_.size();i++)
88  {
89  std::cout<<","<<cellN_[i];
90  }
91  std::cout<<"]"<<std::endl;*/
92  }
93 
94  dimension_ = problem.getSystemDimensions();
95  switch (dimension_)
96  {
97  case 1:
98  length_ = pow((problem.getXMax() - problem.getXMin()) / numParticles, 1);
99  break;
100  case 2:
101  length_ = pow(
102  (problem.getXMax() - problem.getXMin()) * (problem.getYMax() - problem.getYMin()) / numParticles,
103  1.0 / 2.0);
104  break;
105  case 3:
106  length_ = pow((problem.getXMax() - problem.getXMin()) * (problem.getYMax() - problem.getYMin()) *
107  (problem.getZMax() - problem.getZMin()) / numParticles, 1.0 / 3.0);
108  break;
109  }
110 }
BaseParticle * getLargestParticle() const
Returns the pointer of the largest particle in the particle handler. When mercury is running in paral...
double cell2Min(unsigned int i)
Computes the left bound of the cell with given ordinal number.
unsigned int getSize() const
Gets the size of the particleHandler (including mpi and periodic particles)
Definition: BaseHandler.h:655
unsigned int getSystemDimensions() const
Returns the system dimensionality.
Definition: DPMBase.cc:1390
Mdouble getMaxInteractionRadius() const
Returns the particle's interaction radius, which might be different from radius_ (e.g., when dealing with wet particles)
Definition: BaseParticle.h:359
Mdouble getZMax() const
If the length of the problem domain in z-direction is ZMax - ZMin, then getZMax() returns ZMax...
Definition: DPMBase.h:617
Mdouble getXMin() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMin() returns XMin...
Definition: DPMBase.h:586
double length_
The weighted length of the domain.
double intCell2Min(unsigned int i)
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
Mdouble getYMin() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMin() returns YMin...
Definition: DPMBase.h:599
const std::vector< T * >::const_iterator end() const
Gets the end of the const_iterator over all BaseBoundary in this BaseHandler.
Definition: BaseHandler.h:704
double cell2Max(unsigned int i)
Computes the right bound of the cell with given ordinal number.
double intCell2Max(unsigned int i)
const std::vector< T * >::const_iterator begin() const
Gets the begin of the const_iterator over all Object in this BaseHandler.
Definition: BaseHandler.h:690
Mdouble getXMax() const
If the length of the problem domain in x-direction is XMax - XMin, then getXMax() returns XMax...
Definition: DPMBase.h:593
BaseParticle * getSmallestParticle() const
Gets a pointer to the smallest BaseParticle (by interactionRadius) in this ParticleHandler of the loc...
double rMax_
Radius of the largest particle, "rounded" to the smallest double that is larger than the radius of ea...
ParticleHandler particleHandler
An object of the class ParticleHandler, contains the pointers to all the particles created...
Definition: DPMBase.h:1329
std::vector< double > cellN_
Mdouble getYMax() const
If the length of the problem domain in y-direction is YMax - YMin, then getYMax() returns XMax...
Definition: DPMBase.h:605
double rMin_
Radius of the smallest particle, "rounded" to the largest double that is smaller than the radius of e...
std::vector< double > intCellN
unsigned int dimension_
The dimension of the system, usually 3, sometimes 2 or 1.
Mdouble getZMin() const
If the length of the problem domain in z-direction is ZMax - ZMin, then getZMin() returns ZMin...
Definition: DPMBase.h:611
double cellCheckOverContactCheckRatio_
The ratio of the time required for a single geometric contact detection over the time required to ret...
unsigned int radius2Cell(double r)
Assigns a BaseParticle of given radius to a certain cell.
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 113 of file HGridOptimiser.cc.

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

114 {
115  numCells_ = numberOfCells;
117  cellN_.resize(numCells_);
118  for (double& it : cellN_)
119  {
120  it = 0;
121  }
122 
123  rMin_ = nextafter(1, 0.0);
124  rMax_ = nextafter(omega, std::numeric_limits<double>::max());
125  for (unsigned int i = 0; i < cellN_.size(); i++)
126  {
127  double start = cell2Min(i);
128  double end = cell2Max(i);
129  for (unsigned int j = 0; j < coeff.size(); j++)
130  {
131  cellN_[i] += coeff[j] / (1.0 + j) * (std::pow(end, j + 1) - std::pow(start, j + 1));
132  }
133  }
134 
135  intCellN.push_back(1.5 * cellN_[0] - 0.5 * cellN_[1]);
136  for (double& it : cellN_)
137  {
138  intCellN.push_back(it);
139  }
140  intCellN.push_back(1.5 * cellN_[numCells_ - 1] - 0.5 * cellN_[numCells_ - 2]);
141 
142  dimension_ = 1;
143  length_ = 1.0;
144 
145  if (verbosity > 0)
146  {
147  std::cout << "rMin=" << rMin_ << " rMax=" << rMax_ << std::endl;
148  for (unsigned int i = 0; i < cellN_.size(); i++)
149  {
150  std::cout << "From" << std::setw(8) << cell2Min(i) << " to" << std::setw(8) << cell2Max(i) << " number="
151  << std::setw(5) << cellN_[i] << std::endl;
152  }
153 
154  for (unsigned int i = 0; i < intCellN.size() - 1; i++)
155  {
156  std::cout << "From" << std::setw(8) << intCell2Min(i) << " to" << std::setw(8) << intCell2Max(i)
157  << " linear from " << std::setw(8) << intCellN[i] << " to " << std::setw(8) << intCellN[i + 1]
158  << std::endl;
159  }
160  }
161 }
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)
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
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 194 of file HGridOptimiser.cc.

References numCells_, rMax_, and rMin_.

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

195 {
196  if (i == numCells_)
197  return rMax_;
198  else
199  return rMin_ + (rMax_ - rMin_) * (i + 0.5) / numCells_;
200 }
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 186 of file HGridOptimiser.cc.

References numCells_, rMax_, and rMin_.

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

187 {
188  if (i == 0)
189  return rMin_;
190  else
191  return rMin_ + (rMax_ - rMin_) * (i - 0.5) / numCells_;
192 }
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 229 of file HGridOptimiser.cc.

References constants::i, intCell2Max(), intCell2Min(), numCells_, pdfIntCell(), and radius2IntCell().

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

230 {
231  //std::cout<<"pdfInit("<<start<<","<<end<<","<<p<<");"<<std::endl;
232  unsigned int startCell = radius2IntCell(start);
233  unsigned int endCell = radius2IntCell(end);
234 
235  double numerator = 0;
236  if (startCell == endCell)
237  {
238  numerator = pdfIntCell(start, end, startCell, p);
239  }
240  else if (endCell < startCell)
241  {
242  numerator = 0;
243  }
244  else
245  {
246  numerator = pdfIntCell(start, intCell2Max(startCell), startCell, p);
247  //std::cout<<"Teller+="<<pdfIntCell(start,intCell2Max(startCell),startCell,p)<<std::endl;
248  for (unsigned int i = startCell + 1; i < endCell; i++)
249  {
250  numerator += pdfIntCell(intCell2Min(i), intCell2Max(i), i, p);
251  //std::cout<<"Teller+="<<pdfIntCell(intCell2Min(i),intCell2Max(i),i,p)<<std::endl;
252  }
253  numerator += pdfIntCell(intCell2Min(endCell), end, endCell, p);
254  //std::cout<<"Teller+="<<pdfIntCell(intCell2Min(endCell),end,endCell,p)<<std::endl;
255  }
256  double denominator = 0;
257  for (unsigned int i = 0; i <= numCells_; i++)
258  {
259  denominator += pdfIntCell(intCell2Min(i), intCell2Max(i), i, 0);
260  //std::cout<<"Noemer+="<<pdfIntCell(intCell2Min(i),intCell2Max(i),i,0)<<std::endl;
261  }
262  return numerator / denominator;
263 }
double intCell2Min(unsigned int i)
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
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 218 of file HGridOptimiser.cc.

References constants::i, intCell2Max(), intCell2Min(), and intCellN.

Referenced by diffPdfInt(), and pdfInt().

219 {
220  double a = (intCellN[i + 1] - intCellN[i]) / (intCell2Max(i) - intCell2Min(i));
221  double b = (intCellN[i] * intCell2Max(i) - intCellN[i + 1] * intCell2Min(i)) / (intCell2Max(i) - intCell2Min(i));
222  return (a / (p + 2) * (std::pow(end, p + 2) - std::pow(start, p + 2)) +
223  b / (p + 1) * (std::pow(end, p + 1) - std::pow(start, p + 1)));
224 }
double intCell2Min(unsigned int i)
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
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 170 of file HGridOptimiser.cc.

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

Referenced by initialise().

171 {
172  if (r < rMin_ || r >= rMax_)
173  logger(ERROR, "The radius % is not in the range [%,%]", r, rMin_, rMax_);
174  unsigned int y = static_cast<unsigned int> (floor(numCells_ * (r - rMin_) / (rMax_ - rMin_)));
175  return y;
176 }
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 178 of file HGridOptimiser.cc.

References numCells_, rMax_, and rMin_.

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

179 {
180  if (r < rMin_ || r >= rMax_)
181  throw;
182  unsigned int y = static_cast<unsigned int> (floor(numCells_ * (r - rMin_) / (rMax_ - rMin_) + 0.5));
183  return y;
184 }
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 202 of file HGridOptimiser.h.

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

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

Definition at line 211 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 206 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 195 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 177 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 187 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 182 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: