Commit 131f0f6f authored by Tomas Sobotik's avatar Tomas Sobotik
Browse files

separation of synchronization cycles and few other changes

parent 08a9c8bf

examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver.h

+1 −1
Original line number Diff line number Diff line
@@ -70,7 +70,7 @@ private: 


	void insertSubgrid( VectorType u, const int i );



	VectorType runSubgrid( int boundaryCondition, VectorType u);

	VectorType runSubgrid( int boundaryCondition, VectorType u, int subGridID);





	VectorType u0, work_u;

examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h

+232 −167
Original line number Diff line number Diff line
@@ -30,26 +30,11 @@ template< typename Scheme> 
void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::test()

{

/*

	VectorType savedMesh;

	savedMesh.setSize(this->work_u.getSize());

	savedMesh = this->work_u;

	for(int i =0; i < this->subgridValues.getSize(); i++ )

	{

		insertSubgrid(getSubgrid(i), i);

		for(int j = 0; j < this->work_u.getSize(); j++)

		{

		if(savedMesh[j] != this->work_u[j])

			cout << "Error on subMesh " << i << " : " << j << endl;

		}

		this->work_u = savedMesh;

	}

*/



/*

	contractGrid();

	this->u0.save("u-00000.tnl");

	stretchGrid();

*/

}



template< typename Scheme>
@@ -68,7 +53,6 @@ bool tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::init( const tnlPara 


	this->subMesh.save("submesh.tnl");





	const tnlString& initialCondition = parameters.GetParameter <tnlString>("initial-condition");

	this->u0.load( initialCondition );
@@ -88,8 +72,8 @@ bool tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::init( const tnlPara 
	cout << "Setting CFL to " << this->cflCondition << endl;



	stretchGrid();

	//this->stopTime /= (double)(this->gridCols);

	//this->stopTime *= (1.0+1.0/((double)(this->n) - 1.0));

	this->stopTime /= (double)(this->gridCols);

	this->stopTime *= 2.0*(1.0+1.0/((double)(this->n) - 1.0));

	cout << "Setting stopping time to " << this->stopTime << endl;



	cout << "Initializating scheme..." << endl;
@@ -125,7 +109,7 @@ bool tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::init( const tnlPara 
		if(containsCurve)

		{

			//cout << "Computing initial SDF on subgrid " << i << "." << endl;

			insertSubgrid(runSubgrid(0, tmp[i]), i);

			insertSubgrid(runSubgrid(0, tmp[i],i), i);

			setSubgridValue(i, 4);

			//cout << "Computed initial SDF on subgrid " << i  << "." << endl;

		}
@@ -136,39 +120,46 @@ bool tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::init( const tnlPara 
	this->currentStep = 1;

	synchronize();

	cout << "Solver initialized." << endl;



	return true;

}



template< typename Scheme >

void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::run()

{

	while (this->boundaryConditions.max() > 0  && this->currentStep < 3.0*this->gridCols)

	bool end = false;

	while ((this->boundaryConditions.max() > 0 ) || !end)

	{

		if(this->boundaryConditions.max() == 0 )

			end=true;

		else

			end=false;



		for(int i = 0; i < this->subgridValues.getSize(); i++)

		{

			if(getSubgridValue(i) != INT_MAX)

			{

				//cout << "subMesh: " << i << ", BC: " << getBoundaryCondition(i) << endl;



				if(getBoundaryCondition(i) & 1)

				{

					insertSubgrid( runSubgrid(1, getSubgrid(i)), i);

					insertSubgrid( runSubgrid(1, getSubgrid(i),i), i);

				}

				if(getBoundaryCondition(i) & 2)

				{

					insertSubgrid( runSubgrid(2, getSubgrid(i)), i);

					insertSubgrid( runSubgrid(2, getSubgrid(i),i), i);

				}

				if(getBoundaryCondition(i) & 4)

				{

					insertSubgrid( runSubgrid(4, getSubgrid(i)), i);

					insertSubgrid( runSubgrid(4, getSubgrid(i),i), i);

				}

				if(getBoundaryCondition(i) & 8)

				{

					insertSubgrid( runSubgrid(8, getSubgrid(i)), i);

					insertSubgrid( runSubgrid(8, getSubgrid(i),i), i);

				}

				/*if(getBoundaryCondition(i))

				{

					insertSubgrid( runSubgrid(15, getSubgrid(i)), i);

					insertSubgrid( runSubgrid(15, getSubgrid(i),i), i);

				}*/



				setBoundaryCondition(i, 0);
@@ -195,6 +186,8 @@ void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::synchronize() //nee 
	int tmp1, tmp2;

	int grid1, grid2;



	if(this->currentStep & 1)

	{

		for(int j = 0; j < this->gridRows - 1; j++)

		{

			for (int i = 0; i < this->gridCols*this->n; i++)
@@ -203,6 +196,8 @@ void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::synchronize() //nee 
				tmp2 = this->gridCols*this->n*((this->n)+j*this->n) + i;

				grid1 = getSubgridValue(getOwner(tmp1));

				grid2 = getSubgridValue(getOwner(tmp2));

				if(getOwner(tmp1)==getOwner(tmp2))

					cout << "i, j" << i << "," << j << endl;

				if ((fabs(this->work_u[tmp1]) < fabs(this->work_u[tmp2]) - this->delta || grid2 == INT_MAX || grid2 == -INT_MAX) && (grid1 != INT_MAX && grid1 != -INT_MAX))

				{

					this->work_u[tmp2] = this->work_u[tmp1];
@@ -226,10 +221,11 @@ void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::synchronize() //nee 
						setBoundaryCondition(getOwner(tmp1), getBoundaryCondition(getOwner(tmp1))+1);

				}

			}



		}





	}

	else

	{

		for(int i = 1; i < this->gridCols; i++)

		{

			for (int j = 0; j < this->gridRows*this->n; j++)
@@ -238,6 +234,8 @@ void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::synchronize() //nee 
				tmp2 = this->gridCols*this->n*j + i*this->n ;

				grid1 = getSubgridValue(getOwner(tmp1));

				grid2 = getSubgridValue(getOwner(tmp2));

				if(getOwner(tmp1)==getOwner(tmp2))

					cout << "i, j" << i << "," << j << endl;

				if ((fabs(this->work_u[tmp1]) < fabs(this->work_u[tmp2]) - this->delta || grid2 == INT_MAX || grid2 == -INT_MAX) && (grid1 != INT_MAX && grid1 != -INT_MAX))

				{

					this->work_u[tmp2] = this->work_u[tmp1];
@@ -261,9 +259,10 @@ void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::synchronize() //nee 
						setBoundaryCondition(getOwner(tmp1), getBoundaryCondition(getOwner(tmp1))+2);

				}

			}



		}

	}





	this->currentStep++;

	int stepValue = this->currentStep + 4;

	for (int i = 0; i < this->subgridValues.getSize(); i++)
@@ -314,8 +313,11 @@ void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::stretchGrid() 
	cout << "Stretching grid..." << endl;





	this->gridCols = ceil( ((double)(this->mesh.getDimensions().x()-1)) / ((double)(this->n-1)) );

	this->gridRows = ceil( ((double)(this->mesh.getDimensions().y()-1)) / ((double)(this->n-1)) );

	//this->gridCols = ceil( ((double)(this->mesh.getDimensions().x()-1)) / ((double)(this->n-1)) );

	//this->gridRows = ceil( ((double)(this->mesh.getDimensions().y()-1)) / ((double)(this->n-1)) );



	this->gridCols = (this->mesh.getDimensions().x()-1) / (this->n-1) ;

	this->gridRows = (this->mesh.getDimensions().y()-1) / (this->n-1) ;



	cout << "Setting gridCols to " << this->gridCols << "." << endl;

	cout << "Setting gridRows to " << this->gridRows << "." << endl;
@@ -336,11 +338,11 @@ void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::stretchGrid() 
	if(!this->unusedCell.setSize(stretchedSize))

		cerr << "Could not allocate memory for supporting stretched grid." << endl;





	for(int i = 0; i < stretchedSize; i++)

	{

		this->unusedCell[i] = 1;

		int k = i/this->n - i/(this->n*this->gridCols) + this->mesh.getDimensions().x()*(i/(this->n*this->n*this->gridCols));

		//cout << "i = " << i << " : i-k = " << i-k << endl;

		/*int j=(i % (this->n*this->gridCols)) - ( (this->mesh.getDimensions().x() - this->n)/(this->n - 1) + this->mesh.getDimensions().x() - 1)

				+ (this->n*this->gridCols - this->mesh.getDimensions().x())*(i/(this->n*this->n*this->gridCols)) ;
@@ -354,6 +356,7 @@ void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::stretchGrid() 
			k+= l + ( (l / this->mesh.getDimensions().x()) + 1 )*this->mesh.getDimensions().x() - (l % this->mesh.getDimensions().x());*/



		this->work_u[i] = this->u0[i-k];

		//cout << (i-k) <<endl;

	}
@@ -417,7 +420,7 @@ void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::insertSubgrid( Vect 


template< typename Scheme >

typename tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::VectorType

tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::runSubgrid( int boundaryCondition, VectorType u)

tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::runSubgrid( int boundaryCondition, VectorType u, int subGridID)

{



	VectorType fu;
@@ -429,123 +432,189 @@ tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::runSubgrid( int boundary 
 *          Insert Euler-Solver Here

 */



   double time = 0.0;

   double currentTau = this->tau0;

   //bool skip = false;

   if( time + currentTau > this -> stopTime ) currentTau = this -> stopTime - time;

   /*this->resetIterations();

   this->setResidue( this->getConvergenceResidue() + 1.0 );



   this -> refreshSolverMonitor();*/

	/**/



   /****

    * Start the main loop

    */

   //cout << time << flush;

   double maxResidue( 0.0 );

   while( time < this->stopTime/* || maxResidue > (10.0* this -> cflCondition)*/)

	/*for(int i = 0; i < u.getSize(); i++)

	{

      /****

       * Compute the RHS

       */

      //this->problem->getExplicitRHS( time, currentTau, u, fu );

		int x = this->subMesh.getCellCoordinates(i).x();

		int y = this->subMesh.getCellCoordinates(i).y();



      for( int i = 0; i < fu.getSize(); i ++ )

		if(x == 0 && (boundaryCondition & 4) && y ==0)

		{

    	  fu[ i ] = scheme.getValue( this->subMesh, i, this->subMesh.getCellCoordinates(i), u, time, boundaryCondition );

			if((u[subMesh.getCellYSuccessor( i )] - u[i])/subMesh.getHy() > 1.0)

			{

				//cout << "x = 0; y = 0" << endl;

				u[i] = u[subMesh.getCellYSuccessor( i )] - subMesh.getHy();

			}

		}

		else if(x == 0 && (boundaryCondition & 4) && y == subMesh.getDimensions().y() - 1)

		{

			if((u[subMesh.getCellYPredecessor( i )] - u[i])/subMesh.getHy() > 1.0)

			{

				//cout << "x = 0; y = n" << endl;

				u[i] = u[subMesh.getCellYPredecessor( i )] - subMesh.getHy();

			}

		}





      //RealType lastResidue = this->getResidue();

		else if(x == subMesh.getDimensions().x() - 1 && (boundaryCondition & 2) && y ==0)

		{

			if((u[subMesh.getCellYSuccessor( i )] - u[i])/subMesh.getHy() > 1.0)

			{

				//cout << "x = n; y = 0" << endl;

				u[i] = u[subMesh.getCellYSuccessor( i )] - subMesh.getHy();

			}

		}

		else if(x == subMesh.getDimensions().x() - 1 && (boundaryCondition & 2) && y == subMesh.getDimensions().y() - 1)

		{

			if((u[subMesh.getCellYPredecessor( i )] - u[i])/subMesh.getHy() > 1.0)

			{

				//cout << "x = n; y = n" << endl;

				u[i] = u[subMesh.getCellYPredecessor( i )] - subMesh.getHy();

			}

		}



      maxResidue = fu. absMax();

      //cout << fu. absMax() << endl;



      if( this -> cflCondition != 0.0 /*&& currentTau*maxResidue > this -> cflCondition*/)

    	  currentTau =  this -> cflCondition / maxResidue;

      //cout << "\b\b\b\b\b\b\b\b\b\b\b\b\r" << currentTau << flush;

      /*

      if( this -> cflCondition != 0.0 )

		else if(y == 0 && (boundaryCondition & 8) && x ==0)

		{

         maxResidue = fu. absMax();

         if( currentTau * maxResidue > this -> cflCondition )

			if((u[subMesh.getCellXSuccessor( i )] - u[i])/subMesh.getHx() > 1.0)

			{

        	 //cout << '\r' << "For " << currentTau * maxResidue << " : " << currentTau << " --> ";

        	 currentTau =  0.95 * this -> cflCondition / maxResidue;

        	 //cout << currentTau << " at " << this -> cflCondition   << flush;

            //currentTau *= 0.9;

//            continue;

				//cout << "y = 0; x = 0" << endl;

				u[i] = u[subMesh.getCellXSuccessor( i )] - subMesh.getHx();

			}

		}

		else if(y == 0 && (boundaryCondition & 8) && x == subMesh.getDimensions().x() - 1)

		{

			if((u[subMesh.getCellXPredecessor( i )] - u[i])/subMesh.getHx() > 1.0)

			{

				//cout << "y = 0; x = n" << endl;

				u[i] = u[subMesh.getCellXPredecessor( i )] - subMesh.getHx();

			}

		}

      }*/



      //RealType newResidue( 0.0 );

      //computeNewTimeLevel( u, currentTau, newResidue );



      //skip = false;

      if( time + currentTau > this -> stopTime )

    	  currentTau = this -> stopTime - time;

      for( int i = 0; i < fu.getSize(); i ++ )

		else if(y == subMesh.getDimensions().y() - 1 && (boundaryCondition & 1) && x ==0)

		{

    	  if((u[i]+currentTau * fu[ i ])*Sign(u[i]) < 0.0 && fu[i] != 0.0 && u[i] != 0.0 )

			if((u[subMesh.getCellXSuccessor( i )] - u[i])/subMesh.getHx() > 1.0)			{

				//cout << "y = n; x = 0" << endl;

				u[i] = u[subMesh.getCellXSuccessor( i )] - subMesh.getHx();

			}

		}

		else if(y == subMesh.getDimensions().y() - 1 && (boundaryCondition & 1) && x == subMesh.getDimensions().x() - 1)

		{

//    		  skip = true;

    		  currentTau = fabs(u[i]/(2.0*fu[i]));

    		  //currentTau *= 0.9;

//    		  i=fu.getSize();

			if((u[subMesh.getCellXPredecessor( i )] - u[i])/subMesh.getHx() > 1.0)

			{

				//cout << "y = n; x = n" << endl;

				u[i] = u[subMesh.getCellXPredecessor( i )] - subMesh.getHx();

			}

		}

	}*/



     // if(skip)

      //{

    //	  continue;

      //}

	/**/



      for( int i = 0; i < fu.getSize(); i ++ )



	bool tmp = false;

	for(int i = 0; i < u.getSize(); i++)

	{

    	  //double add = currentTau * fu[ i ];

   		  u[ i ] += currentTau * fu[ i ];

          //localResidue += fabs( add );

		if(u[0]*u[i] <= 0.0)

			tmp=true;

	}



      //this->setResidue( newResidue );



      /****

       * When time is close to stopTime the new residue

       * may be inaccurate significantly.

       */

      //if( currentTau + time == this -> stopTime ) this->setResidue( lastResidue );

      time += currentTau;

	if(tmp)

	{}

	else if(boundaryCondition == 4)

	{

		for(int i = 0; i < u.getSize(); i=subMesh.getCellYSuccessor(i))

		{

			for(int j = i; j < subMesh.getDimensions().x(); j=subMesh.getCellXSuccessor(j))

			{

				u[j] = u[i];

			}

		}

	}

	else if(boundaryCondition == 8)

	{

		for(int i = 0; i < subMesh.getDimensions().x(); i=subMesh.getCellXSuccessor(i))

		{

			for(int j = i; j < u.getSize(); j=subMesh.getCellYSuccessor(j))

			{

				u[j] = u[i];

			}

		}

	}

	else if(boundaryCondition == 2)

	{

		for(int i = subMesh.getDimensions().x() - 1; i < u.getSize(); i=subMesh.getCellYSuccessor(i))

		{

			for(int j = i; j > (i-1)*subMesh.getDimensions().x(); j=subMesh.getCellXPredecessor(j))

			{

				u[j] = u[i];

			}

		}

	}

	else if(boundaryCondition == 1)

	{

		for(int i = (subMesh.getDimensions().y() - 1)*subMesh.getDimensions().x(); i < u.getSize(); i=subMesh.getCellXSuccessor(i))

		{

			for(int j = i; j > 0; j=subMesh.getCellYPredecessor(j))

			{

				u[j] = u[i];

			}

		}

	}



      /*if( ! this->nextIteration() )

         return false;*/

	/**/



      /****

       * Compute the new time step.

       */

      /*if( this -> cflCondition != 0.0 )

        currentTau *= 1.1;*/

	/*if (u.max() > 0.0)

		this->stopTime *=(double) this->gridCols;*/



//      if( time + currentTau > this -> stopTime )

//      {

//         currentTau = this -> stopTime - time; //we don't want to keep such tau

//     }

      //else this -> tau = currentTau;



      //this -> refreshSolverMonitor();

   double time = 0.0;

   double currentTau = this->tau0;

   double finalTime = this->stopTime;// + 3.0*(u.max() - u.min());

   if( time + currentTau > finalTime ) currentTau = finalTime - time;



   double maxResidue( 0.0 );

   while( time < finalTime)

   {

      /****

       * Check stop conditions.

       * Compute the RHS

       */

      /*if( time >= this->getStopTime()  ||

          ( this -> getConvergenceResidue() != 0.0 && this->getResidue() < this -> getConvergenceResidue() ) )



      for( int i = 0; i < fu.getSize(); i ++ )

      {

         //this -> refreshSolverMonitor();

         return true;

      }*/

    	  fu[ i ] = scheme.getValue( this->subMesh, i, this->subMesh.getCellCoordinates(i), u, time, boundaryCondition );

      }

      maxResidue = fu. absMax();





      if( this -> cflCondition * maxResidue != 0.0)

    	  currentTau =  this -> cflCondition / maxResidue;





      if( time + currentTau > finalTime ) currentTau = finalTime - time;

      for( int i = 0; i < fu.getSize(); i ++ )

      {

    	  if((u[i]+currentTau * fu[ i ])*u[i] < 0.0 && fu[i] != 0.0 && u[i] != 0.0 )

    		  currentTau = fabs(u[i]/(2.0*fu[i]));

      }





      for( int i = 0; i < fu.getSize(); i ++ )

      {

    	  double add = u[i] + currentTau * fu[ i ];

    	  //if( fabs(u[i]) < fabs(add) or (this->subgridValues[subGridID] == this->currentStep +4) )

    		  u[ i ] = add;

      }

      time += currentTau;



      //cout << '\r' << flush;

     // cout << maxResidue << flush;

      //cout << maxResidue << "   " << currentTau << " @ " << time << flush;

   }



	/*if (u.max() > 0.0)

		this->stopTime /=(double) this->gridCols;*/



	VectorType solution;

	solution.setLike(u);
@@ -553,10 +622,6 @@ tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::runSubgrid( int boundary 
  	{

		solution[i]=u[i];

   	}



    cout << '\r' << flush;



    //cout << "RES: " << maxResidue  << " @ " << time << endl;

	return solution;

}

src/implementation/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h

+16 −3
Original line number Diff line number Diff line
@@ -130,18 +130,28 @@ Real parallelGodunovEikonalScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Re 
          	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	 const IndexType boundaryCondition ) const

{



	if ( (coordinates.x() == 0 && (boundaryCondition & 4)) or

	if ( ((coordinates.x() == 0 && (boundaryCondition & 4)) or

		 (coordinates.x() == mesh.getDimensions().x() - 1 && (boundaryCondition & 2)) or

		 (coordinates.y() == 0 && (boundaryCondition & 8)) or

		 (coordinates.y() == mesh.getDimensions().y() - 1  && (boundaryCondition & 1)))

		 /*and

		 !(		 (coordinates.y() == 0 or coordinates.y() == mesh.getDimensions().y() - 1)

				 and

				 ( coordinates.x() == 0 or coordinates.x() == mesh.getDimensions().x() - 1)

		  )*/

		)

	{

		return 0.0;

	}



	RealType acc = hx*hy*hx*hy;



	RealType nabla, xb, xf, yb, yf, signui;



	signui = sign(u[cellIndex],epsilon);



	if(fabs(u[cellIndex]) < acc) return 0.0;



	   if(signui > 0.0)

	   {

	/**/ /*  if(boundaryCondition & 2)
@@ -191,7 +201,8 @@ Real parallelGodunovEikonalScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Re 
			   yb = 0.0;



		   nabla = sqrt (xf*xf + xb*xb + yf*yf + yb*yb );



		   if(fabs(1.0-nabla) < acc)

			   return 0.0;

		   return signui*(1.0 - nabla);

	   }

	   else if (signui < 0.0)
@@ -246,6 +257,8 @@ Real parallelGodunovEikonalScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Re 


		   nabla = sqrt (xf*xf + xb*xb + yf*yf + yb*yb );



		   if(fabs(1.0-nabla) < acc)

			   return 0.0;

		   return signui*(1.0 - nabla);

	   }

	   else