Commit a8d4fd6d authored by Tomas Sobotik's avatar Tomas Sobotik
Browse files

Minor improvements

parent 28491294

examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h

+36 −25
Original line number Diff line number Diff line
@@ -1151,35 +1151,54 @@ void tnlParallelEikonalSolver<SchemeHost, SchemeDevice, Device, double, int>::ru 
  	}*/

  	//__syncthreads();

  																													//end reduction

      sharedTau[l]=this -> cflCondition/sharedRes[l];

      if((u[l]+sharedTau[l] * fu)*u[l] < 0.0)

    	  sharedTau[l] = fabs(u[l]/(2.0*fu));



      if(l == 0)

    	  if(sharedTau[l] > 1.0 * this->subMesh.getHx())

    		  sharedTau[l] = 1.0 * this->subMesh.getHx();



      if(l == 1)

    	  if( time + sharedTau[l] > finalTime ) sharedTau[l] = finalTime - time;

      __syncthreads();



   //   __syncthreads();

  	/*for(unsigned int s = blockDim.x*blockDim.y/2; s>0; s>>=1)

  	{

  		if( l < s )

  			sharedTau[l] = Min(sharedTau[l],sharedTau[l+s]) ;

  		__syncthreads();

  	}

  	if(l==0) currentTau=sharedTau[l];

  	__syncthreads();*/





	for(unsigned int s = blockDim.x*blockDim.y/2; s>0; s>>=1)

	{

		if( l < s )

			sharedRes[l] = Max(sharedRes[l],sharedRes[l+s]);

		if(l >= blockDim.x*blockDim.y - s)

			sharedTau[l] = Min(sharedTau[l],sharedTau[l-s]);

		__syncthreads();

	}

	if(l==0) maxResidue=sharedRes[l];



    sharedTau[l]=currentTau;

	if(l==0)

	{

		maxResidue=sharedRes[l];

		currentTau=sharedTau[blockDim.x*blockDim.y - 1];

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

			currentTau = Min(this -> cflCondition / maxResidue, currentTau);*/

	}

	__syncthreads();

     // if(l == 0)

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

    		  sharedTau[l] =  this -> cflCondition / maxResidue;



      if(l == 0)

    	  if(sharedTau[l] > 1.0 * this->subMesh.getHx())

    		  sharedTau[l] = 1.0 * this->subMesh.getHx();



      if(l == 1)

    	  if( time + sharedTau[l] > finalTime ) sharedTau[l] = finalTime - time;



      //__syncthreads();

 //



      //double tau2 = finalTime;



      if((u[l]+sharedTau[l] * fu)*u[l] < 0.0 && fu != 0.0 && u[l] != 0.0 )

    	  sharedTau[l] = fabs(u[l]/(2.0*fu));





      	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  //start reduction
@@ -1196,15 +1215,7 @@ void tnlParallelEikonalSolver<SchemeHost, SchemeDevice, Device, double, int>::ru 






    __syncthreads();

	for(unsigned int s = blockDim.x*blockDim.y/2; s>0; s>>=1)

	{

		if( l < s )

			sharedTau[l] = Min(sharedTau[l],sharedTau[l+s]) ;

		__syncthreads();

	}

	if(l==0) currentTau=sharedTau[l];

	__syncthreads();







//

src/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h

+42 −86
Original line number Diff line number Diff line
@@ -315,63 +315,53 @@ Real parallelGodunovEikonalScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Re 
		 (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;

	//RealType acc = hx*hy*hx*hy;



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

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

#ifdef HAVE_CUDA

	//printf("%d   :    %d ;;;; %d   :   %d  , %f \n",threadIdx.x, mesh.getDimensions().x() , threadIdx.y,mesh.getDimensions().y(), epsilon );

#endif

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



	   if(signui > 0.0)

	   {

	/**/ /*  if(boundaryCondition & 2)

			   xf = (u[mesh.getCellXSuccessor( cellIndex )] - u[cellIndex])/hx;

		   else *//*if(boundaryCondition & 4)

			   xf = 0.0;

		   else /**/if(coordinates.x() == mesh.getDimensions().x() - 1)

			   xf = negativePart((u[mesh.template getCellNextToCell<-1,0>( cellIndex )] - u[cellIndex])/hx);





	   if(coordinates.x() == mesh.getDimensions().x() - 1)

		   xf = ((u[mesh.template getCellNextToCell<-1,0>( cellIndex )] - u[cellIndex])/hx);

	   else

			   xf = negativePart((u[mesh.template getCellNextToCell<1,0>( cellIndex )] - u[cellIndex])/hx);

		   xf = ((u[mesh.template getCellNextToCell<1,0>( cellIndex )] - u[cellIndex])/hx);



	/**/ /*  if(boundaryCondition & 4)

			   xb = (u[cellIndex] - u[mesh.getCellXPredecessor( cellIndex )])/hx;

		   else *//*if(boundaryCondition & 2)

			   xb = 0.0;

		   else /**/if(coordinates.x() == 0)

			   xb = positivePart((u[cellIndex] - u[mesh.template getCellNextToCell<1,0>( cellIndex )])/hx);

	   if(coordinates.x() == 0)

		   xb = ((u[cellIndex] - u[mesh.template getCellNextToCell<1,0>( cellIndex )])/hx);

	   else

			   xb = positivePart((u[cellIndex] - u[mesh.template getCellNextToCell<-1,0>( cellIndex )])/hx);

		   xb = ((u[cellIndex] - u[mesh.template getCellNextToCell<-1,0>( cellIndex )])/hx);



	/**/  /* if(boundaryCondition & 1)

			   yf = (u[mesh.getCellYSuccessor( cellIndex )] - u[cellIndex])/hy;

		   else *//*if(boundaryCondition & 8)

			   yf = 0.0;

		   else /**/if(coordinates.y() == mesh.getDimensions().y() - 1)

			   yf = negativePart((u[mesh.template getCellNextToCell<0,-1>( cellIndex )] - u[cellIndex])/hy);

	   if(coordinates.y() == mesh.getDimensions().y() - 1)

		   yf = ((u[mesh.template getCellNextToCell<0,-1>( cellIndex )] - u[cellIndex])/hy);

	   else

			   yf = negativePart((u[mesh.template getCellNextToCell<0,1>( cellIndex )] - u[cellIndex])/hy);

		   yf = ((u[mesh.template getCellNextToCell<0,1>( cellIndex )] - u[cellIndex])/hy);



	/**/  /* if(boundaryCondition & 8)

			   yb = (u[cellIndex] - u[mesh.getCellYPredecessor( cellIndex )])/hy;

		   else *//*if(boundaryCondition & 1)

			   yb = 0.0;

		   else /**/if(coordinates.y() == 0)

			   yb = positivePart((u[cellIndex] - u[mesh.template getCellNextToCell<0,1>( cellIndex )])/hy);

	   if(coordinates.y() == 0)

		   yb = ((u[cellIndex] - u[mesh.template getCellNextToCell<0,1>( cellIndex )])/hy);

	   else

			   yb = positivePart((u[cellIndex] - u[mesh.template getCellNextToCell<0,-1>( cellIndex )])/hy);

		   yb = ((u[cellIndex] - u[mesh.template getCellNextToCell<0,-1>( cellIndex )])/hy);







	   if(signui > 0.0)

	   {

		   xf = negativePart(xf);



		   xb = positivePart(xb);



		   yf = negativePart(yf);



		   yb = positivePart(yb);



		   if(xb + xf > 0.0)

			   xf = 0.0;
@@ -383,50 +373,17 @@ Real parallelGodunovEikonalScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Re 
		   else

			   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)

	   {



	/**/  /* if(boundaryCondition & 2)

			   xf = (u[mesh.getCellXSuccessor( cellIndex )] - u[cellIndex])/hx;

		   else*//* if(boundaryCondition & 4)

			   xf = 0.0;

		   else /**/if(coordinates.x() == mesh.getDimensions().x() - 1)

			   xf = positivePart((u[mesh.template getCellNextToCell<-1,0>( cellIndex )] - u[cellIndex])/hx);

		   else

			   xf = positivePart((u[mesh.template getCellNextToCell<1,0>( cellIndex )] - u[cellIndex])/hx);

		   xb = negativePart(xb);



	/**/  /* if(boundaryCondition & 4)

			   xb = (u[cellIndex] - u[mesh.getCellXPredecessor( cellIndex )])/hx;

		   else*//* if(boundaryCondition & 2)

			   xb = 0.0;

		   else /**/if(coordinates.x() == 0)

			   xb = negativePart((u[cellIndex] - u[mesh.template getCellNextToCell<1,0>( cellIndex )])/hx);

		   else

			   xb = negativePart((u[cellIndex] - u[mesh.template getCellNextToCell<-1,0>( cellIndex )])/hx);

		   xf = positivePart(xf);



	/**/ /*  if(boundaryCondition & 1)

			   yf = (u[mesh.getCellYSuccessor( cellIndex )] - u[cellIndex])/hy;

		   else *//*if(boundaryCondition & 8)

			   yf = 0.0;

		   else /**/if(coordinates.y() == mesh.getDimensions().y() - 1)

			   yf = positivePart((u[mesh.template getCellNextToCell<0,-1>( cellIndex )] - u[cellIndex])/hy);

		   else

			   yf = positivePart((u[mesh.template getCellNextToCell<0,1>( cellIndex )] - u[cellIndex])/hy);



	/**/  /* if(boundaryCondition & 8)

			   yb = (u[cellIndex] - u[mesh.getCellYPredecessor( cellIndex )])/hy;

		   else*//* if(boundaryCondition & 1)

			   yb = 0.0;

		   else /**/if(coordinates.y() == 0)

			   yb = negativePart((u[cellIndex] - u[mesh.template getCellNextToCell<0,1>( cellIndex )])/hy);

		   else

			   yb = negativePart((u[cellIndex] - u[mesh.template getCellNextToCell<0,-1>( cellIndex )])/hy);

		   yb = negativePart(yb);



		   yf = positivePart(yf);



		   if(xb + xf > 0.0)

			   xb = 0.0;
@@ -438,16 +395,15 @@ Real parallelGodunovEikonalScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Re 
		   else

			   yf = 0.0;



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

	   }



		  // if(fabs(1.0-nabla) < acc)

		//	   return 0.0;





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

	   return signui*(1.0 - nabla);

	   }

	   else

	   {

		   return 0.0;

	   }









}