Commit 6fee7bf9 authored by Tomas Sobotik's avatar Tomas Sobotik
Browse files

Added constraint for Tau & few other changes

parent 54b44d9d

examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h

+55 −10
Original line number Diff line number Diff line
@@ -59,7 +59,7 @@ bool tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::init( const tnlPara 
	//cout << this->mesh.getCellCenter(0) << endl;



	this->delta = parameters.GetParameter <double>("delta");

	this->delta *= sqrt(this->mesh.getHx()*this->mesh.getHx() + this->mesh.getHy()*this->mesh.getHy());

	this->delta *= this->mesh.getHx()*this->mesh.getHy();



	cout << "Setting delta to " << this->delta << endl;
@@ -75,6 +75,8 @@ bool tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::init( const tnlPara 
	this->stopTime /= (double)(this->gridCols);

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

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

	this->stopTime = ((double)(this->n) - 1.0)*parameters.GetParameter <double>("stop-time")*this->mesh.getHx();

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



	cout << "Initializating scheme..." << endl;

	if(!this->scheme.init(parameters))
@@ -134,7 +136,7 @@ void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::run() 
			end=true;

		else

			end=false;



#pragma omp parallel for num_threads(3) schedule(dynamic)

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

		{

			if(getSubgridValue(i) != INT_MAX)
@@ -162,28 +164,32 @@ void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::run() 
					this->calculationsCount[i]++;

				}



/*

				if( (getBoundaryCondition(i) & 1) || (getBoundaryCondition(i) & 2) )



				if( ((getBoundaryCondition(i) & 2) )//|| (getBoundaryCondition(i) & 1))

					/*	&&(!(getBoundaryCondition(i) & 5) && !(getBoundaryCondition(i) & 10)) */)

				{

					//cout << "3 @ " << getBoundaryCondition(i) << endl;

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

				}

				if( (getBoundaryCondition(i) & 1) || (getBoundaryCondition(i) & 4) )

				if( ((getBoundaryCondition(i) & 4) )//|| (getBoundaryCondition(i) & 1))

					/*	&&(!(getBoundaryCondition(i) & 3) && !(getBoundaryCondition(i) & 12)) */)

				{

					//cout << "5 @ " << getBoundaryCondition(i) << endl;

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

				}

				if( (getBoundaryCondition(i) & 8) || (getBoundaryCondition(i) & 2) )

				if( ((getBoundaryCondition(i) & 2) )//|| (getBoundaryCondition(i) & 8))

					/*	&&(!(getBoundaryCondition(i) & 12) && !(getBoundaryCondition(i) & 3))*/ )

				{

					//cout << "10 @ " << getBoundaryCondition(i) << endl;

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

				}

				if( (getBoundaryCondition(i) & 8) || (getBoundaryCondition(i) & 4) )

				if(   ((getBoundaryCondition(i) & 4) )//|| (getBoundaryCondition(i) & 8))

					/*&&(!(getBoundaryCondition(i) & 10) && !(getBoundaryCondition(i) & 5)) */)

				{

					//cout << "12 @ " << getBoundaryCondition(i) << endl;

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

				}

*/





				/*if(getBoundaryCondition(i))

				{
@@ -683,7 +689,38 @@ tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::runSubgrid( int boundary 
				u[i*this->n + j] = value;// u[j];

	}



/*



	else if(boundaryCondition == 5)

	{

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

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

				//if(fabs(u[i*this->n + j]) <  fabs(u[i*this->n]))

				u[i*this->n + j] = value;// u[i*this->n];

	}

	else if(boundaryCondition == 10)

	{

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

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

				//if(fabs(u[i*this->n + j]) < fabs(u[(i+1)*this->n - 1]))

				u[i*this->n + j] = value;// u[(i+1)*this->n - 1];

	}

	else if(boundaryCondition == 3)

	{

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

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

				//if(fabs(u[i*this->n + j]) < fabs(u[j + this->n*(this->n - 1)]))

				u[i*this->n + j] = value;// u[j + this->n*(this->n - 1)];

	}

	else if(boundaryCondition == 12)

	{

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

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

				//if(fabs(u[i*this->n + j]) < fabs(u[j]))

				u[i*this->n + j] = value;// u[j];

	}



*/



	/**/
@@ -713,6 +750,14 @@ tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::runSubgrid( int boundary 


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

    	  currentTau =  this -> cflCondition / maxResidue;



     /* if (maxResidue < 0.05)

    	  cout << "Max < 0.05" << endl;*/

      if(currentTau > 1.0 * this->subMesh.getHx())

      {

    	  //cout << currentTau << " >= " << 2.0 * this->subMesh.getHx() << endl;

    	  currentTau = 1.0 * this->subMesh.getHx();

      }

      /*if(maxResidue > lastResidue)

    	  currentTau *=(1.0/10.0);*/
@@ -720,7 +765,7 @@ tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::runSubgrid( int boundary 
      if( time + currentTau > finalTime ) currentTau = finalTime - time;

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

      {

    	  cout << "Too big RHS! i = " << i << ", fu = " << fu[i] << ", u = " << u[i] << endl;

    	  //cout << "Too big RHS! i = " << i << ", fu = " << fu[i] << ", u = " << u[i] << endl;

    	  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]));
@@ -739,7 +784,7 @@ tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::runSubgrid( int boundary 
     //cout << maxResidue << "   " << currentTau << " @ " << time << flush;

     //lastResidue = maxResidue;

   }

   cout << "Time: " << time << ", Res: " << maxResidue <<endl;

   //cout << "Time: " << time << ", Res: " << maxResidue <<endl;

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

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