From a8d4fd6de30a375de631f8382e28fae2d6e6576f Mon Sep 17 00:00:00 2001
From: Tomas Sobotik <sobotto4@fjfi.cvut.cz>
Date: Tue, 14 Apr 2015 18:20:21 +0200
Subject: [PATCH] Minor improvements
---
.../tnlParallelEikonalSolver_impl.h | 61 +++++----
.../parallelGodunovEikonal2D_impl.h | 128 ++++++------------
2 files changed, 78 insertions(+), 111 deletions(-)
diff --git a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h
index ec84e8c0a6..11272d8ecb 100644
--- a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h
+++ b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h
@@ -1151,20 +1151,9 @@ void tnlParallelEikonalSolver<SchemeHost, SchemeDevice, Device, double, int>::ru
}*/
//__syncthreads();
//end reduction
- __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]) ;
- __syncthreads();
- }
- if(l==0) maxResidue=sharedRes[l];
-
- sharedTau[l]=currentTau;
- __syncthreads();
- // if(l == 0)
- if( this -> cflCondition * maxResidue != 0.0)
- sharedTau[l] = this -> cflCondition / maxResidue;
+ 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())
@@ -1172,14 +1161,44 @@ void tnlParallelEikonalSolver<SchemeHost, SchemeDevice, Device, double, int>::ru
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];
+ currentTau=sharedTau[blockDim.x*blockDim.y - 1];
+ /*if( this -> cflCondition * maxResidue != 0.0)
+ currentTau = Min(this -> cflCondition / maxResidue, currentTau);*/
+ }
+ __syncthreads();
+
+
//__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();
+
//
diff --git a/src/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h b/src/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h
index c37134862f..31d6ca7982 100644
--- a/src/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h
+++ b/src/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h
@@ -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(coordinates.x() == mesh.getDimensions().x() - 1)
+ xf = ((u[mesh.template getCellNextToCell<-1,0>( cellIndex )] - u[cellIndex])/hx);
+ else
+ xf = ((u[mesh.template getCellNextToCell<1,0>( cellIndex )] - u[cellIndex])/hx);
+
+ if(coordinates.x() == 0)
+ xb = ((u[cellIndex] - u[mesh.template getCellNextToCell<1,0>( cellIndex )])/hx);
+ else
+ xb = ((u[cellIndex] - u[mesh.template getCellNextToCell<-1,0>( cellIndex )])/hx);
+
+ if(coordinates.y() == mesh.getDimensions().y() - 1)
+ yf = ((u[mesh.template getCellNextToCell<0,-1>( cellIndex )] - u[cellIndex])/hy);
+ else
+ yf = ((u[mesh.template getCellNextToCell<0,1>( cellIndex )] - u[cellIndex])/hy);
+
+ if(coordinates.y() == 0)
+ yb = ((u[cellIndex] - u[mesh.template getCellNextToCell<0,1>( cellIndex )])/hy);
+ else
+ yb = ((u[cellIndex] - u[mesh.template getCellNextToCell<0,-1>( cellIndex )])/hy);
+
+
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);
- else
- xf = negativePart((u[mesh.template getCellNextToCell<1,0>( cellIndex )] - u[cellIndex])/hx);
+ xf = negativePart(xf);
- /**/ /* 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);
- else
- xb = positivePart((u[cellIndex] - u[mesh.template getCellNextToCell<-1,0>( cellIndex )])/hx);
+ xb = positivePart(xb);
- /**/ /* 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);
- else
- yf = negativePart((u[mesh.template getCellNextToCell<0,1>( cellIndex )] - u[cellIndex])/hy);
+ yf = negativePart(yf);
- /**/ /* 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);
- else
- yb = positivePart((u[cellIndex] - u[mesh.template getCellNextToCell<0,-1>( cellIndex )])/hy);
+ 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,17 +395,16 @@ 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;
- return signui*(1.0 - nabla);
- }
- else
- {
- return 0.0;
}
+
+
+ nabla = sqrt (xf*xf + xb*xb + yf*yf + yb*yb );
+ return signui*(1.0 - nabla);
+
+
+
+
}
--
GitLab