From 131f0f6f7217513068367fb1305c655c3195a45d Mon Sep 17 00:00:00 2001
From: Tomas Sobotik <sobotto4@fjfi.cvut.cz>
Date: Wed, 4 Mar 2015 08:15:40 +0100
Subject: [PATCH] separation of synchronization cycles and few other changes
---
.../tnlParallelEikonalSolver.h | 2 +-
.../tnlParallelEikonalSolver_impl.h | 399 ++++++++++--------
.../parallelGodunovEikonal2D_impl.h | 19 +-
3 files changed, 249 insertions(+), 171 deletions(-)
diff --git a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver.h b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver.h
index 58babd8ce2..992335a929 100644
--- a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver.h
+++ b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver.h
@@ -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;
diff --git a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h
index fb1c810503..07435e3629 100644
--- a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h
+++ b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h
@@ -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,75 +186,83 @@ void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::synchronize() //nee
int tmp1, tmp2;
int grid1, grid2;
- for(int j = 0; j < this->gridRows - 1; j++)
+ if(this->currentStep & 1)
{
- for (int i = 0; i < this->gridCols*this->n; i++)
+ for(int j = 0; j < this->gridRows - 1; j++)
{
- 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 ((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 i = 0; i < this->gridCols*this->n; i++)
{
- this->work_u[tmp2] = this->work_u[tmp1];
- this->unusedCell[tmp2] = 0;
- if(grid2 == INT_MAX)
+ 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))
{
- 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)) & 8) )
+ setBoundaryCondition(getOwner(tmp2), getBoundaryCondition(getOwner(tmp2))+8);
}
- 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];
+ this->unusedCell[tmp1] = 0;
+ if(grid1 == INT_MAX)
+ {
+ setSubgridValue(getOwner(tmp1), -INT_MAX);
+ }
+ if(! (getBoundaryCondition(getOwner(tmp1)) & 1) )
+ setBoundaryCondition(getOwner(tmp1), getBoundaryCondition(getOwner(tmp1))+1);
}
- if(! (getBoundaryCondition(getOwner(tmp1)) & 1) )
- setBoundaryCondition(getOwner(tmp1), getBoundaryCondition(getOwner(tmp1))+1);
}
}
}
-
-
- for(int i = 1; i < this->gridCols; i++)
+ else
{
- for (int j = 0; j < this->gridRows*this->n; j++)
+ for(int i = 1; i < this->gridCols; i++)
{
- 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 ((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 j = 0; j < this->gridRows*this->n; j++)
{
- 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;
+ 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))
{
- 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);
}
- 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)
+ 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)) & 2) )
+ setBoundaryCondition(getOwner(tmp1), getBoundaryCondition(getOwner(tmp1))+2);
}
- if(! (getBoundaryCondition(getOwner(tmp1)) & 2) )
- 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
*/
+ /**/
+
+ /*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.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();
+ }
+ }
+
+
+ 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();
+ }
+ }
+
+
+ else if(y == 0 && (boundaryCondition & 8) && x ==0)
+ {
+ if((u[subMesh.getCellXSuccessor( i )] - u[i])/subMesh.getHx() > 1.0)
+ {
+ //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();
+ }
+ }
+
+
+ else if(y == subMesh.getDimensions().y() - 1 && (boundaryCondition & 1) && x ==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)
+ {
+ 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();
+ }
+ }
+ }*/
+
+ /**/
+
+
+ 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)
+ {
+ 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 (u.max() > 0.0)
+ this->stopTime *=(double) this->gridCols;*/
+
+
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 );
+ double finalTime = this->stopTime;// + 3.0*(u.max() - u.min());
+ if( time + currentTau > finalTime ) currentTau = finalTime - time;
- this -> refreshSolverMonitor();*/
-
- /****
- * Start the main loop
- */
- //cout << time << flush;
double maxResidue( 0.0 );
- while( time < this->stopTime/* || maxResidue > (10.0* this -> cflCondition)*/)
+ while( time < finalTime)
{
/****
* Compute the RHS
*/
- //this->problem->getExplicitRHS( time, currentTau, u, fu );
for( int i = 0; i < fu.getSize(); i ++ )
{
fu[ i ] = scheme.getValue( this->subMesh, i, this->subMesh.getCellCoordinates(i), u, time, boundaryCondition );
}
+ maxResidue = fu. absMax();
- //RealType lastResidue = this->getResidue();
+ if( this -> cflCondition * maxResidue != 0.0)
+ currentTau = this -> cflCondition / maxResidue;
- 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 )
- {
- maxResidue = fu. absMax();
- if( currentTau * maxResidue > this -> cflCondition )
- {
- //cout << '\r' << "For " << currentTau * maxResidue << " : " << currentTau << " --> ";
- currentTau = 0.95 * this -> cflCondition / maxResidue;
- //cout << currentTau << " at " << this -> cflCondition << flush;
- //currentTau *= 0.9;
-// continue;
- }
- }*/
-
- //RealType newResidue( 0.0 );
- //computeNewTimeLevel( u, currentTau, newResidue );
-
- //skip = false;
- if( time + currentTau > this -> stopTime )
- currentTau = this -> stopTime - time;
+ if( time + currentTau > finalTime ) currentTau = finalTime - time;
for( int i = 0; i < fu.getSize(); i ++ )
{
- if((u[i]+currentTau * fu[ i ])*Sign(u[i]) < 0.0 && fu[i] != 0.0 && u[i] != 0.0 )
- {
-// skip = true;
+ 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]));
- //currentTau *= 0.9;
-// i=fu.getSize();
- }
}
- // if(skip)
- //{
- // continue;
- //}
for( int i = 0; i < fu.getSize(); i ++ )
{
- //double add = currentTau * fu[ i ];
- u[ i ] += currentTau * fu[ i ];
- //localResidue += fabs( add );
+ double add = u[i] + currentTau * fu[ i ];
+ //if( fabs(u[i]) < fabs(add) or (this->subgridValues[subGridID] == this->currentStep +4) )
+ u[ i ] = add;
}
-
- //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( ! this->nextIteration() )
- return false;*/
-
- /****
- * Compute the new time step.
- */
- /*if( this -> cflCondition != 0.0 )
- currentTau *= 1.1;*/
-
-// if( time + currentTau > this -> stopTime )
-// {
-// currentTau = this -> stopTime - time; //we don't want to keep such tau
-// }
- //else this -> tau = currentTau;
-
- //this -> refreshSolverMonitor();
-
- /****
- * Check stop conditions.
- */
- /*if( time >= this->getStopTime() ||
- ( this -> getConvergenceResidue() != 0.0 && this->getResidue() < this -> getConvergenceResidue() ) )
- {
- //this -> refreshSolverMonitor();
- return true;
- }*/
- // cout << '\r' << flush;
- // cout << maxResidue << flush;
+ //cout << '\r' << 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;
}
diff --git a/src/implementation/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h b/src/implementation/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h
index 22078dd728..489e0a2ad1 100644
--- a/src/implementation/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h
+++ b/src/implementation/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h
@@ -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)) )
+ (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
--
GitLab