From c99c65a1cee4c8b0019b4ae7d3ee662d8c3dd572 Mon Sep 17 00:00:00 2001
From: Tomas Sobotik <sobotto4@fjfi.cvut.cz>
Date: Tue, 29 Mar 2016 22:53:21 +0200
Subject: [PATCH] CPU versions updated for parallel solvers
---
.../tnlParallelMapSolver.h | 2 +-
.../tnlParallelMapSolver2D_impl.h | 412 ++++++----------
.../tnlParallelEikonalSolver2D_impl.h | 44 +-
.../tnlParallelEikonalSolver3D_impl.h | 465 ++++++++++++------
.../parallelGodunovEikonal3D_impl.h | 14 +-
.../godunov-eikonal/parallelGodunovMap.h | 3 +-
.../parallelGodunovMap2D_impl.h | 171 ++-----
7 files changed, 550 insertions(+), 561 deletions(-)
diff --git a/examples/hamilton-jacobi-parallel-map/tnlParallelMapSolver.h b/examples/hamilton-jacobi-parallel-map/tnlParallelMapSolver.h
index 5c3c5ca2fd..fc61312308 100644
--- a/examples/hamilton-jacobi-parallel-map/tnlParallelMapSolver.h
+++ b/examples/hamilton-jacobi-parallel-map/tnlParallelMapSolver.h
@@ -90,7 +90,7 @@ public:
void insertSubgrid( VectorType u, const int i );
- VectorType runSubgrid( int boundaryCondition, VectorType u, int subGridID);
+ VectorType runSubgrid( int boundaryCondition, VectorType u, int subGridID,VectorType map);
tnlMeshFunction<MeshType> u0;
diff --git a/examples/hamilton-jacobi-parallel-map/tnlParallelMapSolver2D_impl.h b/examples/hamilton-jacobi-parallel-map/tnlParallelMapSolver2D_impl.h
index cc611b9501..6c0e48681a 100644
--- a/examples/hamilton-jacobi-parallel-map/tnlParallelMapSolver2D_impl.h
+++ b/examples/hamilton-jacobi-parallel-map/tnlParallelMapSolver2D_impl.h
@@ -32,7 +32,7 @@ template< typename SchemeHost, typename SchemeDevice, typename Device>
tnlParallelMapSolver<2,SchemeHost, SchemeDevice, Device, double, int>::tnlParallelMapSolver()
{
cout << "a" << endl;
- this->device = tnlCudaDevice; /////////////// tnlCuda Device --- vypocet na GPU, tnlHostDevice --- vypocet na CPU
+ this->device = tnlHostDevice; /////////////// tnlCuda Device --- vypocet na GPU, tnlHostDevice --- vypocet na CPU
#ifdef HAVE_CUDA
if(this->device == tnlCudaDevice)
@@ -162,33 +162,43 @@ bool tnlParallelMapSolver<2,SchemeHost, SchemeDevice, Device, double, int>::init
if(this->device == tnlHostDevice)
{
- for(int i = 0; i < this->subgridValues.getSize(); i++)
- {
+ VectorType tmp_map;
+ tmp_map.setSize(this->n * this->n);
+ for(int i = 0; i < this->subgridValues.getSize(); i++)
+ {
- if(! tmp[i].setSize(this->n * this->n))
- cout << "Could not allocate tmp["<< i <<"] array." << endl;
- tmp[i] = getSubgrid(i);
- containsCurve = false;
+ if(! tmp[i].setSize(this->n * this->n))
+ cout << "Could not allocate tmp["<< i <<"] array." << endl;
+ tmp[i] = getSubgrid(i);
+ containsCurve = false;
- for(int j = 0; j < tmp[i].getSize(); j++)
- {
- if(tmp[i][0]*tmp[i][j] <= 0.0)
+ for(int j = 0; j < tmp[i].getSize(); j++)
+ {
+ if(tmp[i][0]*tmp[i][j] <= 0.0)
+ {
+ containsCurve = true;
+ j=tmp[i].getSize();
+ }
+
+ }
+ if(containsCurve)
{
- containsCurve = true;
- j=tmp[i].getSize();
+ for( int j = 0; j < tmp_map.getSize(); j++)
+ {
+ tmp_map[j] = this->map_stretched[ (i / this->gridCols) * this->n*this->n*this->gridCols
+ + (i % this->gridCols) * this->n
+ + (j/this->n) * this->n*this->gridCols
+ + (j % this->n) ];
+ }
+ //cout << "Computing initial SDF on subgrid " << i << "." << endl;
+ tmp[i] = runSubgrid(0, tmp[i],i,tmp_map);
+ insertSubgrid(tmp[i], i);
+ setSubgridValue(i, 4);
+ //cout << "Computed initial SDF on subgrid " << i << "." << endl;
}
+ containsCurve = false;
}
- if(containsCurve)
- {
- //cout << "Computing initial SDF on subgrid " << i << "." << endl;
- insertSubgrid(runSubgrid(0, tmp[i],i), i);
- setSubgridValue(i, 4);
- //cout << "Computed initial SDF on subgrid " << i << "." << endl;
- }
- containsCurve = false;
-
- }
}
#ifdef HAVE_CUDA
else if(this->device == tnlCudaDevice)
@@ -249,13 +259,13 @@ void tnlParallelMapSolver<2,SchemeHost, SchemeDevice, Device, double, int>::run(
if(this->device == tnlHostDevice)
{
- bool end = false;
- while ((this->boundaryConditions.max() > 0 ) || !end)
+// bool end = false;
+ while ((this->boundaryConditions.max() > 0 )/* || !end*/)
{
- if(this->boundaryConditions.max() == 0 )
- end=true;
- else
- end=false;
+// if(this->boundaryConditions.max() == 0 )
+// end=true;
+// else
+// end=false;
#ifdef HAVE_OPENMP
#pragma omp parallel for num_threads(4) schedule(dynamic)
#endif
@@ -263,56 +273,110 @@ void tnlParallelMapSolver<2,SchemeHost, SchemeDevice, Device, double, int>::run(
{
if(getSubgridValue(i) != INT_MAX)
{
+ VectorType tmp, tmp_map;
+ tmp.setSize(this->n * this->n);
+ tmp_map.setSize(this->n * this->n);
+ for( int j = 0; j < tmp_map.getSize(); j++)
+ {
+ tmp_map[j] = this->map_stretched[ (i / this->gridCols) * this->n*this->n*this->gridCols
+ + (i % this->gridCols) * this->n
+ + (j/this->n) * this->n*this->gridCols
+ + (j % this->n) ];
+ }
//cout << "subMesh: " << i << ", BC: " << getBoundaryCondition(i) << endl;
if(getSubgridValue(i) == currentStep+4)
{
- if(getBoundaryCondition(i) & 1)
- {
- insertSubgrid( runSubgrid(1, getSubgrid(i),i), i);
- this->calculationsCount[i]++;
- }
- if(getBoundaryCondition(i) & 2)
- {
- insertSubgrid( runSubgrid(2, getSubgrid(i),i), i);
- this->calculationsCount[i]++;
- }
- if(getBoundaryCondition(i) & 4)
- {
- insertSubgrid( runSubgrid(4, getSubgrid(i),i), i);
- this->calculationsCount[i]++;
+ if(getBoundaryCondition(i) & 1)
+ {
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(1, tmp ,i,tmp_map);
+ insertSubgrid( tmp, i);
+ this->calculationsCount[i]++;
+ }
+ if(getBoundaryCondition(i) & 2)
+ {
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(2, tmp ,i,tmp_map);
+ insertSubgrid( tmp, i);
+ this->calculationsCount[i]++;
+ }
+ if(getBoundaryCondition(i) & 4)
+ {
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(4, tmp ,i,tmp_map);
+ insertSubgrid( tmp, i);
+ this->calculationsCount[i]++;
+ }
+ if(getBoundaryCondition(i) & 8)
+ {
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(8, tmp ,i,tmp_map);
+ insertSubgrid( tmp, i);
+ this->calculationsCount[i]++;
+ }
}
- if(getBoundaryCondition(i) & 8)
+ else
{
- insertSubgrid( runSubgrid(8, getSubgrid(i),i), i);
- this->calculationsCount[i]++;
- }
+
+ if(getBoundaryCondition(i) == 1)
+ {
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(1, tmp ,i,tmp_map);
+ insertSubgrid( tmp, i);
+ this->calculationsCount[i]++;
+ }
+ if(getBoundaryCondition(i) == 2)
+ {
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(2, tmp ,i,tmp_map);
+ insertSubgrid( tmp, i);
+ this->calculationsCount[i]++;
+ }
+ if(getBoundaryCondition(i) == 4)
+ {
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(4, tmp ,i,tmp_map);
+ insertSubgrid( tmp, i);
+ this->calculationsCount[i]++;
+ }
+ if(getBoundaryCondition(i) == 8)
+ {
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(8, tmp ,i,tmp_map);
+ insertSubgrid( tmp, i);
+ this->calculationsCount[i]++;
+ }
}
- if( ((getBoundaryCondition(i) & 2) )|| (getBoundaryCondition(i) & 1)//)
- /* &&(!(getBoundaryCondition(i) & 5) && !(getBoundaryCondition(i) & 10)) */)
+ if(getBoundaryCondition(i) & 3)
{
//cout << "3 @ " << getBoundaryCondition(i) << endl;
- insertSubgrid( runSubgrid(3, getSubgrid(i),i), i);
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(3, tmp ,i,tmp_map);
+ insertSubgrid( tmp, i);
}
- if( ((getBoundaryCondition(i) & 4) )|| (getBoundaryCondition(i) & 1)//)
- /* &&(!(getBoundaryCondition(i) & 3) && !(getBoundaryCondition(i) & 12)) */)
+ if(getBoundaryCondition(i) & 5)
{
//cout << "5 @ " << getBoundaryCondition(i) << endl;
- insertSubgrid( runSubgrid(5, getSubgrid(i),i), i);
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(5, tmp ,i,tmp_map);
+ insertSubgrid( tmp, i);
}
- if( ((getBoundaryCondition(i) & 2) )|| (getBoundaryCondition(i) & 8)//)
- /* &&(!(getBoundaryCondition(i) & 12) && !(getBoundaryCondition(i) & 3))*/ )
+ if(getBoundaryCondition(i) & 10)
{
//cout << "10 @ " << getBoundaryCondition(i) << endl;
- insertSubgrid( runSubgrid(10, getSubgrid(i),i), i);
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(10, tmp ,i,tmp_map);
+ insertSubgrid( tmp, i);
}
- if( ((getBoundaryCondition(i) & 4) )|| (getBoundaryCondition(i) & 8)//)
- /*&&(!(getBoundaryCondition(i) & 10) && !(getBoundaryCondition(i) & 5)) */)
+ if(getBoundaryCondition(i) & 12)
{
//cout << "12 @ " << getBoundaryCondition(i) << endl;
- insertSubgrid( runSubgrid(12, getSubgrid(i),i), i);
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(12, tmp ,i,tmp_map);
+ insertSubgrid( tmp, i);
}
@@ -428,8 +492,8 @@ void tnlParallelMapSolver<2,SchemeHost, SchemeDevice, Device, double, int>::sync
int tmp1, tmp2;
int grid1, grid2;
- if(this->currentStep & 1)
- {
+// if(this->currentStep & 1)
+// {
for(int j = 0; j < this->gridRows - 1; j++)
{
for (int i = 0; i < this->gridCols*this->n; i++)
@@ -465,9 +529,9 @@ void tnlParallelMapSolver<2,SchemeHost, SchemeDevice, Device, double, int>::sync
}
}
- }
- else
- {
+// }
+// else
+// {
for(int i = 1; i < this->gridCols; i++)
{
for (int j = 0; j < this->gridRows*this->n; j++)
@@ -502,7 +566,7 @@ void tnlParallelMapSolver<2,SchemeHost, SchemeDevice, Device, double, int>::sync
}
}
}
- }
+// }
this->currentStep++;
@@ -619,8 +683,10 @@ void tnlParallelMapSolver<2,SchemeHost, SchemeDevice, Device, double, int>::stre
if(l>0)
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];
+ if(abs(this->u0[i-k]) < mesh.template getSpaceStepsProducts< 1, 0 >()+mesh.template getSpaceStepsProducts< 0, 1 >() )
+ this->work_u[i] = this->u0[i-k];
+ else
+ this->work_u[i] = Sign(this->u0[i-k])*MAP_SOLVER_MAX_VALUE;
///**/cout << "bla " << i << " " << this->map_stretched.getSize() << " " << i-k << " " << this->map.getData().getSize() << endl;
this->map_stretched[i] = this->map[i-k];
///**/cout << "bla " << i << " " << this->map_stretched.getSize() << " " << i-k << " " << this->map.getData().getSize() << endl;
@@ -692,7 +758,7 @@ void tnlParallelMapSolver<2,SchemeHost, SchemeDevice, Device, double, int>::inse
template< typename SchemeHost, typename SchemeDevice, typename Device>
typename tnlParallelMapSolver<2,SchemeHost, SchemeDevice, Device, double, int>::VectorType
-tnlParallelMapSolver<2,SchemeHost, SchemeDevice, Device, double, int>::runSubgrid( int boundaryCondition, VectorType u, int subGridID)
+tnlParallelMapSolver<2,SchemeHost, SchemeDevice, Device, double, int>::runSubgrid( int boundaryCondition, VectorType u, int subGridID,VectorType map)
{
VectorType fu;
@@ -700,181 +766,6 @@ tnlParallelMapSolver<2,SchemeHost, SchemeDevice, Device, double, int>::runSubgri
fu.setLike(u);
fu.setValue( 0.0 );
-/*
- * Insert Euler-Solver Here
- */
-
- /**/
-
- /*for(int i = 0; i < u.getSize(); i++)
- {
- int x = this->subMesh.getCellCoordinates(i).x();
- int y = this->subMesh.getCellCoordinates(i).y();
-
- if(x == 0 && (boundaryCondition & 4) && y ==0)
- {
- if((u[subMesh.getCellYSuccessor( i )] - u[i])/subMesh.template getSpaceStepsProducts< 0, 1 >() > 1.0)
- {
- //cout << "x = 0; y = 0" << endl;
- u[i] = u[subMesh.getCellYSuccessor( i )] - subMesh.template getSpaceStepsProducts< 0, 1 >();
- }
- }
- else if(x == 0 && (boundaryCondition & 4) && y == subMesh.getDimensions().y() - 1)
- {
- if((u[subMesh.getCellYPredecessor( i )] - u[i])/subMesh.template getSpaceStepsProducts< 0, 1 >() > 1.0)
- {
- //cout << "x = 0; y = n" << endl;
- u[i] = u[subMesh.getCellYPredecessor( i )] - subMesh.template getSpaceStepsProducts< 0, 1 >();
- }
- }
-
-
- else if(x == subMesh.getDimensions().x() - 1 && (boundaryCondition & 2) && y ==0)
- {
- if((u[subMesh.getCellYSuccessor( i )] - u[i])/subMesh.template getSpaceStepsProducts< 0, 1 >() > 1.0)
- {
- //cout << "x = n; y = 0" << endl;
- u[i] = u[subMesh.getCellYSuccessor( i )] - subMesh.template getSpaceStepsProducts< 0, 1 >();
- }
- }
- else if(x == subMesh.getDimensions().x() - 1 && (boundaryCondition & 2) && y == subMesh.getDimensions().y() - 1)
- {
- if((u[subMesh.getCellYPredecessor( i )] - u[i])/subMesh.template getSpaceStepsProducts< 0, 1 >() > 1.0)
- {
- //cout << "x = n; y = n" << endl;
- u[i] = u[subMesh.getCellYPredecessor( i )] - subMesh.template getSpaceStepsProducts< 0, 1 >();
- }
- }
-
-
- else if(y == 0 && (boundaryCondition & 8) && x ==0)
- {
- if((u[subMesh.getCellXSuccessor( i )] - u[i])/subMesh.template getSpaceStepsProducts< 1, 0 >() > 1.0)
- {
- //cout << "y = 0; x = 0" << endl;
- u[i] = u[subMesh.getCellXSuccessor( i )] - subMesh.template getSpaceStepsProducts< 1, 0 >();
- }
- }
- else if(y == 0 && (boundaryCondition & 8) && x == subMesh.getDimensions().x() - 1)
- {
- if((u[subMesh.getCellXPredecessor( i )] - u[i])/subMesh.template getSpaceStepsProducts< 1, 0 >() > 1.0)
- {
- //cout << "y = 0; x = n" << endl;
- u[i] = u[subMesh.getCellXPredecessor( i )] - subMesh.template getSpaceStepsProducts< 1, 0 >();
- }
- }
-
-
- else if(y == subMesh.getDimensions().y() - 1 && (boundaryCondition & 1) && x ==0)
- {
- if((u[subMesh.getCellXSuccessor( i )] - u[i])/subMesh.template getSpaceStepsProducts< 1, 0 >() > 1.0) {
- //cout << "y = n; x = 0" << endl;
- u[i] = u[subMesh.getCellXSuccessor( i )] - subMesh.template getSpaceStepsProducts< 1, 0 >();
- }
- }
- else if(y == subMesh.getDimensions().y() - 1 && (boundaryCondition & 1) && x == subMesh.getDimensions().x() - 1)
- {
- if((u[subMesh.getCellXPredecessor( i )] - u[i])/subMesh.template getSpaceStepsProducts< 1, 0 >() > 1.0)
- {
- //cout << "y = n; x = n" << endl;
- u[i] = u[subMesh.getCellXPredecessor( i )] - subMesh.template getSpaceStepsProducts< 1, 0 >();
- }
- }
- }*/
-
- /**/
-
-
-/* bool tmp = false;
- for(int i = 0; i < u.getSize(); i++)
- {
- if(u[0]*u[i] <= 0.0)
- tmp=true;
- }
-
-
- if(tmp)
- {}
- else if(boundaryCondition == 4)
- {
- int i;
- for(i = 0; i < u.getSize() - subMesh.getDimensions().x() ; i=subMesh.getCellYSuccessor(i))
- {
- int j;
- for(j = i; j < subMesh.getDimensions().x() - 1; j=subMesh.getCellXSuccessor(j))
- {
- u[j] = u[i];
- }
- u[j] = u[i];
- }
- int j;
- for(j = i; j < subMesh.getDimensions().x() - 1; j=subMesh.getCellXSuccessor(j))
- {
- u[j] = u[i];
- }
- u[j] = u[i];
- }
- else if(boundaryCondition == 8)
- {
- int i;
- for(i = 0; i < subMesh.getDimensions().x() - 1; i=subMesh.getCellXSuccessor(i))
- {
- int j;
- for(j = i; j < u.getSize() - subMesh.getDimensions().x(); j=subMesh.getCellYSuccessor(j))
- {
- u[j] = u[i];
- }
- u[j] = u[i];
- }
- int j;
- for(j = i; j < u.getSize() - subMesh.getDimensions().x(); j=subMesh.getCellYSuccessor(j))
- {
- u[j] = u[i];
- }
- u[j] = u[i];
-
- }
- else if(boundaryCondition == 2)
- {
- int i;
- for(i = subMesh.getDimensions().x() - 1; i < u.getSize() - subMesh.getDimensions().x() ; i=subMesh.getCellYSuccessor(i))
- {
- int j;
- for(j = i; j > (i-1)*subMesh.getDimensions().x(); j=subMesh.getCellXPredecessor(j))
- {
- u[j] = u[i];
- }
- u[j] = u[i];
- }
- int j;
- for(j = i; j > (i-1)*subMesh.getDimensions().x(); j=subMesh.getCellXPredecessor(j))
- {
- u[j] = u[i];
- }
- u[j] = u[i];
- }
- else if(boundaryCondition == 1)
- {
- int i;
- for(i = (subMesh.getDimensions().y() - 1)*subMesh.getDimensions().x(); i < u.getSize() - 1; i=subMesh.getCellXSuccessor(i))
- {
- int j;
- for(j = i; j >=subMesh.getDimensions().x(); j=subMesh.getCellYPredecessor(j))
- {
- u[j] = u[i];
- }
- u[j] = u[i];
- }
- int j;
- for(j = i; j >=subMesh.getDimensions().x(); j=subMesh.getCellYPredecessor(j))
- {
- u[j] = u[i];
- }
- u[j] = u[i];
- }
-*/
- /**/
-
bool tmp = false;
@@ -886,8 +777,6 @@ tnlParallelMapSolver<2,SchemeHost, SchemeDevice, Device, double, int>::runSubgri
if(this->unusedCell[centreGID] == 0 || boundaryCondition == 0)
tmp = true;
}
- //if(this->currentStep + 3 < getSubgridValue(subGridID))
- //tmp = true;
double value = Sign(u[0]) * u.absMax();
@@ -974,6 +863,15 @@ tnlParallelMapSolver<2,SchemeHost, SchemeDevice, Device, double, int>::runSubgri
//double lastResidue( 10000.0 );
tnlGridEntity<MeshType, 2, tnlGridEntityNoStencilStorage > Entity(subMesh);
tnlNeighbourGridEntityGetter<tnlGridEntity< MeshType, 2, tnlGridEntityNoStencilStorage >,2> neighbourEntities(Entity);
+
+ for( int i = 0; i < u.getSize(); i ++ )
+ {
+ if(map[i] == 0.0)
+ {
+ u[i] = /*Sign(u[l])**/MAP_SOLVER_MAX_VALUE;
+ }
+ }
+
while( time < finalTime /*|| maxResidue > subMesh.template getSpaceStepsProducts< 1, 0 >()*/)
{
/****
@@ -985,13 +883,14 @@ tnlParallelMapSolver<2,SchemeHost, SchemeDevice, Device, double, int>::runSubgri
Entity.setCoordinates(tnlStaticVector<2,int>(i % subMesh.getDimensions().x(),i / subMesh.getDimensions().x()));
Entity.refresh();
neighbourEntities.refresh(subMesh,Entity.getIndex());
- fu[ i ] = schemeHost.getValue( this->subMesh, i, tnlStaticVector<2,int>(i % subMesh.getDimensions().x(),i / subMesh.getDimensions().x()), u, time, boundaryCondition,neighbourEntities);
+ if(map[i] != 0.0)
+ fu[ i ] = schemeHost.getValue( this->subMesh, i, tnlStaticVector<2,int>(i % subMesh.getDimensions().x(),i / subMesh.getDimensions().x()), u, time, boundaryCondition,neighbourEntities,map);
}
maxResidue = fu. absMax();
- if( this -> cflCondition * maxResidue != 0.0)
- currentTau = this -> cflCondition / maxResidue;
+ if(maxResidue != 0.0)
+ currentTau = abs(this -> cflCondition / maxResidue);
/* if (maxResidue < 0.05)
cout << "Max < 0.05" << endl;*/
@@ -1016,9 +915,10 @@ tnlParallelMapSolver<2,SchemeHost, SchemeDevice, Device, double, int>::runSubgri
for( int i = 0; i < fu.getSize(); i ++ )
{
- double add = u[i] + currentTau * fu[ i ];
+// double add = u[i] + currentTau * fu[ i ];
//if( fabs(u[i]) < fabs(add) or (this->subgridValues[subGridID] == this->currentStep +4) )
- u[ i ] = add;
+ if(map[i] != 0.0)
+ u[ i ] += currentTau * fu[ i ];
}
time += currentTau;
@@ -1030,13 +930,13 @@ tnlParallelMapSolver<2,SchemeHost, SchemeDevice, Device, double, int>::runSubgri
/*if (u.max() > 0.0)
this->stopTime /=(double) this->gridCols;*/
- VectorType solution;
- solution.setLike(u);
- for( int i = 0; i < u.getSize(); i ++ )
- {
- solution[i]=u[i];
- }
- return solution;
+// VectorType solution;
+// solution.setLike(u);
+// for( int i = 0; i < u.getSize(); i ++ )
+// {
+// solution[i]=u[i];
+// }
+ return u;
}
diff --git a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver2D_impl.h b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver2D_impl.h
index 9912bed9a8..1a0869ab66 100644
--- a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver2D_impl.h
+++ b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver2D_impl.h
@@ -25,7 +25,7 @@ template< typename SchemeHost, typename SchemeDevice, typename Device>
tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>::tnlParallelEikonalSolver()
{
cout << "a" << endl;
- this->device = tnlCudaDevice; /////////////// tnlCuda Device --- vypocet na GPU, tnlHostDevice --- vypocet na CPU
+ this->device = tnlHostDevice; /////////////// tnlCuda Device --- vypocet na GPU, tnlHostDevice --- vypocet na CPU
#ifdef HAVE_CUDA
if(this->device == tnlCudaDevice)
@@ -168,7 +168,8 @@ bool tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>::
if(containsCurve)
{
//cout << "Computing initial SDF on subgrid " << i << "." << endl;
- insertSubgrid(runSubgrid(0, tmp[i],i), i);
+ tmp[i] = runSubgrid(0, tmp[i],i);
+ insertSubgrid(tmp[i], i);
setSubgridValue(i, 4);
//cout << "Computed initial SDF on subgrid " << i << "." << endl;
}
@@ -249,6 +250,8 @@ void tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>::
{
if(getSubgridValue(i) != INT_MAX)
{
+ VectorType tmp;
+ tmp.setSize(this->n * this->n);
//cout << "subMesh: " << i << ", BC: " << getBoundaryCondition(i) << endl;
if(getSubgridValue(i) == currentStep+4)
@@ -256,22 +259,30 @@ void tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>::
if(getBoundaryCondition(i) & 1)
{
- insertSubgrid( runSubgrid(1, getSubgrid(i),i), i);
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(1, tmp ,i);
+ insertSubgrid( tmp, i);
this->calculationsCount[i]++;
}
if(getBoundaryCondition(i) & 2)
{
- insertSubgrid( runSubgrid(2, getSubgrid(i),i), i);
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(1, tmp ,i);
+ insertSubgrid( tmp, 2);
this->calculationsCount[i]++;
}
if(getBoundaryCondition(i) & 4)
{
- insertSubgrid( runSubgrid(4, getSubgrid(i),i), i);
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(4, tmp ,i);
+ insertSubgrid( tmp, i);
this->calculationsCount[i]++;
}
if(getBoundaryCondition(i) & 8)
{
- insertSubgrid( runSubgrid(8, getSubgrid(i),i), i);
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(8, tmp ,i);
+ insertSubgrid( tmp, i);
this->calculationsCount[i]++;
}
}
@@ -280,25 +291,33 @@ void tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>::
/* &&(!(getBoundaryCondition(i) & 5) && !(getBoundaryCondition(i) & 10)) */)
{
//cout << "3 @ " << getBoundaryCondition(i) << endl;
- insertSubgrid( runSubgrid(3, getSubgrid(i),i), i);
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(1, tmp ,i);
+ insertSubgrid( tmp, 3);
}
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);
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(5, tmp ,i);
+ insertSubgrid( tmp, i);
}
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);
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(10, tmp ,i);
+ insertSubgrid( tmp, i);
}
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);
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(12, tmp ,i);
+ insertSubgrid( tmp, i);
}
@@ -659,7 +678,10 @@ void tnlParallelEikonalSolver<2,SchemeHost, SchemeDevice, Device, double, int>::
for( int j = 0; j < this->n*this->n; j++)
{
- int index = (i / this->gridCols)*this->n*this->n*this->gridCols + (i % this->gridCols)*this->n + (j/this->n)*this->n*this->gridCols + (j % this->n);
+ int index = (i / this->gridCols)*this->n*this->n*this->gridCols
+ + (i % this->gridCols)*this->n
+ + (j/this->n)*this->n*this->gridCols
+ + (j % this->n);
//OMP LOCK index
if( (fabs(this->work_u[index]) > fabs(u[j])) || (this->unusedCell[index] == 1) )
{
diff --git a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver3D_impl.h b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver3D_impl.h
index 3ea8d2e4e7..b7468205bb 100644
--- a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver3D_impl.h
+++ b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver3D_impl.h
@@ -148,33 +148,43 @@ bool tnlParallelEikonalSolver<3,SchemeHost, SchemeDevice, Device, double, int>::
if(this->device == tnlHostDevice)
{
- for(int i = 0; i < this->subgridValues.getSize(); i++)
- {
+ for(int i = 0; i < this->subgridValues.getSize(); i++)
+ {
+// cout << "Working on subgrid " << i <<" --- check 1" << endl;
- if(! tmp[i].setSize(this->n * this->n*this->n))
- cout << "Could not allocate tmp["<< i <<"] array." << endl;
- tmp[i] = getSubgrid(i);
- containsCurve = false;
+ if(! tmp[i].setSize(this->n*this->n*this->n))
+ cout << "Could not allocate tmp["<< i <<"] array." << endl;
+// cout << "Working on subgrid " << i <<" --- check 2" << endl;
- for(int j = 0; j < tmp[i].getSize(); j++)
- {
- if(tmp[i][0]*tmp[i][j] <= 0.0)
+ tmp[i] = getSubgrid(i);
+ containsCurve = false;
+// cout << "Working on subgrid " << i <<" --- check 3" << endl;
+
+
+ for(int j = 0; j < tmp[i].getSize(); j++)
{
- containsCurve = true;
- j=tmp[i].getSize();
+ if(tmp[i][0]*tmp[i][j] <= 0.0)
+ {
+ containsCurve = true;
+ j=tmp[i].getSize();
+// cout << tmp[i][0] << " " << tmp[i][j] << endl;
+ }
+
}
+// cout << "Working on subgrid " << i <<" --- check 4" << endl;
- }
- if(containsCurve)
- {
- //cout << "Computing initial SDF on subgrid " << i << "." << endl;
- insertSubgrid(runSubgrid(0, tmp[i],i), i);
- setSubgridValue(i, 4);
- //cout << "Computed initial SDF on subgrid " << i << "." << endl;
- }
- containsCurve = false;
+ if(containsCurve)
+ {
+// cout << "Computing initial SDF on subgrid " << i << "." << endl;
+ tmp[i] = runSubgrid(0, tmp[i] ,i);
+ insertSubgrid( tmp[i], i);
+ setSubgridValue(i, 4);
+// cout << "Computed initial SDF on subgrid " << i << "." << endl;
+ }
+ containsCurve = false;
- }
+ }
+// cout << "CPU: Curve found" << endl;
}
#ifdef HAVE_CUDA
else if(this->device == tnlCudaDevice)
@@ -247,6 +257,8 @@ void tnlParallelEikonalSolver<3,SchemeHost, SchemeDevice, Device, double, int>::
#endif
for(int i = 0; i < this->subgridValues.getSize(); i++)
{
+ VectorType tmp;
+ tmp.setSize(this->n*this->n*this->n);
if(getSubgridValue(i) != INT_MAX)
{
//cout << "subMesh: " << i << ", BC: " << getBoundaryCondition(i) << endl;
@@ -254,61 +266,102 @@ void tnlParallelEikonalSolver<3,SchemeHost, SchemeDevice, Device, double, int>::
if(getSubgridValue(i) == currentStep+4)
{
- if(getBoundaryCondition(i) & 1)
- {
- insertSubgrid( runSubgrid(1, getSubgrid(i),i), i);
- this->calculationsCount[i]++;
+ if(getBoundaryCondition(i) & 1)
+ {
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(1, tmp ,i);
+ insertSubgrid( tmp, i);
+ this->calculationsCount[i]++;
+ }
+ if(getBoundaryCondition(i) & 2)
+ {
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(2, tmp ,i);
+ insertSubgrid( tmp, i);
+ this->calculationsCount[i]++;
+ }
+ if(getBoundaryCondition(i) & 4)
+ {
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(4, tmp ,i);
+ insertSubgrid( tmp, i);
+ this->calculationsCount[i]++;
+ }
+ if(getBoundaryCondition(i) & 8)
+ {
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(8, tmp ,i);
+ insertSubgrid( tmp, i);
+ this->calculationsCount[i]++;
+ }
+ if(getBoundaryCondition(i) & 16)
+ {
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(16, tmp ,i);
+ insertSubgrid( tmp, i);
+ this->calculationsCount[i]++;
+ }
+ if(getBoundaryCondition(i) & 32)
+ {
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(32, tmp ,i);
+ insertSubgrid( tmp, i);
+ this->calculationsCount[i]++;
+ }
}
- if(getBoundaryCondition(i) & 2)
+
+ if( getBoundaryCondition(i) & 19)
{
- insertSubgrid( runSubgrid(2, getSubgrid(i),i), i);
- this->calculationsCount[i]++;
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(19, tmp ,i);
+ insertSubgrid( tmp, i);
}
- if(getBoundaryCondition(i) & 4)
+ if( getBoundaryCondition(i) & 21)
{
- insertSubgrid( runSubgrid(4, getSubgrid(i),i), i);
- this->calculationsCount[i]++;
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(21, tmp ,i);
+ insertSubgrid( tmp, i);
}
- if(getBoundaryCondition(i) & 8)
+ if( getBoundaryCondition(i) & 26)
{
- insertSubgrid( runSubgrid(8, getSubgrid(i),i), i);
- this->calculationsCount[i]++;
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(26, tmp ,i);
+ insertSubgrid( tmp, i);
}
+ if( getBoundaryCondition(i) & 28)
+ {
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(28, tmp ,i);
+ insertSubgrid( tmp, i);
}
- if( ((getBoundaryCondition(i) & 2) )|| (getBoundaryCondition(i) & 1)//)
- /* &&(!(getBoundaryCondition(i) & 5) && !(getBoundaryCondition(i) & 10)) */)
+ if( getBoundaryCondition(i) & 35)
{
- //cout << "3 @ " << getBoundaryCondition(i) << endl;
- insertSubgrid( runSubgrid(3, getSubgrid(i),i), i);
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(35, tmp ,i);
+ insertSubgrid( tmp, i);
}
- if( ((getBoundaryCondition(i) & 4) )|| (getBoundaryCondition(i) & 1)//)
- /* &&(!(getBoundaryCondition(i) & 3) && !(getBoundaryCondition(i) & 12)) */)
+ if( getBoundaryCondition(i) & 37)
{
- //cout << "5 @ " << getBoundaryCondition(i) << endl;
- insertSubgrid( runSubgrid(5, getSubgrid(i),i), i);
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(37, tmp ,i);
+ insertSubgrid( tmp, i);
}
- if( ((getBoundaryCondition(i) & 2) )|| (getBoundaryCondition(i) & 8)//)
- /* &&(!(getBoundaryCondition(i) & 12) && !(getBoundaryCondition(i) & 3))*/ )
+ if( getBoundaryCondition(i) & 42)
{
- //cout << "10 @ " << getBoundaryCondition(i) << endl;
- insertSubgrid( runSubgrid(10, getSubgrid(i),i), i);
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(42, tmp ,i);
+ insertSubgrid( tmp, i);
}
- if( ((getBoundaryCondition(i) & 4) )|| (getBoundaryCondition(i) & 8)//)
- /*&&(!(getBoundaryCondition(i) & 10) && !(getBoundaryCondition(i) & 5)) */)
+ if( getBoundaryCondition(i) & 44)
{
- //cout << "12 @ " << getBoundaryCondition(i) << endl;
- insertSubgrid( runSubgrid(12, getSubgrid(i),i), i);
+ tmp = getSubgrid(i);
+ tmp = runSubgrid(44, tmp ,i);
+ insertSubgrid( tmp, i);
}
- /*if(getBoundaryCondition(i))
- {
- insertSubgrid( runSubgrid(15, getSubgrid(i),i), i);
- }*/
-
setBoundaryCondition(i, 0);
-
setSubgridValue(i, getSubgridValue(i)-1);
}
@@ -413,81 +466,147 @@ void tnlParallelEikonalSolver<3,SchemeHost, SchemeDevice, Device, double, int>::
int tmp1, tmp2;
int grid1, grid2;
- if(this->currentStep & 1)
- {
+// if(this->currentStep & 1)
+// {
for(int j = 0; j < this->gridRows - 1; j++)
{
for (int i = 0; i < this->gridCols*this->n; i++)
{
- tmp1 = this->gridCols*this->n*((this->n-1)+j*this->n) + i;
- 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))
+ for (int k = 0; k < this->gridLevels*this->n; k++)
{
- this->work_u[tmp2] = this->work_u[tmp1];
- this->unusedCell[tmp2] = 0;
- if(grid2 == INT_MAX)
+// cout << "a" << endl;
+ tmp1 = this->gridCols*this->n*((this->n-1)+j*this->n) + i + k*this->gridCols*this->n*this->gridRows*this->n;
+// cout << "b" << endl;
+ tmp2 = this->gridCols*this->n*((this->n)+j*this->n) + i + k*this->gridCols*this->n*this->gridRows*this->n;
+// cout << "c" << endl;
+ if(tmp1 > work_u.getSize())
+ cout << "tmp1: " << tmp1 << " x: " << j <<" y: " << i <<" z: " << k << endl;
+ if(tmp2 > work_u.getSize())
+ cout << "tmp2: " << tmp2 << " x: " << j <<" y: " << i <<" z: " << k << endl;
+ grid1 = getSubgridValue(getOwner(tmp1));
+// cout << "d" << endl;
+ grid2 = getSubgridValue(getOwner(tmp2));
+// cout << "e" << endl;
+ if(getOwner(tmp1)==getOwner(tmp2))
+ cout << "i, j, k" << i << "," << j << "," << k << 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))
{
- setSubgridValue(getOwner(tmp2), -INT_MAX);
+ this->work_u[tmp2] = this->work_u[tmp1];
+// cout << "f" << endl;
+ this->unusedCell[tmp2] = 0;
+// cout << "g" << endl;
+ if(grid2 == INT_MAX)
+ {
+ setSubgridValue(getOwner(tmp2), -INT_MAX);
+ }
+// cout << "h" << endl;
+ if(! (getBoundaryCondition(getOwner(tmp2)) & 8) )
+ setBoundaryCondition(getOwner(tmp2), getBoundaryCondition(getOwner(tmp2))+8);
+// cout << "i" << endl;
}
- if(! (getBoundaryCondition(getOwner(tmp2)) & 8) )
- setBoundaryCondition(getOwner(tmp2), getBoundaryCondition(getOwner(tmp2))+8);
- }
- else if ((fabs(this->work_u[tmp1]) > fabs(this->work_u[tmp2]) + this->delta || grid1 == INT_MAX || grid1 == -INT_MAX) && (grid2 != INT_MAX && grid2 != -INT_MAX))
- {
- this->work_u[tmp1] = this->work_u[tmp2];
- this->unusedCell[tmp1] = 0;
- if(grid1 == INT_MAX)
+ else if ((fabs(this->work_u[tmp1]) > fabs(this->work_u[tmp2]) + this->delta || grid1 == INT_MAX || grid1 == -INT_MAX) && (grid2 != INT_MAX && grid2 != -INT_MAX))
{
- setSubgridValue(getOwner(tmp1), -INT_MAX);
+ this->work_u[tmp1] = this->work_u[tmp2];
+// cout << "j" << endl;
+ this->unusedCell[tmp1] = 0;
+// cout << "k" << endl;
+ if(grid1 == INT_MAX)
+ {
+ setSubgridValue(getOwner(tmp1), -INT_MAX);
+ }
+// cout << "l" << endl;
+ if(! (getBoundaryCondition(getOwner(tmp1)) & 1) )
+ setBoundaryCondition(getOwner(tmp1), getBoundaryCondition(getOwner(tmp1))+1);
+// cout << "m" << endl;
}
- if(! (getBoundaryCondition(getOwner(tmp1)) & 1) )
- setBoundaryCondition(getOwner(tmp1), getBoundaryCondition(getOwner(tmp1))+1);
}
}
}
- }
- else
- {
+// }
+// else
+// {
+
+ cout << "sync 2" << endl;
for(int i = 1; i < this->gridCols; i++)
{
for (int j = 0; j < this->gridRows*this->n; j++)
{
- tmp1 = this->gridCols*this->n*j + i*this->n - 1;
- 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))
+ for (int k = 0; k < this->gridLevels*this->n; k++)
{
- this->work_u[tmp2] = this->work_u[tmp1];
- this->unusedCell[tmp2] = 0;
- if(grid2 == INT_MAX)
+ tmp1 = this->gridCols*this->n*j + i*this->n - 1 + k*this->gridCols*this->n*this->gridRows*this->n;
+ tmp2 = this->gridCols*this->n*j + i*this->n + k*this->gridCols*this->n*this->gridRows*this->n;
+ grid1 = getSubgridValue(getOwner(tmp1));
+ grid2 = getSubgridValue(getOwner(tmp2));
+ if(getOwner(tmp1)==getOwner(tmp2))
+ cout << "i, j, k" << i << "," << j << "," << k << 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))
{
- setSubgridValue(getOwner(tmp2), -INT_MAX);
+ this->work_u[tmp2] = this->work_u[tmp1];
+ this->unusedCell[tmp2] = 0;
+ if(grid2 == INT_MAX)
+ {
+ setSubgridValue(getOwner(tmp2), -INT_MAX);
+ }
+ if(! (getBoundaryCondition(getOwner(tmp2)) & 4) )
+ setBoundaryCondition(getOwner(tmp2), getBoundaryCondition(getOwner(tmp2))+4);
+ }
+ else if ((fabs(this->work_u[tmp1]) > fabs(this->work_u[tmp2]) + this->delta || grid1 == INT_MAX || grid1 == -INT_MAX) && (grid2 != INT_MAX && grid2 != -INT_MAX))
+ {
+ this->work_u[tmp1] = this->work_u[tmp2];
+ this->unusedCell[tmp1] = 0;
+ if(grid1 == INT_MAX)
+ {
+ setSubgridValue(getOwner(tmp1), -INT_MAX);
+ }
+ if(! (getBoundaryCondition(getOwner(tmp1)) & 2) )
+ setBoundaryCondition(getOwner(tmp1), getBoundaryCondition(getOwner(tmp1))+2);
}
- if(! (getBoundaryCondition(getOwner(tmp2)) & 4) )
- setBoundaryCondition(getOwner(tmp2), getBoundaryCondition(getOwner(tmp2))+4);
}
- else if ((fabs(this->work_u[tmp1]) > fabs(this->work_u[tmp2]) + this->delta || grid1 == INT_MAX || grid1 == -INT_MAX) && (grid2 != INT_MAX && grid2 != -INT_MAX))
+ }
+ }
+
+ cout << "sync 3" << endl;
+
+ for(int k = 1; k < this->gridLevels; k++)
+ {
+ for (int j = 0; j < this->gridRows*this->n; j++)
+ {
+ for (int i = 0; i < this->gridCols*this->n; i++)
{
- this->work_u[tmp1] = this->work_u[tmp2];
- this->unusedCell[tmp1] = 0;
- if(grid1 == INT_MAX)
+ tmp1 = this->gridCols*this->n*j + i + (k*this->n-1)*this->gridCols*this->n*this->gridRows*this->n;
+ tmp2 = this->gridCols*this->n*j + i + k*this->n*this->gridCols*this->n*this->gridRows*this->n;
+ grid1 = getSubgridValue(getOwner(tmp1));
+ grid2 = getSubgridValue(getOwner(tmp2));
+ if(getOwner(tmp1)==getOwner(tmp2))
+ cout << "i, j, k" << i << "," << j << "," << k << 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];
+ this->unusedCell[tmp2] = 0;
+ if(grid2 == INT_MAX)
+ {
+ setSubgridValue(getOwner(tmp2), -INT_MAX);
+ }
+ if(! (getBoundaryCondition(getOwner(tmp2)) & 32) )
+ setBoundaryCondition(getOwner(tmp2), getBoundaryCondition(getOwner(tmp2))+32);
+ }
+ else if ((fabs(this->work_u[tmp1]) > fabs(this->work_u[tmp2]) + this->delta || grid1 == INT_MAX || grid1 == -INT_MAX) && (grid2 != INT_MAX && grid2 != -INT_MAX))
{
- setSubgridValue(getOwner(tmp1), -INT_MAX);
+ this->work_u[tmp1] = this->work_u[tmp2];
+ this->unusedCell[tmp1] = 0;
+ if(grid1 == INT_MAX)
+ {
+ setSubgridValue(getOwner(tmp1), -INT_MAX);
+ }
+ if(! (getBoundaryCondition(getOwner(tmp1)) & 16) )
+ setBoundaryCondition(getOwner(tmp1), getBoundaryCondition(getOwner(tmp1))+16);
}
- if(! (getBoundaryCondition(getOwner(tmp1)) & 2) )
- setBoundaryCondition(getOwner(tmp1), getBoundaryCondition(getOwner(tmp1))+2);
}
}
}
- }
+// }
+
this->currentStep++;
@@ -507,7 +626,11 @@ template< typename SchemeHost, typename SchemeDevice, typename Device>
int tnlParallelEikonalSolver<3,SchemeHost, SchemeDevice, Device, double, int>::getOwner(int i) const
{
- return (i / (this->gridCols*this->n*this->n))*this->gridCols + (i % (this->gridCols*this->n))/this->n;
+ int j = i % (this->gridCols*this->gridRows*this->n*this->n);
+
+ return ( (i / (this->gridCols*this->gridRows*this->n*this->n*this->n))*this->gridCols*this->gridRows
+ + (j / (this->gridCols*this->n*this->n))*this->gridCols
+ + (j % (this->gridCols*this->n))/this->n);
}
template< typename SchemeHost, typename SchemeDevice, typename Device>
@@ -648,15 +771,31 @@ template< typename SchemeHost, typename SchemeDevice, typename Device>
typename tnlParallelEikonalSolver<3,SchemeHost, SchemeDevice, Device, double, int>::VectorType
tnlParallelEikonalSolver<3,SchemeHost, SchemeDevice, Device, double, int>::getSubgrid( const int i ) const
{
+
VectorType u;
- u.setSize(this->n*this->n);
+ u.setSize(this->n*this->n*this->n);
+
+ int idx, idy, idz;
+ idz = i / (gridRows*this->gridCols);
+ idy = (i % (this->gridRows*this->gridCols)) / this->gridCols;
+ idx = i % (this->gridCols);
- for( int j = 0; j < u.getSize(); j++)
+ for( int j = 0; j < this->n; j++)
{
- u[j] = this->work_u[ (i / this->gridCols) * this->n*this->n*this->gridCols
- + (i % this->gridCols) * this->n
- + (j/this->n) * this->n*this->gridCols
- + (j % this->n) ];
+ // int index = (i / this->gridCols)*this->n*this->n*this->gridCols + (i % this->gridCols)*this->n + (j/this->n)*this->n*this->gridCols + (j % this->n);
+ for( int k = 0; k < this->n; k++)
+ {
+ for( int l = 0; l < this->n; l++)
+ {
+ int index = (idz*this->n + l) * this->n*this->n*this->gridCols*this->gridRows
+ + (idy) * this->n*this->n*this->gridCols
+ + (idx) * this->n
+ + k * this->n*this->gridCols
+ + j;
+
+ u[j + k*this->n + l*this->n*this->n] = this->work_u[ index ];
+ }
+ }
}
return u;
}
@@ -664,17 +803,35 @@ tnlParallelEikonalSolver<3,SchemeHost, SchemeDevice, Device, double, int>::getSu
template< typename SchemeHost, typename SchemeDevice, typename Device>
void tnlParallelEikonalSolver<3,SchemeHost, SchemeDevice, Device, double, int>::insertSubgrid( VectorType u, const int i )
{
+ int idx, idy, idz;
+ idz = i / (this->gridRows*this->gridCols);
+ idy = (i % (this->gridRows*this->gridCols)) / this->gridCols;
+ idx = i % (this->gridCols);
- for( int j = 0; j < this->n*this->n; j++)
+ for( int j = 0; j < this->n; j++)
{
- int index = (i / this->gridCols)*this->n*this->n*this->gridCols + (i % this->gridCols)*this->n + (j/this->n)*this->n*this->gridCols + (j % this->n);
- //OMP LOCK index
- if( (fabs(this->work_u[index]) > fabs(u[j])) || (this->unusedCell[index] == 1) )
+ // int index = (i / this->gridCols)*this->n*this->n*this->gridCols + (i % this->gridCols)*this->n + (j/this->n)*this->n*this->gridCols + (j % this->n);
+ for( int k = 0; k < this->n; k++)
{
- this->work_u[index] = u[j];
- this->unusedCell[index] = 0;
+ for( int l = 0; l < this->n; l++)
+ {
+
+ int index = (idz*this->n + l) * this->n*this->n*this->gridCols*this->gridRows
+ + (idy) * this->n*this->n*this->gridCols
+ + (idx) * this->n
+ + k * this->n*this->gridCols
+ + j;
+
+ //OMP LOCK index
+// cout<< idx << " " << idy << " " << idz << " " << j << " " << k << " " << l << " " << idz << " " << unusedCell.getSize() << " " << u.getSize() << " " << index <<endl;
+ if( (fabs(this->work_u[index]) > fabs(u[j + k*this->n + l*this->n*this->n])) || (this->unusedCell[index] == 1) )
+ {
+ this->work_u[index] = u[j + k*this->n + l*this->n*this->n];
+ this->unusedCell[index] = 0;
+ }
+ //OMP UNLOCK index
+ }
}
- //OMP UNLOCK index
}
}
@@ -694,10 +851,15 @@ tnlParallelEikonalSolver<3,SchemeHost, SchemeDevice, Device, double, int>::runSu
{
if(u[0]*u[i] <= 0.0)
tmp=true;
- int centreGID = (this->n*(subGridID / this->gridRows)+ (this->n >> 1))*(this->n*this->gridCols) + this->n*(subGridID % this->gridRows) + (this->n >> 1);
- if(this->unusedCell[centreGID] == 0 || boundaryCondition == 0)
- tmp = true;
}
+ int idx,idy,idz;
+ idz = subGridID / (this->gridRows*this->gridCols);
+ idy = (subGridID % (this->gridRows*this->gridCols)) / this->gridCols;
+ idx = subGridID % (this->gridCols);
+ int centreGID = (this->n*idy + (this->n>>1) )*(this->n*this->gridCols) + this->n*idx + (this->n>>1)
+ + ((this->n>>1)+this->n*idz)*this->n*this->n*this->gridRows*this->gridCols;
+ if(this->unusedCell[centreGID] == 0 || boundaryCondition == 0)
+ tmp = true;
//if(this->currentStep + 3 < getSubgridValue(subGridID))
//tmp = true;
@@ -713,35 +875,56 @@ tnlParallelEikonalSolver<3,SchemeHost, SchemeDevice, Device, double, int>::runSu
{
for(int i = 0; i < this->n; i++)
for(int j = 1;j < this->n; j++)
+ for(int k = 0;k < this->n; k++)
//if(fabs(u[i*this->n + j]) < fabs(u[i*this->n]))
- u[i*this->n + j] = value;// u[i*this->n];
+ u[k*this->n*this->n + i*this->n + j] = value;// u[i*this->n];
}
else if(boundaryCondition == 2)
{
for(int i = 0; i < this->n; i++)
for(int j =0 ;j < this->n -1; j++)
+ for(int k = 0;k < this->n; k++)
//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];
+ u[k*this->n*this->n + i*this->n + j] = value;// u[(i+1)*this->n - 1];
}
else if(boundaryCondition == 1)
{
for(int j = 0; j < this->n; j++)
for(int i = 0;i < this->n - 1; i++)
+ for(int k = 0;k < this->n; k++)
//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)];
+ u[k*this->n*this->n + i*this->n + j] = value;// u[j + this->n*(this->n - 1)];
}
else if(boundaryCondition == 8)
{
for(int j = 0; j < this->n; j++)
for(int i = 1;i < this->n; i++)
+ for(int k = 0;k < this->n; k++)
+ //if(fabs(u[i*this->n + j]) < fabs(u[j]))
+ u[k*this->n*this->n + i*this->n + j] = value;// u[j];
+ }
+ else if(boundaryCondition == 16)
+ {
+ for(int j = 0; j < this->n; j++)
+ for(int i = 0;i < this->n ; i++)
+ for(int k = 0;k < this->n-1; k++)
+ //if(fabs(u[i*this->n + j]) < fabs(u[j + this->n*(this->n - 1)]))
+ u[k*this->n*this->n + i*this->n + j] = value;// u[j + this->n*(this->n - 1)];
+ }
+ else if(boundaryCondition == 32)
+ {
+ for(int j = 0; j < this->n; j++)
+ for(int i = 0;i < this->n; i++)
+ for(int k = 1;k < this->n; k++)
//if(fabs(u[i*this->n + j]) < fabs(u[j]))
- u[i*this->n + j] = value;// u[j];
+ u[k*this->n*this->n + i*this->n + j] = value;// u[j];
}
double time = 0.0;
double currentTau = this->tau0;
double finalTime = this->stopTime;// + 3.0*(u.max() - u.min());
+ if(boundaryCondition == 0) finalTime *= 2.0;
if( time + currentTau > finalTime ) currentTau = finalTime - time;
double maxResidue( 1.0 );
@@ -756,13 +939,19 @@ tnlParallelEikonalSolver<3,SchemeHost, SchemeDevice, Device, double, int>::runSu
for( int i = 0; i < fu.getSize(); i ++ )
{
+// cout << "i: " << i << ", time: " << time <<endl;
tnlStaticVector<3,int> coords(i % subMesh.getDimensions().x(),
- (i % subMesh.getDimensions().x()*subMesh.getDimensions().y())/ subMesh.getDimensions().x(),
- i / subMesh.getDimensions().x()*subMesh.getDimensions().y());
+ (i % (subMesh.getDimensions().x()*subMesh.getDimensions().y())) / subMesh.getDimensions().x(),
+ i / (subMesh.getDimensions().x()*subMesh.getDimensions().y()));
+// cout << "b " << i << " " << i % subMesh.getDimensions().x() << " " << (i % (subMesh.getDimensions().x()*subMesh.getDimensions().y())) << " " << (i % subMesh.getDimensions().x()*subMesh.getDimensions().y()) / subMesh.getDimensions().x() << " " << subMesh.getDimensions().x()*subMesh.getDimensions().y() << " " <<endl;
Entity.setCoordinates(coords);
+// cout <<"c" << coords << endl;
Entity.refresh();
+// cout << "d" <<endl;
neighbourEntities.refresh(subMesh,Entity.getIndex());
+// cout << "e" <<endl;
fu[ i ] = schemeHost.getValue( this->subMesh, i, coords,u, time, boundaryCondition, neighbourEntities );
+// cout << "f" <<endl;
}
maxResidue = fu. absMax();
@@ -772,11 +961,8 @@ tnlParallelEikonalSolver<3,SchemeHost, SchemeDevice, Device, double, int>::runSu
/* if (maxResidue < 0.05)
cout << "Max < 0.05" << endl;*/
- if(currentTau > 1.0 * this->subMesh.template getSpaceStepsProducts< 1, 0, 0 >())
- {
- //cout << currentTau << " >= " << 2.0 * this->subMesh.template getSpaceStepsProducts< 1, 0, 0 >() << endl;
- currentTau = 1.0 * this->subMesh.template getSpaceStepsProducts< 1, 0, 0 >();
- }
+ if(currentTau > 0.5 * this->subMesh.template getSpaceStepsProducts< 1, 0, 0 >())
+ currentTau = 0.5 * this->subMesh.template getSpaceStepsProducts< 1, 0, 0 >();
/*if(maxResidue > lastResidue)
currentTau *=(1.0/10.0);*/
@@ -807,13 +993,14 @@ tnlParallelEikonalSolver<3,SchemeHost, SchemeDevice, Device, double, int>::runSu
/*if (u.max() > 0.0)
this->stopTime /=(double) this->gridCols;*/
- VectorType solution;
- solution.setLike(u);
- for( int i = 0; i < u.getSize(); i ++ )
- {
- solution[i]=u[i];
- }
- return solution;
+// VectorType solution;
+// solution.setLike(u);
+// for( int i = 0; i < u.getSize(); i ++ )
+// {
+// solution[i]=u[i];
+// }
+// return solution;
+ return u;
}
diff --git a/src/operators/godunov-eikonal/parallelGodunovEikonal3D_impl.h b/src/operators/godunov-eikonal/parallelGodunovEikonal3D_impl.h
index 7206cc9a40..49cb75cd4b 100644
--- a/src/operators/godunov-eikonal/parallelGodunovEikonal3D_impl.h
+++ b/src/operators/godunov-eikonal/parallelGodunovEikonal3D_impl.h
@@ -145,9 +145,7 @@ Real parallelGodunovEikonalScheme< tnlGrid< 3, MeshReal, Device, MeshIndex >, Re
//-----------------------------------
RealType signui;
- signui = sign(u[cellIndex],this->epsilon);
-
-
+ signui = sign(u[cellIndex], this->epsilon);
RealType xb = u[cellIndex];
@@ -190,12 +188,6 @@ Real parallelGodunovEikonalScheme< tnlGrid< 3, MeshReal, Device, MeshIndex >, Re
else
zb -= u[neighbourEntities.template getEntityIndex< 0, 0, -1 >()];
-
- //xb *= ihx;
- //xf *= ihx;
- // yb *= ihy;
- //yf *= ihy;
-
if(signui > 0.0)
{
xf = negativePart(xf);
@@ -243,7 +235,9 @@ Real parallelGodunovEikonalScheme< tnlGrid< 3, MeshReal, Device, MeshIndex >, Re
else
c = zf;
- d =(1.0 - sqrt(a*a+b*b+c*c)*ihx );
+ d = ( 1.0 - sqrt(a*a + b*b + c*c)*ihx );
+
+// d = 1.0 - sqrt(xf*xf + xb*xb + yf*yf + yb*yb + zf*zf + zb*zb)*ihx; /*upwind*/
if(Sign(d) > 0.0 )
return Sign(u[cellIndex])*d;
diff --git a/src/operators/godunov-eikonal/parallelGodunovMap.h b/src/operators/godunov-eikonal/parallelGodunovMap.h
index 19a764d9ba..50313540ab 100644
--- a/src/operators/godunov-eikonal/parallelGodunovMap.h
+++ b/src/operators/godunov-eikonal/parallelGodunovMap.h
@@ -150,7 +150,8 @@ public:
const Vector& u,
const RealType& time,
const IndexType boundaryCondition,
- const tnlNeighbourGridEntityGetter<tnlGridEntity< MeshType, 2, tnlGridEntityNoStencilStorage >,2> neighbourEntities) const;
+ const tnlNeighbourGridEntityGetter<tnlGridEntity< MeshType, 2, tnlGridEntityNoStencilStorage >,2> neighbourEntities,
+ const Vector& map) const;
#ifdef HAVE_CUDA
__device__
diff --git a/src/operators/godunov-eikonal/parallelGodunovMap2D_impl.h b/src/operators/godunov-eikonal/parallelGodunovMap2D_impl.h
index 3be04db3fd..fc26c698b6 100644
--- a/src/operators/godunov-eikonal/parallelGodunovMap2D_impl.h
+++ b/src/operators/godunov-eikonal/parallelGodunovMap2D_impl.h
@@ -145,7 +145,8 @@ Real parallelGodunovMapScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real,
const Vector& u,
const Real& time,
const IndexType boundaryCondition,
- const tnlNeighbourGridEntityGetter<tnlGridEntity< MeshType, 2, tnlGridEntityNoStencilStorage >,2> neighbourEntities ) const
+ const tnlNeighbourGridEntityGetter<tnlGridEntity< MeshType, 2, tnlGridEntityNoStencilStorage >,2> neighbourEntities,
+ const Vector& map) const
{
@@ -164,15 +165,22 @@ Real parallelGodunovMapScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real,
}
- //-----------------------------------
-
RealType signui;
- signui = sign(u[cellIndex],this->epsilon);
+// if(boundaryCondition == 0)
+ signui = sign(u[cellIndex],/*(boundaryCondition != 0) * */this->epsilon);
+// else
+// signui = Sign(u[cellIndex]);
+ RealType value;
+// if(map[cellIndex] == 0.0)
+// {
+//// value = INT_MAX;
+// u[cellIndex] = Sign(u[cellIndex])*INT_MAX;
+// return 0.0;
+// }
+// else
+ value = 50.0/ map[cellIndex];
-#ifdef HAVE_CUDA
- //printf("%d : %d ;;;; %d : %d , %f \n",threadIdx.x, mesh.getDimensions().x() , threadIdx.y,mesh.getDimensions().y(), epsilon );
-#endif
RealType xb = u[cellIndex];
RealType xf = -u[cellIndex];
@@ -202,11 +210,6 @@ Real parallelGodunovMapScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real,
yb -= u[neighbourEntities.template getEntityIndex< 0, -1 >()];
- //xb *= ihx;
- //xf *= ihx;
- // yb *= ihy;
- //yf *= ihy;
-
if(signui > 0.0)
{
xf = negativePart(xf);
@@ -230,6 +233,15 @@ Real parallelGodunovMapScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real,
yf = positivePart(yf);
}
+ /*if(boundaryCondition &1)
+ yb = 0.0;
+ else
+ yf = 0.0;
+
+ if(boundaryCondition & 2)
+ xb = 0.0;
+ else
+ xf = 0.0;*/
if(xb - xf > 0.0)
a = xb;
@@ -241,140 +253,13 @@ Real parallelGodunovMapScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real,
else
b = yf;
- c =(1.0 - sqrt(a*a+b*b)*ihx );
+ c =(value - sqrt(a*a+b*b)*ihx );
- if(Sign(c) > 0.0 )
+ if(c > 0.0 )
return Sign(u[cellIndex])*c;
else
- return signui*c;
-
- //--------------------------------------------------
-
-//
-//
-//
-// RealType acc = hx*hy*hx*hy;
-//
-// RealType nabla, xb, xf, yb, yf, signui;
-//
-// signui = sign(u[cellIndex],this->epsilon);
-//
-//
-// //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[neighbourEntities.template getEntityIndex< -1, 0 >()] - u[cellIndex])*ihx);
-// else
-// xf = negativePart((u[neighbourEntities.template getEntityIndex< 1, 0 >()] - u[cellIndex])*ihx);
-//
-// /**/ /* if(boundaryCondition & 4)
-// xb = (u[cellIndex] - u[mesh.getCellXPredecessor( cellIndex )])*ihx;
-// else *//*if(boundaryCondition & 2)
-// xb = 0.0;
-// else /**/if(coordinates.x() == 0)
-// xb = positivePart((u[cellIndex] - u[mesh.template getCellNextToCell<+1,0>( cellIndex )])*ihx);
-// else
-// xb = positivePart((u[cellIndex] - u[neighbourEntities.template getEntityIndex< -1, 0 >()])*ihx);
-//
-// /**/ /* if(boundaryCondition & 1)
-// yf = (u[mesh.getCellYSuccessor( cellIndex )] - u[cellIndex])*ihy;
-// else *//*if(boundaryCondition & 8)
-// yf = 0.0;
-// else /**/if(coordinates.y() == mesh.getDimensions().y() - 1)
-// yf = negativePart((u[neighbourEntities.template getEntityIndex< 0, -1 >()] - u[cellIndex])*ihy);
-// else
-// yf = negativePart((u[neighbourEntities.template getEntityIndex< 0, 1 >()] - u[cellIndex])*ihy);
-//
-// /**/ /* if(boundaryCondition & 8)
-// yb = (u[cellIndex] - u[mesh.getCellYPredecessor( cellIndex )])*ihy;
-// else *//*if(boundaryCondition & 1)
-// yb = 0.0;
-// else /**/if(coordinates.y() == 0)
-// yb = positivePart((u[cellIndex] - u[neighbourEntities.template getEntityIndex< 0, 1 >()])*ihy);
-// else
-// yb = positivePart((u[cellIndex] - u[neighbourEntities.template getEntityIndex< 0, -1 >()])*ihy);
-//
-// if(xb - xf > 0.0)
-// xf = 0.0;
-// else
-// xb = 0.0;
-//
-// if(yb - yf > 0.0)
-// yf = 0.0;
-// 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])*ihx;
-// else*//* if(boundaryCondition & 4)
-// xf = 0.0;
-// else /**/if(coordinates.x() == mesh.getDimensions().x() - 1)
-// xf = positivePart((u[neighbourEntities.template getEntityIndex< -1, 0 >()] - u[cellIndex])*ihx);
-// else
-// xf = positivePart((u[neighbourEntities.template getEntityIndex< 1, 0 >()] - u[cellIndex])*ihx);
-//
-// /**/ /* if(boundaryCondition & 4)
-// xb = (u[cellIndex] - u[mesh.getCellXPredecessor( cellIndex )])*ihx;
-// else*//* if(boundaryCondition & 2)
-// xb = 0.0;
-// else /**/if(coordinates.x() == 0)
-// xb = negativePart((u[cellIndex] - u[neighbourEntities.template getEntityIndex< 1, 0 >()])*ihx);
-// else
-// xb = negativePart((u[cellIndex] - u[neighbourEntities.template getEntityIndex< -1, 0 >()])*ihx);
-//
-// /**/ /* if(boundaryCondition & 1)
-// yf = (u[mesh.getCellYSuccessor( cellIndex )] - u[cellIndex])*ihy;
-// else *//*if(boundaryCondition & 8)
-// yf = 0.0;
-// else /**/if(coordinates.y() == mesh.getDimensions().y() - 1)
-// yf = positivePart((u[neighbourEntities.template getEntityIndex< 0, -1 >()] - u[cellIndex])*ihy);
-// else
-// yf = positivePart((u[neighbourEntities.template getEntityIndex< 0, 1 >()] - u[cellIndex])*ihy);
-//
-// /**/ /* if(boundaryCondition & 8)
-// yb = (u[cellIndex] - u[mesh.getCellYPredecessor( cellIndex )])*ihy;
-// else*//* if(boundaryCondition & 1)
-// yb = 0.0;
-// else /**/if(coordinates.y() == 0)
-// yb = negativePart((u[cellIndex] - u[neighbourEntities.template getEntityIndex< 0, 1 >()])*ihy);
-// else
-// yb = negativePart((u[cellIndex] - u[neighbourEntities.template getEntityIndex< 0, -1 >()])*ihy);
-//
-//
-// if(xb - xf > 0.0)
-// xf = 0.0;
-// else
-// xb = 0.0;
-//
-// if(yb - yf > 0.0)
-// yf = 0.0;
-// 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
-// {
-// return 0.0;
-// }
+
+ return signui*c;
}
--
GitLab