multirate RK43 and Group of Elements (minEdge, maxEdge)

Solver/TESTCASES/Advection_1.lua

+model = GModel  ()
+print'*** Loading the geo ***'
+model:load ('square_m.geo')
+-- model:load ('square_p.geo')
+print'*** Loading the mesh ***'
+model:load ('square_m.msh')
+-- model:load('square_p.msh')
+print'*** mesh loaded ***'
+dg = dgSystemOfEquations (model)
+order=3
+dg:setOrder(order)
+print("*** order = " .. order .. " ***")
+print'*** Model set ***'
+
+print("*** set initial condition ***")
+-- initial condition
+function initial_condition( xyz , f )
+  for i=0,xyz:size1()-1 do
+    x = xyz:get(i,0)
+    y = xyz:get(i,1)
+    z = xyz:get(i,2)
+    f:set (i, 0, math.exp(-100*((x+0.2)^2 +(y+0.3)^2)))
+  end
+end
+
+-- conservation law
+-- advection speed
+v=fullMatrix(3,1);
+v:set(0,0,0.01)
+v:set(1,0,0.025)
+v:set(2,0,0)
+--print("*** set advection speed to: v_x=" .. v:get(0,0) ..", v_y=" .. v:get(1,0) .. ", v_z=" .. v:get(2,0) .. "  ***")
+law = dgConservationLawAdvectionDiffusion(functionConstant(v):getName(),'')
+dg:setConservationLaw(law)
+
+print("*** set boundary conditions *** ")
+-- boundary condition
+outside=fullMatrix(1,1)
+outside:set(0,0,0.0)
+bndcondition=law:newOutsideValueBoundary(functionConstant(outside):getName())
+law:addBoundaryCondition('Border',bndcondition)
+
+--dg:setup()
+-- create 2 groups of elements; criterion=maxEdge 
+dg:createGroups("maxEdge")
+-- dg:createGroups("minEdge")
+-- dg:createGroups("elementType")
+
+dg:L2Projection(functionLua(1,'initial_condition',{'XYZ'}):getName())
+
+dg:exportSolution('output/Advection_00000')
+
+dg:computeInvSpectralRadius();
+T=40
+dt=0.1
+
+print("***     default time step (dt)=" .. dt .. ",  Period (T)=" ..T .."    ***")
+
+N=T/dt
+print("***     Number of time steps=" .. N .. "          ***")
+-- main loop
+for i=1,N do
+  norm = dg:multirateRK43(dt)
+    print('iter',i,norm)
+  if (i % 20 == 0) then 
+    dg:exportSolution(string.format("output/Advection_%05d", i)) 
+  end
+end
+print '*** adding the mesh and the model ***'

Solver/TESTCASES/square_m.geo

+Point(1) = {-0.5, -0.5, 0, 0.025};
+Point(2) = {0.5, -0.5, 0, 0.025};
+Point(3) = {0.5, 0.5, 0, 0.125};
+Point(4) = {-0.5, 0.5, 0, 0.125};
+Line(1) = {4, 3};
+Line(2) = {3, 2};
+Line(3) = {2, 1};
+Line(4) = {1, 4};
+Line Loop(5) = {2, 3, 4, 1};
+Plane Surface(6) = {5};
+Physical Line("Border") = {1, 2, 3, 4};
+Physical Surface("Inside") = {6};
+//Mesh.CharacteristicLengthExtendFromBoundary=1;
+//Transfinite Line {1, 2, 4, 3} = 20 Using Progression 1;
+//Transfinite Surface {6};

Solver/TESTCASES/square_p.geo

+Point(1) = {-0.5, -0.5, 0, .03};
+Point(2) = {0, -0.5, 0, .03};
+Point(3) = {0, 0.5, 0, .03};
+Point(4) = {-0.5, 0.5, 0, .03};
+
+Point(5) = {0.5, -0.5, 0, .03};
+Point(6) = {0.5, 0.5, 0, .03};
+
+Line(1) = {4, 3};
+Line(2) = {3, 2};
+Line(3) = {2, 1};
+Line(4) = {1, 4};
+
+Line(7) = {3, 6};
+Line(8) = {6, 5};
+Line(9) = {5, 2};
+
+Line Loop(5) = {2, 3, 4, 1};
+Line Loop(6) = {8, 9, -2, 7};
+Plane Surface(7) = {5};
+Plane Surface(8) = {6};
+Transfinite Line {4,8} = 20 ;
+Transfinite Line{2}=20;
+Transfinite Surface {7};
+Recombine Surface {7};
+Transfinite Line {3,1,7,9} = 10;
+Transfinite Surface{8};
+
+Physical Line("Border") = {1, 3, 4, 7, 8, 9};
+Physical Surface("Inside1") = {7};
+Physical Surface("Inside2") = {8};

Solver/dgAlgorithm.cpp

@@ -357,9 +357,11 @@ void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw,			// conse
 //   double b[4]={1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0};
 //   double c[4]={0, 1.0/2.0, 1.0/2.0, 1};
    
-
+	 double A[10][10];
+	 double *b;
+	 double *c;
 // Small step RK43
-   double A[10][10]={
+   double A1[10][10]={
 {0,         0,         0,         0,         0,         0,         0,         0,         0,         0},
 {1.0/4.0   ,0,         0,         0,         0,         0,         0,         0,         0,         0},
 {-1.0/12.0, 1.0/3.0,   0,         0,         0,         0,         0,         0,         0,         0},
@@ -371,24 +373,24 @@ void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw,			// conse
 {1.0/12.0,  1.0/6.0,   1.0/6.0,   1.0/12.0,  0,         1.0/6.0,   -1.0/6.0,  1.0/2.0,   0,         0},
 {1.0/12.0,  1.0/6.0,   1.0/6.0,   1.0/12.0,  0,         1.0/12.0,  1.0/6.0,   1.0/6.0,   1.0/12.0, 0}
    };
-   double b[10]={1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0,1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0};
-   double c[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
+   double b1[10]={1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0,1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0};
+   double c1[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
 
 // Big step RK43
-//   double A[10][10]={
-//{0,         0,         0,         0,         0,         0,         0,         0,         0,         0},
-//{1.0/4.0 ,  0,         0,         0,         0,         0,         0,         0,         0,         0},
-//{1.0/4.0 ,  0,         0,         0,         0,         0,         0,         0,         0,         0},
-//{1.0/2.0 ,  0,         0,         0,         0,         0,         0,         0,         0,         0},
-//{1.0/2.0 ,  0,         0,         0,         0,         0,         0,         0,         0,         0},
-//{-1.0/6.0,  0,         0,         0,         2.0/3.0,   0,         0,         0,         0,         0},
-//{1.0/12.0,  0,         0,         0,         1.0/6.0,   1.0/2.0,   0,         0,         0,         0},
-//{1.0/12.0,  0,         0,         0,         1.0/6.0,   1.0/2.0,   0,         0,         0,         0},
-//{1.0/3.0 ,  0,         0,         0,         -1.0/3.0,  1,         0,         0,         0,         0},
-//{1.0/3.0 ,  0,         0,         0,         -1.0/3.0,  1,         0,         0,         0,         0},
-//   };
-//   double b[10]={1.0/6.0, 0 , 0 , 0 , 1.0/3.0,1.0/3.0, 0 , 0 , 0 , 1.0/6.0};
-//   double c[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
+   double A2[10][10]={
+{0,         0,         0,         0,         0,         0,         0,         0,         0,         0},
+{1.0/4.0 ,  0,         0,         0,         0,         0,         0,         0,         0,         0},
+{1.0/4.0 ,  0,         0,         0,         0,         0,         0,         0,         0,         0},
+{1.0/2.0 ,  0,         0,         0,         0,         0,         0,         0,         0,         0},
+{1.0/2.0 ,  0,         0,         0,         0,         0,         0,         0,         0,         0},
+{-1.0/6.0,  0,         0,         0,         2.0/3.0,   0,         0,         0,         0,         0},
+{1.0/12.0,  0,         0,         0,         1.0/6.0,   1.0/2.0,   0,         0,         0,         0},
+{1.0/12.0,  0,         0,         0,         1.0/6.0,   1.0/2.0,   0,         0,         0,         0},
+{1.0/3.0 ,  0,         0,         0,         -1.0/3.0,  1,         0,         0,         0,         0},
+{1.0/3.0 ,  0,         0,         0,         -1.0/3.0,  1,         0,         0,         0,         0},
+   };
+   double b2[10]={1.0/6.0, 0 , 0 , 0 , 1.0/3.0,1.0/3.0, 0 , 0 , 0 , 1.0/6.0};
+   double c2[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
 
 
    // fullMatrix<double> K(sol);
@@ -410,11 +412,25 @@ void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw,			// conse
 
    Unp.scale(0.0);
    Unp.axpy(sol);
-
+   // Case of 2 groups: first with big steps, second with small steps	 
    for(int j=0; j<nStages;j++){
      tmp.scale(0.0);
      tmp.axpy(sol);
      for(int k=0;k < groups.getNbElementGroups();k++) {
+		 if (k==0) {
+			 for (int ii=0; ii<10; ii++) {
+				 for (int jj=0; jj<10; jj++) {
+					 A[ii][jj]=A2[ii][jj];
+				 }
+			 }
+		 }
+		 else{
+			 for (int ii=0; ii<10; ii++) {
+				 for (int jj=0; jj<10; jj++) {
+					 A[ii][jj]=A1[ii][jj];
+				 }
+			 }
+		 }
        for(int i=0;i<j;i++){
          if(fabs(A[j][i])>1e-12){
            tmp.getGroupProxy(k).add(K[i]->getGroupProxy(k),h*A[j][i]);
@@ -423,6 +439,15 @@ void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw,			// conse
      }
      this->residual(claw,groups,tmp,resd);
      for(int k=0;k < groups.getNbElementGroups(); k++) {
+		 if (k==0) {
+			 b=b2;
+			 c=c2;
+		 }
+		 else{
+			 b=b1;
+			 c=c1;
+		 }
+		 
        dgGroupOfElements *group = groups.getElementGroup(k);
        int nbNodes = group->getNbNodes();
        for(int i=0;i<group->getNbElements();i++) {

Solver/dgAlgorithm.h

@@ -51,6 +51,8 @@ class dgAlgorithm {
 					  fullMatrix<double> &solution, // the solution 
 					  std::vector<double> & DT // elementary time steps
 					   );
+  // works only (for the moment) for two groups of elements, small and big step
+  // Uses RK43 from Schlegel (10 stages)	
   void multirateRungeKutta (
 		   const dgConservationLaw &claw,
        dgGroupCollection &groups,

Solver/dgGroupOfElements.cpp

@@ -838,6 +838,248 @@ void dgGroupCollection::buildGroups(GModel *model, int dim, int order)
   delete []shiftRecv;
 }
 
+
+void dgGroupCollection::buildGroups(GModel *model, int dim, int order,std::string groupType)
+{
+	_model = model;
+	std::map<const std::string,std::set<MVertex*> > boundaryVertices;
+	std::map<const std::string,std::set<MEdge, Less_Edge> > boundaryEdges;
+	std::map<const std::string,std::set<MFace, Less_Face> > boundaryFaces;
+	std::vector<GEntity*> entities;
+	model->getEntities(entities);
+	std::map<int, std::vector<MElement *> >localElements;
+	std::vector<std::map<int, std::vector<MElement *> > >ghostElements(Msg::GetCommSize()); // [partition][elementType]
+	int nlocal=0;
+	int nghosts=0;
+	std::multimap<MElement*, short> &ghostsCells = model->getGhostCells();
+	for(unsigned int i = 0; i < entities.size(); i++){
+		GEntity *entity = entities[i];
+		if(entity->dim() == dim-1){
+			for(unsigned int j = 0; j < entity->physicals.size(); j++){
+				const std::string physicalName = model->getPhysicalName(entity->dim(), entity->physicals[j]);
+				for (int k = 0; k < entity->getNumMeshElements(); k++) {
+					MElement *element = entity->getMeshElement(k);
+					switch(dim) {
+						case 1:
+							boundaryVertices[physicalName].insert( element->getVertex(0) ); 
+							break;
+						case 2:
+							boundaryEdges[physicalName].insert( element->getEdge(0) );
+							break;
+						case 3:
+							boundaryFaces[physicalName].insert( element->getFace(0));
+							break;
+						default :
+							throw;
+					}
+				}
+			}
+		}else if(entity->dim() == dim){
+			double max_edgeM=-1.e22;
+			double max_edgem=1.e22;
+			double min_edgeM=-1.e22;
+			double min_edgem=1.e22;
+			
+			for (int iel=0; iel<entity->getNumMeshElements(); iel++){
+				MElement *el=entity->getMeshElement(iel);
+				if (el->maxEdge()>max_edgeM) {
+					max_edgeM=el->maxEdge();
+				}
+				if (el->maxEdge()<max_edgem) {
+					max_edgem=el->maxEdge();
+				}
+				if (el->minEdge()>min_edgeM) {
+					min_edgeM=el->minEdge();
+				}
+				if (el->minEdge()<min_edgeM) {
+					min_edgem=el->minEdge();
+				}									
+			}
+			//printf("\n a1=%f, a2=%f \n",(min_edgeM+min_edgem)/2.,(max_edgeM+max_edgem)/2.);
+			for (int iel=0; iel<entity->getNumMeshElements(); iel++){
+				MElement *el=entity->getMeshElement(iel);
+				if(el->getPartition()==Msg::GetCommRank()+1 || el->getPartition()==0){
+					
+					if (groupType=="minEdge") {
+							if (el->minEdge()> (min_edgeM+min_edgem)/2.) {
+							 localElements[0].push_back(el);
+							 }
+							 else{
+							 localElements[1].push_back(el);
+							 }
+					}
+					else if(groupType=="maxEdge"){
+							if (el->maxEdge()> (max_edgeM+max_edgem)/2.) {
+							 localElements[0].push_back(el);
+							}
+							else{
+							 localElements[1].push_back(el);
+							}
+						
+					}
+					else{
+							localElements[el->getType()].push_back(el);
+					}
+
+					
+					nlocal++;
+				}else{
+					std::multimap<MElement*, short>::iterator ghost=ghostsCells.lower_bound(el);
+					while(ghost!= ghostsCells.end() && ghost->first==el){
+						if(abs(ghost->second)-1==Msg::GetCommRank()){
+							ghostElements[el->getPartition()-1][el->getType()].push_back(el);
+							nghosts+=1;
+						}
+						ghost++;
+					}
+				}
+			}
+		}
+	}
+	
+	
+	Msg::Info("%d groups of elements",localElements.size());
+	for (int n=0; n<localElements.size(); n++) {
+		printf("\n **** Group %d: %d elements **** \n", n,localElements[n].size());
+	}
+	// do a group of elements for every element type in the mesh
+	for (std::map<int, std::vector<MElement *> >::iterator it = localElements.begin(); it !=localElements.end() ; ++it){
+		_elementGroups.push_back(new dgGroupOfElements(it->second,order));
+	}
+	
+	
+	for (int i=0;i<_elementGroups.size();i++){
+		// create a group of faces on the boundary for every group ef elements
+		switch(dim) {
+			case 1 : {
+				std::map<const std::string, std::set<MVertex*> >::iterator mapIt;
+				for(mapIt=boundaryVertices.begin(); mapIt!=boundaryVertices.end(); mapIt++) {
+					dgGroupOfFaces *gof = new dgGroupOfFaces(*_elementGroups[i],mapIt->first,order,mapIt->second);
+					if (gof->getNbElements())
+						_boundaryGroups.push_back(gof);
+					else
+						delete gof;
+					break;
+				}
+			}
+			case 2 : {
+				std::map<const std::string, std::set<MEdge, Less_Edge> >::iterator mapIt;
+				for(mapIt=boundaryEdges.begin(); mapIt!=boundaryEdges.end(); mapIt++) {
+					dgGroupOfFaces *gof=new dgGroupOfFaces(*_elementGroups[i],mapIt->first,order,mapIt->second);
+					if(gof->getNbElements())
+						_boundaryGroups.push_back(gof);
+					else
+						delete gof;
+				}
+				break;
+			}
+			case 3 : {
+				std::map<const std::string, std::set<MFace, Less_Face> >::iterator mapIt;
+				for(mapIt=boundaryFaces.begin(); mapIt!=boundaryFaces.end(); mapIt++) {
+					dgGroupOfFaces *gof=new dgGroupOfFaces(*_elementGroups[i],mapIt->first,order,mapIt->second);
+					if(gof->getNbElements())
+						_boundaryGroups.push_back(gof);
+					else
+						delete gof;
+				}
+				break;
+			}
+		}
+		// create a group of faces for every face that is common to elements of the same group 
+		_faceGroups.push_back(new dgGroupOfFaces(*_elementGroups[i],order));
+		// create a groupe of d
+		for (int j=i+1;j<_elementGroups.size();j++){
+			dgGroupOfFaces *gof = new dgGroupOfFaces(*_elementGroups[i],*_elementGroups[j],order);
+			if (gof->getNbElements())
+				_faceGroups.push_back(gof);
+			else
+				delete gof;
+		}
+	}
+	//create ghost groups
+	for(int i=0;i<Msg::GetCommSize();i++){
+		for (std::map<int, std::vector<MElement *> >::iterator it = ghostElements[i].begin(); it !=ghostElements[i].end() ; ++it){
+			_ghostGroups.push_back(new dgGroupOfElements(it->second,order,i));
+		}
+	}
+	//create face group for ghostGroups
+	for (int i=0; i<_ghostGroups.size(); i++){
+		for (int j=0;j<_elementGroups.size();j++){
+			dgGroupOfFaces *gof = new dgGroupOfFaces(*_ghostGroups[i],*_elementGroups[j],order);
+			if (gof->getNbElements())
+				_faceGroups.push_back(gof);
+			else
+				delete gof;
+		}
+	}
+	
+	// build the ghosts structure
+	int *nGhostElements = new int[Msg::GetCommSize()];
+	int *nParentElements = new int[Msg::GetCommSize()];
+	for (int i=0;i<Msg::GetCommSize();i++) {
+		nGhostElements[i]=0;
+	}
+	for (size_t i = 0; i< getNbGhostGroups(); i++){
+		dgGroupOfElements *group = getGhostGroup(i);
+		nGhostElements[group->getGhostPartition()] += group->getNbElements();
+	}
+#ifdef HAVE_MPI
+	MPI_Alltoall(nGhostElements,1,MPI_INT,nParentElements,1,MPI_INT,MPI_COMM_WORLD); 
+#else
+	nParentElements[0]=nGhostElements[0];
+#endif
+	int *shiftSend = new int[Msg::GetCommSize()];
+	int *shiftRecv = new int[Msg::GetCommSize()];
+	int *curShiftSend = new int[Msg::GetCommSize()];
+	int totalSend=0,totalRecv=0;
+	for(int i= 0; i<Msg::GetCommSize();i++){
+		shiftSend[i] = (i==0 ? 0 : shiftSend[i-1]+nGhostElements[i-1]);
+		curShiftSend[i] = shiftSend[i];
+		shiftRecv[i] = (i==0 ? 0 : shiftRecv[i-1]+nParentElements[i-1]);
+		totalSend += nGhostElements[i];
+		totalRecv += nParentElements[i];
+	}
+	
+	int *idSend = new int[totalSend];
+	int *idRecv = new int[totalRecv];
+	for (int i = 0; i<getNbGhostGroups(); i++){
+		dgGroupOfElements *group = getGhostGroup(i);
+		int part = group->getGhostPartition();
+		for (int j=0; j< group->getNbElements() ; j++){
+			idSend[curShiftSend[part]++] = group->getElement(j)->getNum();
+		}
+	}
+#ifdef HAVE_MPI
+	MPI_Alltoallv(idSend,nGhostElements,shiftSend,MPI_INT,idRecv,nParentElements,shiftRecv,MPI_INT,MPI_COMM_WORLD);
+#else
+	memcpy(idRecv,idSend,nParentElements[0]*sizeof(int));
+#endif
+	//create a Map elementNum :: group, position in group
+	std::map<int, std::pair<int,int> > elementMap;
+	for(size_t i = 0; i< getNbElementGroups(); i++) {
+		dgGroupOfElements *group = getElementGroup(i);
+		for(int j=0; j< group->getNbElements(); j++){
+			elementMap[group->getElement(j)->getNum()]=std::pair<int,int>(i,j);
+		}
+	}
+	_elementsToSend.resize(Msg::GetCommSize());
+	for(int i = 0; i<Msg::GetCommSize();i++){
+		_elementsToSend[i].resize(nParentElements[i]);
+		for(int j=0;j<nParentElements[i];j++){
+			int num = idRecv[shiftRecv[i]+j];
+			_elementsToSend[i][j]=elementMap[num];
+		}
+	}
+	delete []idRecv;
+	delete []idSend;
+	delete []curShiftSend;
+	delete []shiftSend;
+	delete []shiftRecv;
+}
+
+
+
+
 dgGroupCollection::~dgGroupCollection() {
   for (int i=0; i< _elementGroups.size(); i++)
     delete _elementGroups[i];

Solver/dgGroupOfElements.h

@@ -208,6 +208,7 @@ class dgGroupCollection {
   inline int getImageElementPositionInGroup(int partId, int i) const {return _elementsToSend[partId][i].second;}
 
   void buildGroups (GModel *model,int dimension, int order);
+  void buildGroups (GModel *model,int dimension, int order, std::string groupType);
   ~dgGroupCollection();
 };
 #endif

Solver/dgSystemOfEquations.cpp

@@ -58,6 +58,9 @@ void dgSystemOfEquations::registerBindings(binding *b) {
   cm->setDescription("set the conservation law this system solves");
   cm = cb->addMethod("setup",&dgSystemOfEquations::setup);
   cm->setDescription("allocate and init internal stuff, call this function after setOrder and setLaw and before anything else on this instance");
+  cm = cb->addMethod("createGroups",&dgSystemOfEquations::createGroups);
+  cm->setArgNames("groupType",NULL);
+  cm->setDescription("allocate and init internal stuff, creates groups form criterion groupType");
   cm = cb->addMethod("exportSolution",&dgSystemOfEquations::exportSolution);
   cm->setArgNames("filename",NULL);
   cm->setDescription("Print the solution into a file. This file does not contain the mesh. To visualize the solution in gmsh you have to open the mesh file first.");
@@ -101,6 +104,16 @@ void dgSystemOfEquations::L2Projection (std::string functionName){
   }
 }
 
+// dgSystemOfEquations::setup() + build groups with criterion:
+// - default: elementType
+// - minEdge (based on minimum edges of elements) , for the moment only two groups possible
+// - maxEdge (based on maximum edges of elements) , for the moment only two groups possible
+void dgSystemOfEquations::createGroups(std::string groupType){
+	_groups.buildGroups(_gm,_dimension,_order,groupType);
+	_solution = new dgDofContainer(_groups,_claw->nbFields());
+	_rightHandSide = new dgDofContainer(_groups,_claw->nbFields());
+	}
+	
 // ok, we can setup the groups and create solution vectors
 void dgSystemOfEquations::setup(){
   if (!_claw) throw;

Solver/dgSystemOfEquations.h

@@ -42,6 +42,7 @@ public:
   void setOrder (int order); // set the polynomial order
   void setConservationLaw (dgConservationLaw *law); // set the conservationLaw
   dgSystemOfEquations(GModel *_gm);
+  void createGroups(std::string groupType); // create groups from criterion: minEdge (2 groups), maxEdge (2 groups), elementType, allocate (same as setup())
   void setup (); // setup the groups and allocate
   void exportSolution (std::string filename); // export the solution
   void limitSolution (); // apply the limiter on the solution