diff --git a/examples/heat-equation/heatEquationSolver_impl.h b/examples/heat-equation/heatEquationSolver_impl.h
index 377d4a92b52f21e1b475ece03fcac6b6c6bb4982..de42684c1c8deeb7a20c322e9327c87ba5f760bd 100644
--- a/examples/heat-equation/heatEquationSolver_impl.h
+++ b/examples/heat-equation/heatEquationSolver_impl.h
@@ -125,8 +125,8 @@ bindAuxiliaryDofs( const MeshType& mesh,
template< typename Mesh, typename DifferentialOperator, typename BoundaryCondition, typename RightHandSide >
-bool heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide >
-:: setInitialCondition( const tnlParameterContainer& parameters,
+bool heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide >::
+setInitialCondition( const tnlParameterContainer& parameters,
const MeshType& mesh,
DofVectorType& dofs )
{
@@ -149,13 +149,19 @@ heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide
setupLinearSystem( const MeshType& mesh,
MatrixType& matrix )
{
+ const IndexType dofs = this->getDofs( mesh );
RowLengthsVectorType rowLengths;
+ if( ! rowLengths.setSize( dofs ) )
+ return false;
tnlMatrixSetter< MeshType, DifferentialOperator, BoundaryCondition, RowLengthsVectorType > matrixSetter;
matrixSetter.template getRowLengths< Mesh::Dimensions >( mesh,
differentialOperator,
boundaryCondition,
rowLengths );
- matrix.setRowLengths( rowLengths );
+ matrix.setDimensions( dofs, dofs );
+ if( ! matrix.setRowLengths( rowLengths ) )
+ return false;
+ return true;
}
template< typename Mesh,
@@ -226,7 +232,7 @@ assemblyLinearSystem( const RealType& time,
DofVectorType& b )
{
tnlLinearSystemAssembler< Mesh, DofVectorType, DifferentialOperator, BoundaryCondition, RightHandSide, MatrixType > systemAssembler;
- /*systemAssembler.template assembly< Mesh::Dimensions >( time,
+ systemAssembler.template assembly< Mesh::Dimensions >( time,
tau,
mesh,
this->differentialOperator,
@@ -235,7 +241,8 @@ assemblyLinearSystem( const RealType& time,
u,
matrix,
b );
- */
+ //matrix.print( cout );
+ //abort();
}
template< typename Mesh,
diff --git a/examples/heat-equation/tnlHeatEquationEocRhs.h b/examples/heat-equation/tnlHeatEquationEocRhs.h
index 73a6108284039a8219663c4f92fa271ebe065298..95bd607cc043cdddd129c80a3f1e6d4674f3b7fa 100644
--- a/examples/heat-equation/tnlHeatEquationEocRhs.h
+++ b/examples/heat-equation/tnlHeatEquationEocRhs.h
@@ -37,7 +37,7 @@ class tnlHeatEquationEocRhs
template< typename Vertex,
typename Real >
Real getValue( const Vertex& vertex,
- const Real& time )
+ const Real& time ) const
{
return testFunction.getTimeDerivative( vertex, time )
- exactOperator.getValue( testFunction, vertex, time );
diff --git a/src/config/tnlParameterContainer.h b/src/config/tnlParameterContainer.h
index 606350cf080828ed75bb3227aecde53424279ef5..cc56df2d58601f186fa6866f367996101889f9fd 100644
--- a/src/config/tnlParameterContainer.h
+++ b/src/config/tnlParameterContainer.h
@@ -89,7 +89,8 @@ class tnlParameterContainer
for( i = 0; i < size; i ++ )
if( parameters[ i ] -> name == name )
return ( ( tnlParameter< T >* ) parameters[ i ] ) -> value;
- cerr << "Unknown parameter " << name << endl;
+ cerr << "The program attempts to get unknown parameter " << name << endl;
+ cerr << "Aborting the program." << endl;
abort();
}
diff --git a/src/implementation/matrices/tnlCSRMatrix_impl.h b/src/implementation/matrices/tnlCSRMatrix_impl.h
index b48a238f2b694fbc86708841fbd853b520c109e3..8ae0f660fc1bd137a947dfd8c7cab2b1397c9794 100644
--- a/src/implementation/matrices/tnlCSRMatrix_impl.h
+++ b/src/implementation/matrices/tnlCSRMatrix_impl.h
@@ -77,11 +77,14 @@ bool tnlCSRMatrix< Real, Device, Index >::setRowLengths( const RowLengthsVector&
* necessary length of the vectors this->values
* and this->columnIndexes.
*/
+ tnlAssert( this->getRows() > 0, );
+ tnlAssert( this->getColumns() > 0, );
tnlSharedVector< IndexType, DeviceType, IndexType > rowPtrs;
rowPtrs.bind( this->rowPointers.getData(), this->getRows() );
rowPtrs = rowLengths;
this->rowPointers.setElement( this->rows, 0 );
this->rowPointers.computeExclusivePrefixSum();
+ this->maxRowLength = rowLengths.max();
/****
* Allocate values and column indexes
diff --git a/src/implementation/matrices/tnlChunkedEllpackMatrix_impl.h b/src/implementation/matrices/tnlChunkedEllpackMatrix_impl.h
index 9affb69e2f4c23e10254529b82c692b50d6a9138..4ca1fc091773063d8ec7b3e338ea8d97b3d5d393 100644
--- a/src/implementation/matrices/tnlChunkedEllpackMatrix_impl.h
+++ b/src/implementation/matrices/tnlChunkedEllpackMatrix_impl.h
@@ -219,6 +219,8 @@ template< typename Real,
typename Index >
bool tnlChunkedEllpackMatrix< Real, Device, Index >::setRowLengths( const RowLengthsVector& rowLengths )
{
+ tnlAssert( this->getRows() > 0, );
+ tnlAssert( this->getColumns() > 0, );
IndexType elementsToAllocation( 0 );
@@ -252,6 +254,7 @@ bool tnlChunkedEllpackMatrix< Real, Device, Index >::setRowLengths( const RowLen
this->numberOfSlices = hostMatrix.numberOfSlices;
elementsToAllocation = hostMatrix.values.getSize();
}
+ this->maxRowLength = rowLengths.max();
return tnlSparseMatrix< Real, Device, Index >::allocateMatrixElements( elementsToAllocation );
}
diff --git a/src/implementation/matrices/tnlEllpackMatrix_impl.h b/src/implementation/matrices/tnlEllpackMatrix_impl.h
index 1257741aefed7370877eaebc851251a108e563dd..d82d8704dbbd320687c7532bac02e518fd957e6c 100644
--- a/src/implementation/matrices/tnlEllpackMatrix_impl.h
+++ b/src/implementation/matrices/tnlEllpackMatrix_impl.h
@@ -74,12 +74,12 @@ template< typename Real,
typename Index >
bool tnlEllpackMatrix< Real, Device, Index >::setRowLengths( const RowLengthsVector& rowLengths )
{
- this->rowLengths = rowLengths.max();
+ tnlAssert( this->getRows() > 0, );
+ tnlAssert( this->getColumns() > 0, );
tnlAssert( this->rowLengths > 0,
cerr << "this->rowLengths = " << this->rowLengths );
- if( this->rows > 0 )
- return allocateElements();
- return true;
+ this->rowLengths = this->maxRowLength = rowLengths.max();
+ return allocateElements();
}
template< typename Real,
diff --git a/src/implementation/matrices/tnlMatrixSetter_impl.h b/src/implementation/matrices/tnlMatrixSetter_impl.h
index 97a575ab7ab7cc8f895bd269cd4a14ddc3b43cce..b8a808013f2b9879d5c86b46e4a307fa7cab5462 100644
--- a/src/implementation/matrices/tnlMatrixSetter_impl.h
+++ b/src/implementation/matrices/tnlMatrixSetter_impl.h
@@ -62,13 +62,13 @@ getRowLengths( const MeshType& mesh,
TraversalBoundaryEntitiesProcessor boundaryEntitiesProcessor;
TraversalInteriorEntitiesProcessor interiorEntitiesProcessor;
tnlTraversal< MeshType, EntityDimensions > meshTraversal;
- /*meshTraversal.template processEntities< TraversalUserData,
+ meshTraversal.template processEntities< TraversalUserData,
TraversalBoundaryEntitiesProcessor,
TraversalInteriorEntitiesProcessor >
( mesh,
userData,
boundaryEntitiesProcessor,
- interiorEntitiesProcessor );*/
+ interiorEntitiesProcessor );
}
diff --git a/src/implementation/matrices/tnlMultidiagonalMatrix_impl.h b/src/implementation/matrices/tnlMultidiagonalMatrix_impl.h
index f24c3935a991b5b9ac8d1717f2099c7d65e667b7..0fd8bef475d00041df399bbe53857d618895d47e 100644
--- a/src/implementation/matrices/tnlMultidiagonalMatrix_impl.h
+++ b/src/implementation/matrices/tnlMultidiagonalMatrix_impl.h
@@ -160,6 +160,16 @@ Index tnlMultidiagonalMatrix< Real, Device, Index > :: getNumberOfNonzeroMatrixE
return nonzeroElements;
}
+template< typename Real,
+ typename Device,
+ typename Index >
+Index
+tnlMultidiagonalMatrix< Real, Device, Index >::
+getMaxRowlength() const
+{
+ return diagonalsShift.getSize();
+}
+
template< typename Real,
typename Device,
typename Index >
diff --git a/src/implementation/matrices/tnlSlicedEllpackMatrix_impl.h b/src/implementation/matrices/tnlSlicedEllpackMatrix_impl.h
index 203087460d471eb373e282f70bb5a1e031fb912a..f2579ee6b6f4781cc8759a1b6e6751be4dc9ffa7 100644
--- a/src/implementation/matrices/tnlSlicedEllpackMatrix_impl.h
+++ b/src/implementation/matrices/tnlSlicedEllpackMatrix_impl.h
@@ -71,6 +71,8 @@ template< typename Real,
int SliceSize >
bool tnlSlicedEllpackMatrix< Real, Device, Index, SliceSize >::setRowLengths( const RowLengthsVector& rowLengths )
{
+ tnlAssert( this->getRows() > 0, );
+ tnlAssert( this->getColumns() > 0, );
const IndexType slices = roundUpDivision( this->rows, SliceSize );
if( ! this->sliceRowLengths.setSize( slices ) ||
! this->slicePointers.setSize( slices + 1 ) )
@@ -78,8 +80,9 @@ bool tnlSlicedEllpackMatrix< Real, Device, Index, SliceSize >::setRowLengths( co
DeviceDependentCode::computeMaximalRowLengthInSlices( *this, rowLengths );
- this->slicePointers.computeExclusivePrefixSum();
+ this->maxRowLength = rowLengths.max();
+ this->slicePointers.computeExclusivePrefixSum();
return this->allocateMatrixElements( this->slicePointers.getElement( slices ) );
}
diff --git a/src/implementation/matrices/tnlSparseMatrix_impl.h b/src/implementation/matrices/tnlSparseMatrix_impl.h
index 399b6584024174e77771d9b1945eae8a793e19e6..6ec90d0e9414dd62b33c7f3eacf28a75836bdd87 100644
--- a/src/implementation/matrices/tnlSparseMatrix_impl.h
+++ b/src/implementation/matrices/tnlSparseMatrix_impl.h
@@ -22,6 +22,7 @@ template< typename Real,
typename Device,
typename Index >
tnlSparseMatrix< Real, Device, Index >::tnlSparseMatrix()
+: maxRowLength( 0 )
{
}
@@ -60,6 +61,16 @@ Index tnlSparseMatrix< Real, Device, Index >::getNumberOfNonzeroMatrixElements()
return nonzeroElements;
}
+template< typename Real,
+ typename Device,
+ typename Index >
+Index
+tnlSparseMatrix< Real, Device, Index >::
+getMaxRowLength() const
+{
+ return this->maxRowLength;
+}
+
template< typename Real,
typename Device,
typename Index >
diff --git a/src/implementation/matrices/tnlTridiagonalMatrix_impl.h b/src/implementation/matrices/tnlTridiagonalMatrix_impl.h
index 363391a28263ad7c74a6a3931c75f46da09b42b8..c46a8fdb8cbcc552937197862ddce26eba352fd0 100644
--- a/src/implementation/matrices/tnlTridiagonalMatrix_impl.h
+++ b/src/implementation/matrices/tnlTridiagonalMatrix_impl.h
@@ -133,6 +133,16 @@ Index tnlTridiagonalMatrix< Real, Device, Index > :: getNumberOfNonzeroMatrixEle
return nonzeroElements;
}
+template< typename Real,
+ typename Device,
+ typename Index >
+Index
+tnlTridiagonalMatrix< Real, Device, Index >::
+getMaxRowlength() const
+{
+ return 3;
+}
+
template< typename Real,
typename Device,
typename Index >
diff --git a/src/implementation/operators/diffusion/tnlLinearDiffusion_impl.h b/src/implementation/operators/diffusion/tnlLinearDiffusion_impl.h
index 28e15c000e67357e8949df67709826520d07bc51..2acc343c26ef8f19fb69487ac86ef157a3e0a1b1 100644
--- a/src/implementation/operators/diffusion/tnlLinearDiffusion_impl.h
+++ b/src/implementation/operators/diffusion/tnlLinearDiffusion_impl.h
@@ -59,6 +59,37 @@ getLinearSystemRowLength( const MeshType& mesh,
return 3;
}
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+ template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlLinearDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
+updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ IndexType* columns,
+ RealType* values,
+ IndexType& rowLength ) const
+{
+ columns[ 0 ] = mesh.getCellXPredecessor( index );
+ columns[ 1 ] = index;
+ columns[ 2 ] = mesh.getCellXSuccessor( index );
+ values[ 0 ] = -tau;
+ values[ 1 ] = 1.0 + 2.0 * tau;
+ values[ 2 ] = -tau;
+ rowLength = 3;
+}
+
template< typename MeshReal,
typename Device,
typename MeshIndex,
@@ -117,6 +148,42 @@ getValue( const MeshType& mesh,
+ u[ mesh.getCellYSuccessor( cellIndex ) ] ) * mesh.getHySquareInverse();
}
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+ template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ IndexType* columns,
+ RealType* values,
+ IndexType& rowLength ) const
+{
+ columns[ 0 ] = mesh.getCellYPredecessor( index );
+ columns[ 1 ] = mesh.getCellXPredecessor( index );
+ columns[ 2 ] = index;
+ columns[ 3 ] = mesh.getCellXSuccessor( index );
+ columns[ 4 ] = mesh.getCellYSuccessor( index );
+ values[ 0 ] = -tau;
+ values[ 1 ] = -tau;
+ values[ 2 ] = 1.0 + 4.0 * tau;
+ values[ 3 ] = -tau;
+ values[ 4 ] = -tau;
+ rowLength = 5;
+}
+
+
template< typename MeshReal,
typename Device,
typename MeshIndex,
@@ -177,5 +244,45 @@ getLinearSystemRowLength( const MeshType& mesh,
return 7;
}
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+ template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
+updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ IndexType* columns,
+ RealType* values,
+ IndexType& rowLength ) const
+{
+ columns[ 0 ] = mesh.getCellZPredecessor( index );
+ columns[ 1 ] = mesh.getCellYPredecessor( index );
+ columns[ 2 ] = mesh.getCellXPredecessor( index );
+ columns[ 3 ] = index;
+ columns[ 4 ] = mesh.getCellXSuccessor( index );
+ columns[ 5 ] = mesh.getCellYSuccessor( index );
+ columns[ 6 ] = mesh.getCellZSuccessor( index );
+ values[ 0 ] = -tau;
+ values[ 1 ] = -tau;
+ values[ 2 ] = -tau;
+ values[ 3 ] = 1.0 + 6.0 * tau;
+ values[ 4 ] = -tau;
+ values[ 5 ] = -tau;
+ values[ 6 ] = -tau;
+ rowLength = 7;
+}
+
+
#endif /* TNLLINEARDIFFUSION_IMP_H */
diff --git a/src/implementation/operators/tnlDirichletBoundaryConditions_impl.h b/src/implementation/operators/tnlDirichletBoundaryConditions_impl.h
index ce99c3b4837465bd01ff2b9e9834667509f4af85..6c781e00db627ea43fab0379e82bf9d2f8ff7cda 100644
--- a/src/implementation/operators/tnlDirichletBoundaryConditions_impl.h
+++ b/src/implementation/operators/tnlDirichletBoundaryConditions_impl.h
@@ -35,7 +35,7 @@ setBoundaryConditions( const RealType& time,
DofVectorType& fu )
{
fu[ index ] = 0;
- u[ index ] = function.getValue( mesh.getVertex( coordinates ), time );
+ u[ index ] = function.getValue( mesh.getCellCenter( coordinates ), time );
}
template< typename MeshReal,
@@ -56,6 +56,33 @@ getLinearSystemRowLength( const MeshType& mesh,
return 1;
}
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Function,
+ typename Real,
+ typename Index >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlDirichletBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >::
+updateLinearSystem( const RealType& time,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ DofVectorType& u,
+ DofVectorType& b,
+ IndexType* columns,
+ RealType* values,
+ IndexType& rowLength ) const
+{
+ columns[ 0 ] = index;
+ values[ 0 ] = 1.0;
+ b[ index ] = function.getValue( mesh.getCellCenter( coordinates ), time );
+ rowLength = 1;
+}
+
template< typename MeshReal,
typename Device,
@@ -91,7 +118,7 @@ setBoundaryConditions( const RealType& time,
DofVectorType& fu )
{
fu[ index ] = 0;
- u[ index ] = function.getValue( mesh.getVertex( coordinates ), time );;
+ u[ index ] = function.getValue( mesh.getCellCenter( coordinates ), time );;
}
template< typename MeshReal,
@@ -112,6 +139,34 @@ getLinearSystemRowLength( const MeshType& mesh,
return 1;
}
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Function,
+ typename Real,
+ typename Index >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlDirichletBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >::
+updateLinearSystem( const RealType& time,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ DofVectorType& u,
+ DofVectorType& b,
+ IndexType* columns,
+ RealType* values,
+ IndexType& rowLength ) const
+{
+ columns[ 0 ] = index;
+ values[ 0 ] = 1.0;
+ b[ index ] = function.getValue( mesh.getCellCenter( coordinates ), time );
+ rowLength = 1;
+}
+
+
template< typename MeshReal,
typename Device,
@@ -146,7 +201,7 @@ setBoundaryConditions( const RealType& time,
DofVectorType& fu )
{
fu[ index ] = 0;
- u[ index ] = function.getValue( mesh.getVertex( coordinates ), time );;
+ u[ index ] = function.getValue( mesh.getCellCenter( coordinates ), time );;
}
template< typename MeshReal,
@@ -167,5 +222,33 @@ getLinearSystemRowLength( const MeshType& mesh,
return 1;
}
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Function,
+ typename Real,
+ typename Index >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlDirichletBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >::
+updateLinearSystem( const RealType& time,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ DofVectorType& u,
+ DofVectorType& b,
+ IndexType* columns,
+ RealType* values,
+ IndexType& rowLength ) const
+{
+ columns[ 0 ] = index;
+ values[ 0 ] = 1.0;
+ b[ index ] = function.getValue( mesh.getCellCenter( coordinates ), time );
+ rowLength = 1;
+}
+
+
#endif /* TNLDIRICHLETBOUNDARYCONDITIONS_IMPL_H */
diff --git a/src/implementation/solvers/linear/krylov/tnlBICGStabSolver_impl.h b/src/implementation/solvers/linear/krylov/tnlBICGStabSolver_impl.h
index fe7d92a2d585e7197a7c0ebdc03e851c77f61574..18dee1fbb7c681409fd59e985d038f8d3ee4e2eb 100644
--- a/src/implementation/solvers/linear/krylov/tnlBICGStabSolver_impl.h
+++ b/src/implementation/solvers/linear/krylov/tnlBICGStabSolver_impl.h
@@ -43,6 +43,26 @@ tnlString tnlBICGStabSolver< Matrix, Preconditioner > :: getType() const
this -> preconditioner -> getType() + " >";
}
+template< typename Matrix,
+ typename Preconditioner >
+void
+tnlBICGStabSolver< Matrix, Preconditioner >::
+configSetup( tnlConfigDescription& config,
+ const tnlString& prefix )
+{
+ tnlIterativeSolver< RealType, IndexType >::configSetup( config, prefix );
+}
+
+template< typename Matrix,
+ typename Preconditioner >
+bool
+tnlBICGStabSolver< Matrix, Preconditioner >::
+setup( const tnlParameterContainer& parameters,
+ const tnlString& prefix )
+{
+ tnlIterativeSolver< RealType, IndexType >::setup( parameters, prefix );
+}
+
template< typename Matrix,
typename Preconditioner >
void tnlBICGStabSolver< Matrix, Preconditioner > :: setMatrix( const MatrixType& matrix )
diff --git a/src/implementation/solvers/linear/krylov/tnlCGSolver_impl.h b/src/implementation/solvers/linear/krylov/tnlCGSolver_impl.h
index fd1bb90c1c68faad9de4e2f2b53e032366f149c3..22d5096f53ff7f44c8fc9da6f0f26896da0dd1ce 100644
--- a/src/implementation/solvers/linear/krylov/tnlCGSolver_impl.h
+++ b/src/implementation/solvers/linear/krylov/tnlCGSolver_impl.h
@@ -33,6 +33,26 @@ tnlString tnlCGSolver< Matrix, Preconditioner > :: getType() const
this -> preconditioner -> getType() + " >";
}
+template< typename Matrix,
+ typename Preconditioner >
+void
+tnlCGSolver< Matrix, Preconditioner >::
+configSetup( tnlConfigDescription& config,
+ const tnlString& prefix )
+{
+ tnlIterativeSolver< RealType, IndexType >::configSetup( config, prefix );
+}
+
+template< typename Matrix,
+ typename Preconditioner >
+bool
+tnlCGSolver< Matrix, Preconditioner >::
+setup( const tnlParameterContainer& parameters,
+ const tnlString& prefix )
+{
+ tnlIterativeSolver< RealType, IndexType >::setup( parameters, prefix );
+}
+
template< typename Matrix,
typename Preconditioner >
void tnlCGSolver< Matrix, Preconditioner > :: setMatrix( const MatrixType& matrix )
diff --git a/src/implementation/solvers/linear/krylov/tnlGMRESSolver_impl.h b/src/implementation/solvers/linear/krylov/tnlGMRESSolver_impl.h
index 4d697d0887f85004cdbe054d83d60991b5680f45..e115dace6a495b85f4b1a77dd627d8ef1b30cfca 100644
--- a/src/implementation/solvers/linear/krylov/tnlGMRESSolver_impl.h
+++ b/src/implementation/solvers/linear/krylov/tnlGMRESSolver_impl.h
@@ -38,14 +38,38 @@ tnlGMRESSolver< Matrix, Preconditioner > :: tnlGMRESSolver()
};
template< typename Matrix,
- typename Preconditioner >
-tnlString tnlGMRESSolver< Matrix, Preconditioner > :: getType() const
+ typename Preconditioner >
+tnlString
+tnlGMRESSolver< Matrix, Preconditioner >::
+getType() const
{
return tnlString( "tnlGMRESSolver< " ) +
this -> matrix -> getType() + ", " +
this -> preconditioner -> getType() + " >";
}
+template< typename Matrix,
+ typename Preconditioner >
+void
+tnlGMRESSolver< Matrix, Preconditioner >::
+configSetup( tnlConfigDescription& config,
+ const tnlString& prefix )
+{
+ tnlIterativeSolver< RealType, IndexType >::configSetup( config, prefix );
+ config.addEntry< int >( prefix + "gmres-restarting", "Number of iterations after which the GMRES restarts.", 10 );
+}
+
+template< typename Matrix,
+ typename Preconditioner >
+bool
+tnlGMRESSolver< Matrix, Preconditioner >::
+setup( const tnlParameterContainer& parameters,
+ const tnlString& prefix )
+{
+ tnlIterativeSolver< RealType, IndexType >::setup( parameters, prefix );
+ this->setRestarting( parameters.GetParameter< int >( "gmres-restarting" ) );
+}
+
template< typename Matrix,
typename Preconditioner >
void tnlGMRESSolver< Matrix, Preconditioner > :: setRestarting( IndexType rest )
@@ -155,7 +179,9 @@ bool tnlGMRESSolver< Matrix, Preconditioner > :: solve( const Vector& b, Vector&
- //dbgCout( " ----------- Starting m-loop -----------------" );
+ /****
+ * Starting m-loop
+ */
for( i = 0; i < m && this -> getIterations() <= this -> getMaxIterations(); i++ )
{
vi. bind( &( _v. getData()[ i * size ] ), size );
@@ -196,8 +222,9 @@ bool tnlGMRESSolver< Matrix, Preconditioner > :: solve( const Vector& b, Vector&
vi. bind( &( _v. getData()[ ( i + 1 ) * size ] ), size );
vi. addVector( _w, ( RealType ) 1.0 / normw );
-
- //dbgCout( "Applying rotations" );
+ /****
+ * Applying the Givens rotations
+ */
for( k = 0; k < i; k++ )
applyPlaneRotation( H[ k + i * ( m + 1 )],
H[ k + 1 + i * ( m + 1 ) ],
@@ -217,7 +244,8 @@ bool tnlGMRESSolver< Matrix, Preconditioner > :: solve( const Vector& b, Vector&
cs[ i ],
sn[ i ] );
- this -> setResidue( fabs( s[ i + 1 ] ) / normb );
+ this->setResidue( fabs( s[ i + 1 ] ) / normb );
+ this->refreshSolverMonitor();
if( this -> getResidue() < this -> getConvergenceResidue() )
{
@@ -226,15 +254,14 @@ bool tnlGMRESSolver< Matrix, Preconditioner > :: solve( const Vector& b, Vector&
// this -> printOut();
return true;
}
- //DBG_WAIT;
- this -> refreshSolverMonitor();
if( ! this -> nextIteration() )
return false;
}
update( m - 1, m, _H, _s, _v, x );
- //dbgCout_ARRAY( x, size );
- // r = M.solve(b - A * x);
+ /****
+ * r = M.solve(b - A * x);
+ */
beta = 0.0;
if( preconditioner )
{
@@ -257,7 +284,7 @@ bool tnlGMRESSolver< Matrix, Preconditioner > :: solve( const Vector& b, Vector&
return false;
}
this -> refreshSolverMonitor();
- if( this -> getResidue() >this -> getConvergenceResidue() ) return false;
+ if( this -> getResidue() > this -> getConvergenceResidue() ) return false;
return true;
};
diff --git a/src/implementation/solvers/linear/krylov/tnlTFQMRSolver_impl.h b/src/implementation/solvers/linear/krylov/tnlTFQMRSolver_impl.h
index a9af14cf01a527a9272808975e7cac6fa07bf893..1c3dd5adf345fb8334fbf432a9555e5a4e80a2ed 100644
--- a/src/implementation/solvers/linear/krylov/tnlTFQMRSolver_impl.h
+++ b/src/implementation/solvers/linear/krylov/tnlTFQMRSolver_impl.h
@@ -42,6 +42,26 @@ tnlString tnlTFQMRSolver< Matrix, Preconditioner > :: getType() const
tnlString( getType( ( IndexType ) 0 ) ) + " >";
}
+template< typename Matrix,
+ typename Preconditioner >
+void
+tnlTFQMRSolver< Matrix, Preconditioner >::
+configSetup( tnlConfigDescription& config,
+ const tnlString& prefix )
+{
+ tnlIterativeSolver< RealType, IndexType >::configSetup( config, prefix );
+}
+
+template< typename Matrix,
+ typename Preconditioner >
+bool
+tnlTFQMRSolver< Matrix, Preconditioner >::
+setup( const tnlParameterContainer& parameters,
+ const tnlString& prefix )
+{
+ tnlIterativeSolver< RealType, IndexType >::setup( parameters, prefix );
+}
+
template< typename Matrix,
typename Preconditioner >
void tnlTFQMRSolver< Matrix, Preconditioner > :: setMatrix( const MatrixType& matrix )
diff --git a/src/implementation/solvers/linear/stationary/tnlSORSolver_impl.h b/src/implementation/solvers/linear/stationary/tnlSORSolver_impl.h
index 4e90e3a0c5954a7501a82eec21562f794253e55c..54c121d9eb3dbf12e341dd39abe32f1d13c913db 100644
--- a/src/implementation/solvers/linear/stationary/tnlSORSolver_impl.h
+++ b/src/implementation/solvers/linear/stationary/tnlSORSolver_impl.h
@@ -32,6 +32,29 @@ tnlString tnlSORSolver< Matrix, Preconditioner > :: getType() const
this -> preconditioner -> getType() + " >";
}
+template< typename Matrix,
+ typename Preconditioner >
+void
+tnlSORSolver< Matrix, Preconditioner >::
+configSetup( tnlConfigDescription& config,
+ const tnlString& prefix )
+{
+ tnlIterativeSolver< RealType, IndexType >::configSetup( config, prefix );
+ config.addEntry< double >( prefix + "sor-omega", "Relaxation parameter of the SOR method.", 1.0 );
+}
+
+template< typename Matrix,
+ typename Preconditioner >
+bool
+tnlSORSolver< Matrix, Preconditioner >::
+setup( const tnlParameterContainer& parameters,
+ const tnlString& prefix )
+{
+ tnlIterativeSolver< RealType, IndexType >::setup( parameters, prefix );
+ this->setOmega( parameters.GetParameter< int >( "sor-omega" ) );
+}
+
+
template< typename Matrix, typename Preconditioner >
void tnlSORSolver< Matrix, Preconditioner > :: setOmega( const RealType& omega )
{
@@ -67,12 +90,12 @@ bool tnlSORSolver< Matrix, Preconditioner > :: solve( const Vector& b,
const IndexType size = matrix -> getRows();
this -> resetIterations();
- this -> setResidue( this -> getMaxResidue() + 1.0 );
+ this -> setResidue( this -> getConvergenceResidue() + 1.0 );
RealType bNorm = b. lpNorm( ( RealType ) 2.0 );
while( this -> getIterations() < this -> getMaxIterations() &&
- this -> getResidue() > this -> getMaxResidue() )
+ this -> getResidue() > this -> getConvergenceResidue() )
{
/*matrix -> performSORIteration( this -> getOmega(),
b,
@@ -87,7 +110,7 @@ bool tnlSORSolver< Matrix, Preconditioner > :: solve( const Vector& b,
}
this -> setResidue( ResidueGetter :: getResidue( *matrix, b, x, bNorm ) );
this -> refreshSolverMonitor();
- if( this -> getResidue() > this -> getMaxResidue() ) return false;
+ if( this -> getResidue() > this -> getConvergenceResidue() ) return false;
return true;
};
diff --git a/src/implementation/solvers/pde/CMakeLists.txt b/src/implementation/solvers/pde/CMakeLists.txt
index bebcfbe999c8a0c24264a1fb2fecef7692c538c2..69ce0ff2e828093b24ea0e0f3d8c406f0a67799c 100755
--- a/src/implementation/solvers/pde/CMakeLists.txt
+++ b/src/implementation/solvers/pde/CMakeLists.txt
@@ -1,7 +1,8 @@
SET( headers tnlExplicitTimeStepper_impl.h
tnlPDESolver_impl.h
tnlExplicitUpdater_impl.h
- tnlLinearSystemAssembler_impl.h )
+ tnlLinearSystemAssembler_impl.h
+ tnlSemiImplicitTimeStepper_impl.h )
INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/implementation/solvers/pde )
diff --git a/src/implementation/solvers/pde/tnlExplicitTimeStepper_impl.h b/src/implementation/solvers/pde/tnlExplicitTimeStepper_impl.h
index 68c671ba0451b0b5ce6e698951ea49183bc0f932..d8175fe3d3820104b2d47fdd2de7b118190a36a3 100644
--- a/src/implementation/solvers/pde/tnlExplicitTimeStepper_impl.h
+++ b/src/implementation/solvers/pde/tnlExplicitTimeStepper_impl.h
@@ -96,12 +96,13 @@ bool tnlExplicitTimeStepper< Problem, OdeSolver >::solve( const RealType& time,
const MeshType& mesh,
DofVectorType& dofVector )
{
+ tnlAssert( this->odeSolver, );
this->odeSolver->setTau( this -> timeStep );
this->odeSolver->setProblem( * this );
this->odeSolver->setTime( time );
this->odeSolver->setStopTime( stopTime );
this->mesh = &mesh;
- return this -> odeSolver -> solve( dofVector );
+ return this->odeSolver->solve( dofVector );
}
template< typename Problem,
diff --git a/src/implementation/solvers/pde/tnlLinearSystemAssembler_impl.h b/src/implementation/solvers/pde/tnlLinearSystemAssembler_impl.h
index f7f18c2259b36e11109b19a034ce7879c69566c5..0d366139b3aa4c2d0505e2ad5157c33a4a828cc6 100644
--- a/src/implementation/solvers/pde/tnlLinearSystemAssembler_impl.h
+++ b/src/implementation/solvers/pde/tnlLinearSystemAssembler_impl.h
@@ -76,7 +76,14 @@ assembly( const RealType& time,
MatrixType& matrix,
DofVector& b ) const
{
- TraversalUserData userData( time, tau, differentialOperator, boundaryConditions, rightHandSide, u, matrix, b );
+ const IndexType maxRowLength = matrix.getMaxRowLength();
+ tnlAssert( maxRowLength > 0, );
+ typename TraversalUserData::RowValuesType values;
+ typename TraversalUserData::RowColumnsType columns;
+ values.setSize( maxRowLength );
+ columns.setSize( maxRowLength );
+
+ TraversalUserData userData( time, tau, differentialOperator, boundaryConditions, rightHandSide, u, matrix, b, columns, values );
TraversalBoundaryEntitiesProcessor boundaryEntitiesProcessor;
TraversalInteriorEntitiesProcessor interiorEntitiesProcessor;
tnlTraversal< MeshType, EntityDimensions > meshTraversal;
diff --git a/src/implementation/solvers/pde/tnlSemiImplicitTimeStepper_impl.h b/src/implementation/solvers/pde/tnlSemiImplicitTimeStepper_impl.h
index 4bcfdc106798ee5803c4116e98018d5aeb1e46be..91f775f6e7f71be45a585d9daf64fe658e113171 100644
--- a/src/implementation/solvers/pde/tnlSemiImplicitTimeStepper_impl.h
+++ b/src/implementation/solvers/pde/tnlSemiImplicitTimeStepper_impl.h
@@ -37,6 +37,7 @@ tnlSemiImplicitTimeStepper< Problem, LinearSystemSolver >::
configSetup( tnlConfigDescription& config,
const tnlString& prefix )
{
+ config.addEntry< bool >( "verbose", "Verbose mode.", true );
}
template< typename Problem,
@@ -46,6 +47,7 @@ tnlSemiImplicitTimeStepper< Problem, LinearSystemSolver >::
setup( const tnlParameterContainer& parameters,
const tnlString& prefix )
{
+ this->verbose = parameters.GetParameter< bool >( "verbose" );
return true;
}
@@ -55,7 +57,21 @@ bool
tnlSemiImplicitTimeStepper< Problem, LinearSystemSolver >::
init( const MeshType& mesh )
{
- return this->problem->setupLinearSystem( mesh, this->matrix );
+ if( ! this->problem->setupLinearSystem( mesh, this->matrix ) )
+ return false;
+ if( this->matrix.getRows() == 0 || this->matrix.getColumns() == 0 )
+ {
+ cerr << "The matrix for the semi-implicit time stepping was not set correctly." << endl;
+ if( ! this->matrix.getRows() )
+ cerr << "The matrix dimensions are set to 0 rows." << endl;
+ if( ! this->matrix.getColumns() )
+ cerr << "The matrix dimensions are set to 0 columns." << endl;
+ cerr << "Please check the method 'setupLinearSystem' in your solver." << endl;
+ return false;
+ }
+ if( ! this->rightHandSide.setSize( this->matrix.getRows() ) )
+ return false;
+ return true;
}
template< typename Problem,
@@ -118,17 +134,29 @@ solve( const RealType& time,
{
tnlAssert( this->problem != 0, );
RealType t = time;
+ this->linearSystemSolver->setMatrix( this->matrix );
while( t < stopTime )
{
RealType currentTau = Min( this->timeStep, stopTime - t );
+
+ if( verbose )
+ cout << " Assembling the linear system ... \r" << flush;
this->problem->assemblyLinearSystem( t,
currentTau,
mesh,
dofVector,
this->matrix,
this->rightHandSide );
+ if( verbose )
+ cout << " Solving the linear system for time " << t << " \r" << flush;
+ if( ! this->linearSystemSolver->solve( this->rightHandSide, dofVector ) )
+ {
+ cerr << "The linear system solver did not converge." << endl;
+ return false;
+ }
t += currentTau;
}
+ return true;
}
#endif /* TNLSEMIIMPLICITTIMESTEPPER_IMPL_H_ */
diff --git a/src/implementation/solvers/tnlIterativeSolver_impl.h b/src/implementation/solvers/tnlIterativeSolver_impl.h
index 6583fa53bcd0a5f39d4f463946f8b6052013bc8c..e423b30e5af2a2bb5e5da47aa411f248acb2f966 100644
--- a/src/implementation/solvers/tnlIterativeSolver_impl.h
+++ b/src/implementation/solvers/tnlIterativeSolver_impl.h
@@ -89,10 +89,14 @@ void tnlIterativeSolver< Real, Index> :: resetIterations()
template< typename Real, typename Index >
bool tnlIterativeSolver< Real, Index> :: nextIteration()
{
+ // TODO: fix
+ //tnlAssert( solverMonitor, );
+ //solverMonitor->setIterations( ++this->currentIteration );
+ //solverMonitor->setResidue( this->getResidue() );
if( this->solverMonitor &&
this->currentIteration % this->refreshRate == 0 )
solverMonitor->refresh();
- this->currentIteration ++;
+
if( std::isnan( this->getResidue() ) ||
( this->getResidue() > this->getDivergenceResidue() &&
this->getIterations() > this->minIterations ) ||
diff --git a/src/implementation/solvers/tnlSolverConfig_impl.h b/src/implementation/solvers/tnlSolverConfig_impl.h
index 1cd3558fea3bc0fb9e847d112f77ff3be09f4b64..555385a39e732b0cfd6e6546bd45ac4166bc3df0 100644
--- a/src/implementation/solvers/tnlSolverConfig_impl.h
+++ b/src/implementation/solvers/tnlSolverConfig_impl.h
@@ -106,6 +106,17 @@ bool tnlSolverConfig< ConfigTag, ProblemConfig >::configSetup( tnlConfigDescript
if( tnlConfigTagExplicitSolver< ConfigTag, tnlExplicitMersonSolverTag >::enabled )
config.addEntryEnum( "merson" );
}
+ if( tnlConfigTagTimeDiscretisation< ConfigTag, tnlSemiImplicitTimeDiscretisationTag >::enabled )
+ {
+ if( tnlConfigTagSemiImplicitSolver< ConfigTag, tnlSemiImplicitCGSolverTag >::enabled )
+ config.addEntryEnum( "cg" );
+ if( tnlConfigTagSemiImplicitSolver< ConfigTag, tnlSemiImplicitBICGStabSolverTag >::enabled )
+ config.addEntryEnum( "bicgstab" );
+ if( tnlConfigTagSemiImplicitSolver< ConfigTag, tnlSemiImplicitGMRESSolverTag >::enabled )
+ config.addEntryEnum( "gmres" );
+ if( tnlConfigTagSemiImplicitSolver< ConfigTag, tnlSemiImplicitSORSolverTag >::enabled )
+ config.addEntryEnum( "sor" );
+ }
if( tnlConfigTagTimeDiscretisation< ConfigTag, tnlExplicitTimeDiscretisationTag >::enabled )
{
config.addDelimiter( " === Explicit solvers parameters === " );
@@ -115,6 +126,19 @@ bool tnlSolverConfig< ConfigTag, ProblemConfig >::configSetup( tnlConfigDescript
if( tnlConfigTagExplicitSolver< ConfigTag, tnlExplicitMersonSolverTag >::enabled )
tnlMersonSolver< tnlDummyProblem< double, tnlHost, int > >::configSetup( config );
}
+ if( tnlConfigTagTimeDiscretisation< ConfigTag, tnlSemiImplicitTimeDiscretisationTag >::enabled )
+ {
+ config.addDelimiter( " === Semi-implicit solvers parameters === " );
+ typedef tnlCSRMatrix< double, tnlHost, int > MatrixType;
+ if( tnlConfigTagSemiImplicitSolver< ConfigTag, tnlSemiImplicitCGSolverTag >::enabled )
+ tnlCGSolver< MatrixType >::configSetup( config );
+ if( tnlConfigTagSemiImplicitSolver< ConfigTag, tnlSemiImplicitBICGStabSolverTag >::enabled )
+ tnlBICGStabSolver< MatrixType >::configSetup( config );
+ if( tnlConfigTagSemiImplicitSolver< ConfigTag, tnlSemiImplicitGMRESSolverTag >::enabled )
+ tnlGMRESSolver< MatrixType >::configSetup( config );
+ if( tnlConfigTagSemiImplicitSolver< ConfigTag, tnlSemiImplicitSORSolverTag >::enabled )
+ tnlSORSolver< MatrixType >::configSetup( config );
+ }
config.addDelimiter( " === Logs and messages ===" );
config.addEntry< int >( "verbose", "Set the verbose mode. The higher number the more messages are generated.", 1 );
diff --git a/src/implementation/solvers/tnlSolverInitiator_impl.h b/src/implementation/solvers/tnlSolverInitiator_impl.h
index f6bef3fae120f272ac54cd5a267ccb4836e27fa2..f2773db7b90456e155cb346f5b4a54fc85f3a3ac 100644
--- a/src/implementation/solvers/tnlSolverInitiator_impl.h
+++ b/src/implementation/solvers/tnlSolverInitiator_impl.h
@@ -18,6 +18,10 @@
#include <config/tnlParameterContainer.h>
#include <solvers/tnlMeshTypeResolver.h>
#include <solvers/tnlConfigTags.h>
+#include <solvers/linear/stationary/tnlSORSolver.h>
+#include <solvers/linear/krylov/tnlCGSolver.h>
+#include <solvers/linear/krylov/tnlBICGStabSolver.h>
+#include <solvers/linear/krylov/tnlGMRESSolver.h>
#include <core/tnlHost.h>
#include <core/tnlCuda.h>
diff --git a/src/implementation/solvers/tnlSolverStarter_impl.h b/src/implementation/solvers/tnlSolverStarter_impl.h
index 2f78311e96dd4c8690a744709af36249282eda9b..816baec41bd024134c13b0a6115d9d909c5a69f2 100644
--- a/src/implementation/solvers/tnlSolverStarter_impl.h
+++ b/src/implementation/solvers/tnlSolverStarter_impl.h
@@ -147,6 +147,12 @@ class tnlSolverStarterTimeDiscretisationSetter< Problem, tnlSemiImplicitTimeDisc
const tnlParameterContainer& parameters )
{
const tnlString& discreteSolver = parameters. GetParameter< tnlString>( "discrete-solver" );
+ if( discreteSolver == "sor" )
+ return tnlSolverStarterSemiImplicitSolverSetter< Problem, tnlSemiImplicitSORSolverTag, ConfigTag >::run( problem, parameters );
+ if( discreteSolver == "cg" )
+ return tnlSolverStarterSemiImplicitSolverSetter< Problem, tnlSemiImplicitCGSolverTag, ConfigTag >::run( problem, parameters );
+ if( discreteSolver == "bicgstab" )
+ return tnlSolverStarterSemiImplicitSolverSetter< Problem, tnlSemiImplicitBICGStabSolverTag, ConfigTag >::run( problem, parameters );
if( discreteSolver == "gmres" )
return tnlSolverStarterSemiImplicitSolverSetter< Problem, tnlSemiImplicitGMRESSolverTag, ConfigTag >::run( problem, parameters );
return false;
@@ -236,6 +242,57 @@ class tnlSolverStarterSemiImplicitSolverSetter< Problem, SemiImplicitSolver, Con
}
};
+template< typename Problem,
+ typename ConfigTag >
+class tnlSolverStarterSemiImplicitSolverSetter< Problem, tnlSemiImplicitSORSolverTag, ConfigTag, true >
+{
+ public:
+ static bool run( Problem& problem,
+ const tnlParameterContainer& parameters )
+ {
+ typedef typename Problem::MatrixType MatrixType;
+ typedef tnlSORSolver< MatrixType > LinearSystemSolver;
+ typedef tnlSemiImplicitTimeStepper< Problem, LinearSystemSolver > TimeStepper;
+ return tnlSolverStarterSemiImplicitTimeStepperSetter< Problem,
+ TimeStepper,
+ ConfigTag >::run( problem, parameters );
+ }
+};
+
+template< typename Problem,
+ typename ConfigTag >
+class tnlSolverStarterSemiImplicitSolverSetter< Problem, tnlSemiImplicitCGSolverTag, ConfigTag, true >
+{
+ public:
+ static bool run( Problem& problem,
+ const tnlParameterContainer& parameters )
+ {
+ typedef typename Problem::MatrixType MatrixType;
+ typedef tnlCGSolver< MatrixType > LinearSystemSolver;
+ typedef tnlSemiImplicitTimeStepper< Problem, LinearSystemSolver > TimeStepper;
+ return tnlSolverStarterSemiImplicitTimeStepperSetter< Problem,
+ TimeStepper,
+ ConfigTag >::run( problem, parameters );
+ }
+};
+
+template< typename Problem,
+ typename ConfigTag >
+class tnlSolverStarterSemiImplicitSolverSetter< Problem, tnlSemiImplicitBICGStabSolverTag, ConfigTag, true >
+{
+ public:
+ static bool run( Problem& problem,
+ const tnlParameterContainer& parameters )
+ {
+ typedef typename Problem::MatrixType MatrixType;
+ typedef tnlBICGStabSolver< MatrixType > LinearSystemSolver;
+ typedef tnlSemiImplicitTimeStepper< Problem, LinearSystemSolver > TimeStepper;
+ return tnlSolverStarterSemiImplicitTimeStepperSetter< Problem,
+ TimeStepper,
+ ConfigTag >::run( problem, parameters );
+ }
+};
+
template< typename Problem,
typename ConfigTag >
class tnlSolverStarterSemiImplicitSolverSetter< Problem, tnlSemiImplicitGMRESSolverTag, ConfigTag, true >
diff --git a/src/matrices/tnlMatrixSetter.h b/src/matrices/tnlMatrixSetter.h
index 2535389454dc7d23bee3dcfdd1f5b06f04bc958c..cf04a6e2b389d4f186ea626c4e02866d8536a050 100644
--- a/src/matrices/tnlMatrixSetter.h
+++ b/src/matrices/tnlMatrixSetter.h
@@ -126,11 +126,10 @@ class tnlMatrixSetter< tnlGrid< Dimensions, Real, Device, Index >,
{
public:
- template< int EntityDimension >
- void processEntity( const MeshType& mesh,
- TraversalUserData& userData,
- const IndexType index,
- const CoordinatesType& coordinates )
+ void processCell( const MeshType& mesh,
+ TraversalUserData& userData,
+ const IndexType index,
+ const CoordinatesType& coordinates )
{
userData.rowLengths[ index ] =
userData.boundaryConditions.getLinearSystemRowLength( mesh, index, coordinates );
@@ -142,11 +141,10 @@ class tnlMatrixSetter< tnlGrid< Dimensions, Real, Device, Index >,
{
public:
- template< int EntityDimensions >
- void processEntity( const MeshType& mesh,
- TraversalUserData& userData,
- const IndexType index,
- const CoordinatesType& coordinates )
+ void processCell( const MeshType& mesh,
+ TraversalUserData& userData,
+ const IndexType index,
+ const CoordinatesType& coordinates )
{
userData.rowLengths[ index ] =
userData.differentialOperator.getLinearSystemRowLength( mesh, index, coordinates );
diff --git a/src/matrices/tnlMultidiagonalMatrix.h b/src/matrices/tnlMultidiagonalMatrix.h
index 702b4481002262773ade66a5bf32a662dd13ee5e..40f299d031d5eb74a124434adca5f757240a228c 100644
--- a/src/matrices/tnlMultidiagonalMatrix.h
+++ b/src/matrices/tnlMultidiagonalMatrix.h
@@ -63,6 +63,8 @@ class tnlMultidiagonalMatrix : public tnlMatrix< Real, Device, Index >
IndexType getNumberOfNonzeroMatrixElements() const;
+ IndexType getMaxRowlength() const;
+
void reset();
template< typename Real2, typename Device2, typename Index2 >
diff --git a/src/matrices/tnlSparseMatrix.h b/src/matrices/tnlSparseMatrix.h
index 1f5a06b31136391a510894ed58e4ac14e18af7ac..40a8cdfa754bf7cce24fbc1261446a35d8541d10 100644
--- a/src/matrices/tnlSparseMatrix.h
+++ b/src/matrices/tnlSparseMatrix.h
@@ -43,6 +43,8 @@ class tnlSparseMatrix : public tnlMatrix< Real, Device, Index >
IndexType getNumberOfNonzeroMatrixElements() const;
+ IndexType getMaxRowLength() const;
+
#ifdef HAVE_CUDA
__device__ __host__
#endif
@@ -60,8 +62,9 @@ class tnlSparseMatrix : public tnlMatrix< Real, Device, Index >
bool allocateMatrixElements( const IndexType& numberOfMatrixElements );
-
tnlVector< Index, Device, Index > columnIndexes;
+
+ Index maxRowLength;
};
#include <implementation/matrices/tnlSparseMatrix_impl.h>
diff --git a/src/matrices/tnlTridiagonalMatrix.h b/src/matrices/tnlTridiagonalMatrix.h
index b7f77ba2fdd9dfd802ff66ac61c85cb4b7b6c729..3ed1dfd2858e050231ce7f0b3e0517142f0f05a7 100644
--- a/src/matrices/tnlTridiagonalMatrix.h
+++ b/src/matrices/tnlTridiagonalMatrix.h
@@ -60,6 +60,8 @@ class tnlTridiagonalMatrix : public tnlMatrix< Real, Device, Index >
IndexType getNumberOfNonzeroMatrixElements() const;
+ IndexType getMaxRowlength() const;
+
void reset();
template< typename Real2, typename Device2, typename Index2 >
diff --git a/src/operators/diffusion/tnlLinearDiffusion.h b/src/operators/diffusion/tnlLinearDiffusion.h
index 5798a3b6cb8c2a7891ae0606f665e83d262d1728..246943b32de564b4bba3bc48aa51826861395be0 100644
--- a/src/operators/diffusion/tnlLinearDiffusion.h
+++ b/src/operators/diffusion/tnlLinearDiffusion.h
@@ -5,8 +5,8 @@
#include <mesh/tnlGrid.h>
template< typename Mesh,
- typename Real,// = typename Mesh::RealType,
- typename Index >// = typename Mesh::IndexType >
+ typename Real = typename Mesh::RealType,
+ typename Index = typename Mesh::IndexType >
class tnlLinearDiffusion
{
@@ -47,6 +47,21 @@ class tnlLinearDiffusion< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index
const IndexType& index,
const CoordinatesType& coordinates ) const;
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ IndexType* columns,
+ RealType* values,
+ IndexType& rowLength ) const;
+
};
@@ -84,6 +99,20 @@ class tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index
const IndexType& index,
const CoordinatesType& coordinates ) const;
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ IndexType* columns,
+ RealType* values,
+ IndexType& rowLength ) const;
};
@@ -121,6 +150,21 @@ class tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index
const IndexType& index,
const CoordinatesType& coordinates ) const;
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ IndexType* columns,
+ RealType* values,
+ IndexType& rowLength ) const;
+
};
diff --git a/src/operators/tnlDirichletBoundaryConditions.h b/src/operators/tnlDirichletBoundaryConditions.h
index 3cd5e6057f3be08ec54438904f2d16ab74725301..b631df7c3d135ff7a9b354a0c30e6f782bfe58f5 100644
--- a/src/operators/tnlDirichletBoundaryConditions.h
+++ b/src/operators/tnlDirichletBoundaryConditions.h
@@ -55,6 +55,19 @@ class tnlDirichletBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >,
const IndexType& index,
const CoordinatesType& coordinates ) const;
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void updateLinearSystem( const RealType& time,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ DofVectorType& u,
+ DofVectorType& b,
+ IndexType* columns,
+ RealType* values,
+ IndexType& rowLength ) const;
+
protected:
Function function;
@@ -100,6 +113,18 @@ class tnlDirichletBoundaryConditions< tnlGrid< 2, MeshReal, Device, MeshIndex >,
const IndexType& index,
const CoordinatesType& coordinates ) const;
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void updateLinearSystem( const RealType& time,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ DofVectorType& u,
+ DofVectorType& b,
+ IndexType* columns,
+ RealType* values,
+ IndexType& rowLength ) const;
protected:
@@ -147,6 +172,19 @@ class tnlDirichletBoundaryConditions< tnlGrid< 3, MeshReal, Device, MeshIndex >,
const IndexType& index,
const CoordinatesType& coordinates ) const;
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void updateLinearSystem( const RealType& time,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ DofVectorType& u,
+ DofVectorType& b,
+ IndexType* columns,
+ RealType* values,
+ IndexType& rowLength ) const;
+
protected:
Function function;
diff --git a/src/solvers/linear/krylov/tnlBICGStabSolver.h b/src/solvers/linear/krylov/tnlBICGStabSolver.h
index e4cf84e7cb493e5832de0a0d14ed1f63220a5695..f5fc09107deb3e069c7fbdd54763e45dfb9bbe08 100644
--- a/src/solvers/linear/krylov/tnlBICGStabSolver.h
+++ b/src/solvers/linear/krylov/tnlBICGStabSolver.h
@@ -49,6 +49,12 @@ class tnlBICGStabSolver : public tnlObject,
tnlString getType() const;
+ static void configSetup( tnlConfigDescription& config,
+ const tnlString& prefix = "" );
+
+ bool setup( const tnlParameterContainer& parameters,
+ const tnlString& prefix = "" );
+
void setMatrix( const MatrixType& matrix );
void setPreconditioner( const Preconditioner& preconditioner );
diff --git a/src/solvers/linear/krylov/tnlCGSolver.h b/src/solvers/linear/krylov/tnlCGSolver.h
index 6ef497a2cfd5900ccad77362daf4e44bfba9fbed..3f8757d1f28bb4e52d24f89f038f07d1f959dd6f 100644
--- a/src/solvers/linear/krylov/tnlCGSolver.h
+++ b/src/solvers/linear/krylov/tnlCGSolver.h
@@ -47,6 +47,12 @@ class tnlCGSolver : public tnlObject,
tnlString getType() const;
+ static void configSetup( tnlConfigDescription& config,
+ const tnlString& prefix = "" );
+
+ bool setup( const tnlParameterContainer& parameters,
+ const tnlString& prefix = "" );
+
void setMatrix( const MatrixType& matrix );
void setPreconditioner( const Preconditioner& preconditioner );
diff --git a/src/solvers/linear/krylov/tnlGMRESSolver.h b/src/solvers/linear/krylov/tnlGMRESSolver.h
index 4ccee98ee4981307b96e099461611323f3b7c2a8..e9fc5103ece64554aa0e8b60b76f0c58d13bac89 100644
--- a/src/solvers/linear/krylov/tnlGMRESSolver.h
+++ b/src/solvers/linear/krylov/tnlGMRESSolver.h
@@ -47,6 +47,12 @@ class tnlGMRESSolver : public tnlObject,
tnlString getType() const;
+ static void configSetup( tnlConfigDescription& config,
+ const tnlString& prefix = "" );
+
+ bool setup( const tnlParameterContainer& parameters,
+ const tnlString& prefix = "" );
+
void setRestarting( IndexType rest );
void setMatrix( const MatrixType& matrix );
diff --git a/src/solvers/linear/krylov/tnlTFQMRSolver.h b/src/solvers/linear/krylov/tnlTFQMRSolver.h
index e1e519af52b6cd9ef65ebc85f4533e2d50cefce0..0869a990ad6cc4dde2cc8b08116cd69c451c92e5 100644
--- a/src/solvers/linear/krylov/tnlTFQMRSolver.h
+++ b/src/solvers/linear/krylov/tnlTFQMRSolver.h
@@ -49,6 +49,12 @@ class tnlTFQMRSolver : public tnlObject,
tnlString getType() const;
+ static void configSetup( tnlConfigDescription& config,
+ const tnlString& prefix = "" );
+
+ bool setup( const tnlParameterContainer& parameters,
+ const tnlString& prefix = "" );
+
void setMatrix( const MatrixType& matrix );
void setPreconditioner( const Preconditioner& preconditioner );
diff --git a/src/solvers/linear/stationary/tnlSORSolver.h b/src/solvers/linear/stationary/tnlSORSolver.h
index 10aff7b8c794d57e76d4c751f29439a91efd9333..93f853716893b9bb5a5ee5eedd0cd05758e40181 100644
--- a/src/solvers/linear/stationary/tnlSORSolver.h
+++ b/src/solvers/linear/stationary/tnlSORSolver.h
@@ -33,6 +33,7 @@ class tnlSORSolver : public tnlObject,
public tnlIterativeSolver< typename Matrix :: RealType,
typename Matrix :: IndexType >
{
+ public:
typedef typename Matrix :: RealType RealType;
typedef typename Matrix :: IndexType IndexType;
@@ -40,12 +41,16 @@ class tnlSORSolver : public tnlObject,
typedef Matrix MatrixType;
typedef Preconditioner PreconditionerType;
- public:
tnlSORSolver();
tnlString getType() const;
+ static void configSetup( tnlConfigDescription& config,
+ const tnlString& prefix = "" );
+
+ bool setup( const tnlParameterContainer& parameters,
+ const tnlString& prefix = "" );
void setOmega( const RealType& omega );
diff --git a/src/solvers/pde/tnlLinearSystemAssembler.h b/src/solvers/pde/tnlLinearSystemAssembler.h
index 924572017f8984cb2124ede81ebb429af024e38c..def29d602b6ec07425c6105d9e3788667e601f16 100644
--- a/src/solvers/pde/tnlLinearSystemAssembler.h
+++ b/src/solvers/pde/tnlLinearSystemAssembler.h
@@ -27,6 +27,10 @@ template< typename Real,
class tnlLinearSystemAssemblerTraversalUserData
{
public:
+ typedef Matrix MatrixType;
+ typedef typename Matrix::DeviceType DeviceType;
+ typedef tnlVector< typename MatrixType::RealType, DeviceType, typename MatrixType::IndexType > RowValuesType;
+ typedef tnlVector< typename MatrixType::IndexType, DeviceType, typename MatrixType::IndexType > RowColumnsType;
const Real &time;
@@ -42,6 +46,10 @@ class tnlLinearSystemAssemblerTraversalUserData
Matrix &matrix;
+ RowValuesType& values;
+
+ RowColumnsType& columns;
+
tnlLinearSystemAssemblerTraversalUserData( const Real& time,
const Real& tau,
const DifferentialOperator& differentialOperator,
@@ -49,7 +57,9 @@ class tnlLinearSystemAssemblerTraversalUserData
const RightHandSide& rightHandSide,
DofVector& u,
Matrix& matrix,
- DofVector& b )
+ DofVector& b,
+ RowColumnsType& columns,
+ RowValuesType& values )
: time( time ),
tau( tau ),
differentialOperator( differentialOperator ),
@@ -57,7 +67,9 @@ class tnlLinearSystemAssemblerTraversalUserData
rightHandSide( rightHandSide ),
u( u ),
b( b ),
- matrix( matrix )
+ matrix( matrix ),
+ columns( columns ),
+ values( values )
{};
protected:
@@ -109,12 +121,19 @@ class tnlLinearSystemAssembler
TraversalUserData& userData,
const IndexType index )
{
+ typename MatrixType::IndexType rowLength;
userData.boundaryConditions.updateLinearSystem( userData.time,
mesh,
index,
userData.u,
- userData.matrix,
- userData.b );
+ userData.b,
+ userData.columns.getData(),
+ userData.values.getData(),
+ rowLength );
+ userData.matrix.setRowFast( index,
+ userData.columns.getData(),
+ userData.values.getData(),
+ rowLength );
}
};
@@ -131,14 +150,24 @@ class tnlLinearSystemAssembler
TraversalUserData& userData,
const IndexType index )
{
- userData.differentialOperator.updateLinearSystem( mesh,
+ userData.b[ index ] = userData.u[ index ] +
+ userData.tau * userData.rightHandSide.getValue( mesh.getEntityCenter< EntityDimensions >( index ),
+ userData.time );
+ typename MatrixType::IndexType rowLength;
+ userData.differentialOperator.updateLinearSystem( userData.time,
+ userData.tau,
+ mesh,
index,
userData.u,
- userData.matrix,
userData.b,
- userData.time );
- userData.b[ index ] += userData.tau * userData.rightHandSide.getValue( mesh.getEntityCenter< EntityDimensions >( index ),
- userData.time );
+ userData.columns.getData(),
+ userData.values.getData(),
+ rowLength );
+ userData.matrix.setRowFast( index,
+ userData.columns.getData(),
+ userData.values.getData(),
+ rowLength );
+ userData.matrix.addElement( index, index, 1.0, 1.0 );
}
};
@@ -189,22 +218,28 @@ class tnlLinearSystemAssembler< tnlGrid< Dimensions, Real, Device, Index >,
{
public:
- template< int EntityDimension >
#ifdef HAVE_CUDA
__host__ __device__
#endif
- void processEntity( const MeshType& mesh,
- TraversalUserData& userData,
- const IndexType index,
- const CoordinatesType& coordinates )
+ void processCell( const MeshType& mesh,
+ TraversalUserData& userData,
+ const IndexType index,
+ const CoordinatesType& coordinates )
{
+ typename MatrixType::IndexType rowLength;
userData.boundaryConditions.updateLinearSystem( userData.time,
mesh,
index,
coordinates,
userData.u,
- userData.matrix,
- userData.b );
+ userData.b,
+ userData.columns.getData(),
+ userData.values.getData(),
+ rowLength );
+ userData.matrix.setRowFast( index,
+ userData.columns.getData(),
+ userData.values.getData(),
+ rowLength );
}
};
@@ -213,25 +248,33 @@ class tnlLinearSystemAssembler< tnlGrid< Dimensions, Real, Device, Index >,
{
public:
- template< int EntityDimensions >
#ifdef HAVE_CUDA
__host__ __device__
#endif
- void processEntity( const MeshType& mesh,
- TraversalUserData& userData,
- const IndexType index,
- const CoordinatesType& coordinates )
+ void processCell( const MeshType& mesh,
+ TraversalUserData& userData,
+ const IndexType index,
+ const CoordinatesType& coordinates )
{
- userData.differentialOperator.updateLinearSystem( mesh,
+ userData.b[ index ] = userData.u[ index ] +
+ userData.tau * userData.rightHandSide.getValue( mesh.getCellCenter( index ),
+ userData.time );
+ typename MatrixType::IndexType rowLength;
+ userData.differentialOperator.updateLinearSystem( userData.time,
+ userData.tau,
+ mesh,
index,
coordinates,
userData.u,
- userData.matrix,
userData.b,
- userData.time );
- userData.b[ index ] += userData.tau *
- userData.rightHandSide.getValue( mesh.getEntityCenter< EntityDimensions >( index ),
- userData.time );
+ userData.columns.getData(),
+ userData.values.getData(),
+ rowLength );
+ userData.matrix.setRowFast( index,
+ userData.columns.getData(),
+ userData.values.getData(),
+ rowLength );
+ userData.matrix.addElement( index, index, 1.0, 1.0 );
}
};
diff --git a/src/solvers/pde/tnlSemiImplicitTimeStepper.h b/src/solvers/pde/tnlSemiImplicitTimeStepper.h
index be8e25a93dd01cb8658d0c9ff9918bbcb95f54e1..2735144ac2fdd2560f99eba877b0092ab4794edd 100644
--- a/src/solvers/pde/tnlSemiImplicitTimeStepper.h
+++ b/src/solvers/pde/tnlSemiImplicitTimeStepper.h
@@ -71,6 +71,8 @@ class tnlSemiImplicitTimeStepper
LinearSystemSolver* linearSystemSolver;
RealType timeStep;
+
+ bool verbose;
};
#include <implementation/solvers/pde/tnlSemiImplicitTimeStepper_impl.h>
diff --git a/src/solvers/tnlConfigTags.h b/src/solvers/tnlConfigTags.h
index 7fab4b74e014e9653254c2a52c6460493339ea3e..9f1794dcc0652a99c3340a931b1f592ed8687ecb 100644
--- a/src/solvers/tnlConfigTags.h
+++ b/src/solvers/tnlConfigTags.h
@@ -88,6 +88,9 @@ template< typename ConfigTag, typename ExplicitSolver > struct tnlConfigTagExpli
* All semi-implicit solvers are enabled by default
*/
+class tnlSemiImplicitSORSolverTag{};
+class tnlSemiImplicitCGSolverTag{};
+class tnlSemiImplicitBICGStabSolverTag{};
class tnlSemiImplicitGMRESSolverTag{};
template< typename ConfigTag, typename SemiImplicitSolver > struct tnlConfigTagSemiImplicitSolver{ enum { enabled = true }; };