From 73ac70ad6d738476f7739600b673e69770a679d9 Mon Sep 17 00:00:00 2001
From: Tomas Oberhuber <tomas.oberhuber@fjfi.cvut.cz>
Date: Sun, 2 Nov 2014 18:42:50 +0100
Subject: [PATCH] Debuging the heat equation EOC test.
---
TODO | 16 ++++
examples/heat-equation/tnl-heat-equation.cpp | 9 ++-
.../functions/tnlSinWaveFunction_impl.h | 4 +-
.../matrices/tnlMatrixSetter_impl.h | 4 +-
.../mesh/tnlTraversal_Grid1D_impl.h | 19 +++--
.../mesh/tnlTraversal_Grid2D_impl.h | 31 ++++++--
.../mesh/tnlTraversal_Grid3D_impl.h | 37 +++++++++-
.../diffusion/tnlLinearDiffusion_impl.h | 36 +++++----
.../tnlDirichletBoundaryConditions_impl.h | 2 +-
.../pde/tnlSemiImplicitTimeStepper_impl.h | 9 +++
src/matrices/tnlMatrixSetter.h | 12 +--
.../tnlDirichletBoundaryConditions.h | 1 +
.../matrices/tnlSparseMatrixTester.h | 40 +++++-----
.../diffusion/tnlLinearDiffusionTest.cpp | 32 +++++++-
.../operators/tnlPDEOperatorEocTestResult.h | 1 +
.../operators/tnlPDEOperatorEocTestSetter.h | 7 ++
.../operators/tnlPDEOperatorEocTester.h | 14 ++--
tests/unit-tests/tnlApproximationError.h | 60 ++++++++++++++-
tests/unit-tests/tnlApproximationError_impl.h | 74 ++++++++++++++++++-
19 files changed, 335 insertions(+), 73 deletions(-)
diff --git a/TODO b/TODO
index 601def17fe..333e2c0927 100644
--- a/TODO
+++ b/TODO
@@ -1,4 +1,20 @@
+TODO: test aproximace pro implicitni schemata
+TODO: pokusit se napsat test explicitni aproximace podobne jako implicinit bez specializace pro grid
+TODO: doladit semi-implicitni resic
+TODO: opravit argumenty metod v heat equation - nekde se nepredavaji dofy a pomocne dofy
+TODO: doresit okrajove podminky v heat equation - nacitani ze souboru
+ => soucasne Dir. a Neum. okrajove podminky prejmenovat asi tnlAnalytic... a zavest nove tnlDirichlet... a tnlNeumann...
+ ktere budou mit sve hodnoty ulozene ve vektoru
TODO: neumanovy okrajove podminky
+TODO: doladit vse s CUDA
+TODO: doplnit mesh travelsals pro jine mesh entity nez cell
+TODO: implementace maticovych resicu
+ * Gaussova eliminace
+ * SOR metoda
+ * Jacobiho metoda
+ * TFQMR metoda
+ * IDR metody
+
TODO: implementovat tridu tnlFileName pro generovani jmen souboru
diff --git a/examples/heat-equation/tnl-heat-equation.cpp b/examples/heat-equation/tnl-heat-equation.cpp
index 449eb20a30..d74f517f68 100644
--- a/examples/heat-equation/tnl-heat-equation.cpp
+++ b/examples/heat-equation/tnl-heat-equation.cpp
@@ -32,8 +32,13 @@ class heatEquationConfig
public:
static void configSetup( tnlConfigDescription& config )
{
- //config.addDelimiter( "Heat equation settings:" );
- }
+ config.addDelimiter( "Heat equation settings:" );
+ config.addEntry< tnlString >( "boundary-conditions-type", "Choose the boundary conditions type.", "dirichlet");
+ config.addEntryEnum< tnlString >( "dirichlet" );
+ config.addEntryEnum< tnlString >( "neumann" );
+ config.addEntry< tnlString >( "boundary-conditions", "File with the values of the boundary conditions.", "boundary.tnl" );
+ config.addEntry< tnlString >( "initial-condition", "File with the initial condition.", "initial.tnl");
+ };
};
template< typename Real,
diff --git a/src/implementation/functions/tnlSinWaveFunction_impl.h b/src/implementation/functions/tnlSinWaveFunction_impl.h
index 10f0afea48..90f6f003cd 100644
--- a/src/implementation/functions/tnlSinWaveFunction_impl.h
+++ b/src/implementation/functions/tnlSinWaveFunction_impl.h
@@ -129,11 +129,11 @@ getValue( const Vertex& v,
if( XDiffOrder == 1 && YDiffOrder == 0 )
return 2.0 * M_PI * x / ( this->waveLength * sqrt( x * x + y * y ) ) * this->amplitude * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y ) / this->waveLength );
if( XDiffOrder == 2 && YDiffOrder == 0 )
- return 2.0 * M_PI * y * y / ( this->waveLength * sqrt( x * x + y * y ) * sqrt( x * x + y * y ) * sqrt( x * x + y * y ) ) * this->amplitude * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y ) / this->waveLength ) - 4.0 * M_PI * M_PI * x * x / ( this->waveLength * this->waveLength * ( x * x + y * y ) ) * this->amplitude * sin( this->phase + 2.0 * M_PI * sqrt( x * x + y * y ) / this->waveLength );
+ return 2.0 * M_PI * x * x / ( this->waveLength * sqrt( x * x + y * y ) * sqrt( x * x + y * y ) * sqrt( x * x + y * y ) ) * this->amplitude * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y ) / this->waveLength ) - 4.0 * M_PI * M_PI * x * x / ( this->waveLength * this->waveLength * ( x * x + y * y ) ) * this->amplitude * sin( this->phase + 2.0 * M_PI * sqrt( x * x + y * y ) / this->waveLength );
if( XDiffOrder == 0 && YDiffOrder == 1 )
return 2.0 * M_PI * y / ( this->waveLength * sqrt( x * x + y * y ) ) * this->amplitude * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y ) / this->waveLength );
if( XDiffOrder == 0 && YDiffOrder == 2 )
- return 2.0 * M_PI * x * x / ( this->waveLength * sqrt( x * x + y * y ) * sqrt( x * x + y * y ) * sqrt( x * x + y * y ) ) * this->amplitude * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y ) / this->waveLength ) - 4.0 * M_PI * M_PI * y * y / ( this->waveLength * this->waveLength * ( x * x + y * y ) ) * this->amplitude * sin( this->phase + 2.0 * M_PI * sqrt( x * x + y * y ) / this->waveLength );
+ return 2.0 * M_PI * y * y / ( this->waveLength * sqrt( x * x + y * y ) * sqrt( x * x + y * y ) * sqrt( x * x + y * y ) ) * this->amplitude * cos( this->phase + 2.0 * M_PI * sqrt( x * x + y * y ) / this->waveLength ) - 4.0 * M_PI * M_PI * y * y / ( this->waveLength * this->waveLength * ( x * x + y * y ) ) * this->amplitude * sin( this->phase + 2.0 * M_PI * sqrt( x * x + y * y ) / this->waveLength );
return 0.0;
}
diff --git a/src/implementation/matrices/tnlMatrixSetter_impl.h b/src/implementation/matrices/tnlMatrixSetter_impl.h
index b8a808013f..d1e310eea6 100644
--- a/src/implementation/matrices/tnlMatrixSetter_impl.h
+++ b/src/implementation/matrices/tnlMatrixSetter_impl.h
@@ -54,8 +54,8 @@ template< int Dimensions,
void
tnlMatrixSetter< tnlGrid< Dimensions, Real, Device, Index >, DifferentialOperator, BoundaryConditions, RowLengthsVector >::
getRowLengths( const MeshType& mesh,
- DifferentialOperator& differentialOperator,
- BoundaryConditions& boundaryConditions,
+ const DifferentialOperator& differentialOperator,
+ const BoundaryConditions& boundaryConditions,
RowLengthsVector& rowLengths ) const
{
TraversalUserData userData( differentialOperator, boundaryConditions, rowLengths );
diff --git a/src/implementation/mesh/tnlTraversal_Grid1D_impl.h b/src/implementation/mesh/tnlTraversal_Grid1D_impl.h
index 99259fa78e..b17f5133d4 100644
--- a/src/implementation/mesh/tnlTraversal_Grid1D_impl.h
+++ b/src/implementation/mesh/tnlTraversal_Grid1D_impl.h
@@ -35,13 +35,22 @@ processEntities( const GridType& grid,
*/
CoordinatesType coordinates;
const IndexType& xSize = grid.getDimensions().x();
- for( coordinates.x() = 0; coordinates.x() < xSize; coordinates.x() ++ )
+
+ /****
+ * Boundary conditions
+ */
+ coordinates.x() = 0;
+ boundaryEntitiesProcessor.template processCell( grid, userData, 0, coordinates );
+ coordinates.x() = xSize - 1;
+ boundaryEntitiesProcessor.template processCell( grid, userData, xSize - 1, coordinates );
+
+ /****
+ * Interior cells
+ */
+ for( coordinates.x() = 1; coordinates.x() < xSize-1; coordinates.x() ++ )
{
const IndexType index = grid.getCellIndex( coordinates );
- if( grid.isBoundaryCell( coordinates ) )
- boundaryEntitiesProcessor.template processCell( grid, userData, index, coordinates );
- else
- interiorEntitiesProcessor.template processCell( grid, userData, index, coordinates );
+ interiorEntitiesProcessor.template processCell( grid, userData, index, coordinates );
}
}
diff --git a/src/implementation/mesh/tnlTraversal_Grid2D_impl.h b/src/implementation/mesh/tnlTraversal_Grid2D_impl.h
index 2f1bb1ee02..a5e6f34609 100644
--- a/src/implementation/mesh/tnlTraversal_Grid2D_impl.h
+++ b/src/implementation/mesh/tnlTraversal_Grid2D_impl.h
@@ -37,17 +37,36 @@ processEntities( const GridType& grid,
CoordinatesType coordinates;
const IndexType& xSize = grid.getDimensions().x();
const IndexType& ySize = grid.getDimensions().y();
+
+ /****
+ * Boundary conditions
+ */
+ for( coordinates.x() = 0; coordinates.x() < xSize; coordinates.x() ++ )
+ {
+ coordinates.y() = 0;
+ boundaryEntitiesProcessor.template processCell( grid, userData, grid.getCellIndex( coordinates ), coordinates );
+ coordinates.y() = ySize - 1;
+ boundaryEntitiesProcessor.template processCell( grid, userData, grid.getCellIndex( coordinates ), coordinates );
+ }
+ for( coordinates.y() = 1; coordinates.y() < ySize - 1; coordinates.y() ++ )
+ {
+ coordinates.x() = 0;
+ boundaryEntitiesProcessor.template processCell( grid, userData, grid.getCellIndex( coordinates ), coordinates );
+ coordinates.x() = xSize - 1;
+ boundaryEntitiesProcessor.template processCell( grid, userData, grid.getCellIndex( coordinates ), coordinates );
+ }
+
+ /****
+ * Interior cells
+ */
#ifdef HAVE_OPENMP
//#pragma omp parallel for
#endif
- for( coordinates.y() = 0; coordinates.y() < ySize; coordinates.y() ++ )
- for( coordinates.x() = 0; coordinates.x() < xSize; coordinates.x() ++ )
+ for( coordinates.y() = 1; coordinates.y() < ySize - 1; coordinates.y() ++ )
+ for( coordinates.x() = 1; coordinates.x() < xSize - 1; coordinates.x() ++ )
{
const IndexType index = grid.getCellIndex( coordinates );
- if( grid.isBoundaryCell( coordinates ) )
- boundaryEntitiesProcessor.template processCell( grid, userData, index, coordinates );
- else
- interiorEntitiesProcessor.template processCell( grid, userData, index, coordinates );
+ interiorEntitiesProcessor.template processCell( grid, userData, index, coordinates );
}
}
diff --git a/src/implementation/mesh/tnlTraversal_Grid3D_impl.h b/src/implementation/mesh/tnlTraversal_Grid3D_impl.h
index d0f21814f6..e64353b0f0 100644
--- a/src/implementation/mesh/tnlTraversal_Grid3D_impl.h
+++ b/src/implementation/mesh/tnlTraversal_Grid3D_impl.h
@@ -37,6 +37,38 @@ processEntities( const GridType& grid,
const IndexType& xSize = grid.getDimensions().x();
const IndexType& ySize = grid.getDimensions().y();
const IndexType& zSize = grid.getDimensions().z();
+
+ /****
+ * Boundary conditions
+ */
+ for( coordinates.y() = 0; coordinates.y() < ySize; coordinates.y() ++ )
+ for( coordinates.x() = 0; coordinates.x() < xSize; coordinates.x() ++ )
+ {
+ coordinates.z() = 0;
+ boundaryEntitiesProcessor.template processCell( grid, userData, grid.getCellIndex( coordinates ), coordinates );
+ coordinates.z() = zSize - 1;
+ boundaryEntitiesProcessor.template processCell( grid, userData, grid.getCellIndex( coordinates ), coordinates );
+ }
+
+ for( coordinates.z() = 0; coordinates.z() < zSize; coordinates.z() ++ )
+ for( coordinates.x() = 0; coordinates.x() < xSize; coordinates.x() ++ )
+ {
+ coordinates.y() = 0;
+ boundaryEntitiesProcessor.template processCell( grid, userData, grid.getCellIndex( coordinates ), coordinates );
+ coordinates.y() = ySize - 1;
+ boundaryEntitiesProcessor.template processCell( grid, userData, grid.getCellIndex( coordinates ), coordinates );
+ }
+
+ for( coordinates.z() = 0; coordinates.z() < zSize; coordinates.z() ++ )
+ for( coordinates.y() = 0; coordinates.y() < ySize; coordinates.y() ++ )
+ {
+ coordinates.x() = 0;
+ boundaryEntitiesProcessor.template processCell( grid, userData, grid.getCellIndex( coordinates ), coordinates );
+ coordinates.x() = xSize - 1;
+ boundaryEntitiesProcessor.template processCell( grid, userData, grid.getCellIndex( coordinates ), coordinates );
+ }
+
+
#ifdef HAVE_OPENMP
//#pragma omp parallel for
#endif
@@ -45,10 +77,7 @@ processEntities( const GridType& grid,
for( coordinates.x() = 0; coordinates.x() < xSize; coordinates.x() ++ )
{
const IndexType index = grid.getCellIndex( coordinates );
- if( grid.isBoundaryCell( coordinates ) )
- boundaryEntitiesProcessor.template processCell( grid, userData, index, coordinates );
- else
- interiorEntitiesProcessor.template processCell( grid, userData, index, coordinates );
+ interiorEntitiesProcessor.template processCell( grid, userData, index, coordinates );
}
}
diff --git a/src/implementation/operators/diffusion/tnlLinearDiffusion_impl.h b/src/implementation/operators/diffusion/tnlLinearDiffusion_impl.h
index 2acc343c26..05ade1d5db 100644
--- a/src/implementation/operators/diffusion/tnlLinearDiffusion_impl.h
+++ b/src/implementation/operators/diffusion/tnlLinearDiffusion_impl.h
@@ -81,12 +81,13 @@ updateLinearSystem( const RealType& time,
RealType* values,
IndexType& rowLength ) const
{
+ const RealType lambdaX = tau * mesh.getHxSquareInverse();
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;
+ values[ 0 ] = -lambdaX;
+ values[ 1 ] = 1.0 + 2.0 * lambdaX;
+ values[ 2 ] = -lambdaX;
rowLength = 3;
}
@@ -170,16 +171,18 @@ updateLinearSystem( const RealType& time,
RealType* values,
IndexType& rowLength ) const
{
+ const RealType lambdaX = tau * mesh.getHxSquareInverse();
+ const RealType lambdaY = tau * mesh.getHySquareInverse();
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;
+ values[ 0 ] = -lambdaY;
+ values[ 1 ] = -lambdaX;
+ values[ 2 ] = 1.0 + 2.0 * ( lambdaX + lambdaY );
+ values[ 3 ] = -lambdaX;
+ values[ 4 ] = -lambdaY;
rowLength = 5;
}
@@ -266,6 +269,9 @@ updateLinearSystem( const RealType& time,
RealType* values,
IndexType& rowLength ) const
{
+ const RealType lambdaX = tau * mesh.getHxSquareInverse();
+ const RealType lambdaY = tau * mesh.getHySquareInverse();
+ const RealType lambdaZ = tau * mesh.getHzSquareInverse();
columns[ 0 ] = mesh.getCellZPredecessor( index );
columns[ 1 ] = mesh.getCellYPredecessor( index );
columns[ 2 ] = mesh.getCellXPredecessor( index );
@@ -273,13 +279,13 @@ updateLinearSystem( const RealType& time,
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;
+ values[ 0 ] = -lambdaZ;
+ values[ 1 ] = -lambdaY;
+ values[ 2 ] = -lambdaX;
+ values[ 3 ] = 1.0 + 2.0 * ( lambdaX + lambdaY + lambdaZ );
+ values[ 4 ] = -lambdaX;
+ values[ 5 ] = -lambdaY;
+ values[ 6 ] = -lambdaZ;
rowLength = 7;
}
diff --git a/src/implementation/operators/tnlDirichletBoundaryConditions_impl.h b/src/implementation/operators/tnlDirichletBoundaryConditions_impl.h
index 6c781e00db..a00a250c79 100644
--- a/src/implementation/operators/tnlDirichletBoundaryConditions_impl.h
+++ b/src/implementation/operators/tnlDirichletBoundaryConditions_impl.h
@@ -11,7 +11,7 @@ template< typename MeshReal,
bool
tnlDirichletBoundaryConditions< tnlGrid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >::
setup( const tnlParameterContainer& parameters,
- const tnlString& prefix )
+ const tnlString& prefix )
{
return function.setup( parameters, prefix );
}
diff --git a/src/implementation/solvers/pde/tnlSemiImplicitTimeStepper_impl.h b/src/implementation/solvers/pde/tnlSemiImplicitTimeStepper_impl.h
index 91f775f6e7..8f44de2620 100644
--- a/src/implementation/solvers/pde/tnlSemiImplicitTimeStepper_impl.h
+++ b/src/implementation/solvers/pde/tnlSemiImplicitTimeStepper_impl.h
@@ -147,6 +147,8 @@ solve( const RealType& time,
dofVector,
this->matrix,
this->rightHandSide );
+ cout << "u = " << dofVector << endl;
+ cout << "b = " << this->rightHandSide << endl;
if( verbose )
cout << " Solving the linear system for time " << t << " \r" << flush;
if( ! this->linearSystemSolver->solve( this->rightHandSide, dofVector ) )
@@ -154,6 +156,13 @@ solve( const RealType& time,
cerr << "The linear system solver did not converge." << endl;
return false;
}
+ DofVectorType aux;
+ aux.setLike( dofVector );
+ this->matrix.vectorProduct( dofVector, aux );
+ cout << "tau = " << currentTau << endl;
+ cout << "u = " << dofVector << endl;
+ cout << "Ax = " << aux << endl;
+ cout << "b = " << this->rightHandSide << endl;
t += currentTau;
}
return true;
diff --git a/src/matrices/tnlMatrixSetter.h b/src/matrices/tnlMatrixSetter.h
index cf04a6e2b3..01448ea314 100644
--- a/src/matrices/tnlMatrixSetter.h
+++ b/src/matrices/tnlMatrixSetter.h
@@ -25,14 +25,14 @@ class tnlMatrixSetterTraversalUserData
{
public:
- DifferentialOperator& differentialOperator;
+ const DifferentialOperator& differentialOperator;
- BoundaryConditions& boundaryConditions;
+ const BoundaryConditions& boundaryConditions;
RowLengthsVector& rowLengths;
- tnlMatrixSetterTraversalUserData( DifferentialOperator& differentialOperator,
- BoundaryConditions& boundaryConditions,
+ tnlMatrixSetterTraversalUserData( const DifferentialOperator& differentialOperator,
+ const BoundaryConditions& boundaryConditions,
RowLengthsVector& rowLengths )
: differentialOperator( differentialOperator ),
boundaryConditions( boundaryConditions ),
@@ -117,8 +117,8 @@ class tnlMatrixSetter< tnlGrid< Dimensions, Real, Device, Index >,
template< int EntityDimensions >
void getRowLengths( const MeshType& mesh,
- DifferentialOperator& differentialOperator,
- BoundaryConditions& boundaryConditions,
+ const DifferentialOperator& differentialOperator,
+ const BoundaryConditions& boundaryConditions,
RowLengthsVector& rowLengths ) const;
diff --git a/src/operators/tnlDirichletBoundaryConditions.h b/src/operators/tnlDirichletBoundaryConditions.h
index b631df7c3d..24c9a7379b 100644
--- a/src/operators/tnlDirichletBoundaryConditions.h
+++ b/src/operators/tnlDirichletBoundaryConditions.h
@@ -4,6 +4,7 @@
#include <core/vectors/tnlStaticVector.h>
#include <config/tnlParameterContainer.h>
#include <functions/tnlConstantFunction.h>
+#include <core/vectors/tnlSharedVector.h>
template< typename Mesh,
typename Function = tnlConstantFunction< Mesh::Dimensions,
diff --git a/tests/unit-tests/matrices/tnlSparseMatrixTester.h b/tests/unit-tests/matrices/tnlSparseMatrixTester.h
index 2b25a67d4d..5d12983eb8 100644
--- a/tests/unit-tests/matrices/tnlSparseMatrixTester.h
+++ b/tests/unit-tests/matrices/tnlSparseMatrixTester.h
@@ -259,12 +259,12 @@ class tnlSparseMatrixTester : public CppUnit :: TestCase
rowLengths.setValue( 7 );
m.setRowLengths( rowLengths );
- if( DeviceType::DeviceType == tnlHostDevice )
+ if( DeviceType::DeviceType == ( int ) tnlHostDevice )
{
for( int i = 0; i < 10; i++ )
m.setElementFast( i, i, i );
}
- if( DeviceType::DeviceType == tnlCudaDevice )
+ if( DeviceType::DeviceType == ( int ) tnlCudaDevice )
{
#ifdef HAVE_CUDA
MatrixType* kernel_matrix = tnlCuda::passToDevice( m );
@@ -340,7 +340,7 @@ class tnlSparseMatrixTester : public CppUnit :: TestCase
rowLengths.setValue( 10 );
m.setRowLengths( rowLengths );
- if( DeviceType::DeviceType == tnlHostDevice )
+ if( DeviceType::DeviceType == ( int ) tnlHostDevice )
{
for( int i = 0; i < 10; i++ )
m.setElementFast( i, i, i );
@@ -348,7 +348,7 @@ class tnlSparseMatrixTester : public CppUnit :: TestCase
for( int j = 0; j < 10; j++ )
m.addElementFast( i, j, 1, 0.5 );
}
- if( DeviceType::DeviceType == tnlCudaDevice )
+ if( DeviceType::DeviceType == ( int ) tnlCudaDevice )
{
#ifdef HAVE_CUDA
MatrixType* kernel_matrix = tnlCuda::passToDevice( m );
@@ -377,13 +377,13 @@ class tnlSparseMatrixTester : public CppUnit :: TestCase
m.reset();
m.setDimensions( 10, 10 );
m.setRowLengths( rowLengths );
- if( DeviceType::DeviceType == tnlHostDevice )
+ if( DeviceType::DeviceType == ( int ) tnlHostDevice )
{
for( int i = 9; i >= 0; i-- )
for( int j = 9; j >= 0; j-- )
m.setElementFast( i, j, i+j );
}
- if( DeviceType::DeviceType == tnlCudaDevice )
+ if( DeviceType::DeviceType == ( int ) tnlCudaDevice )
{
#ifdef HAVE_CUDA
MatrixType* kernel_matrix = tnlCuda::passToDevice( m );
@@ -456,13 +456,13 @@ class tnlSparseMatrixTester : public CppUnit :: TestCase
rowLengths.setElement( i, i+1 );
m.setRowLengths( rowLengths );
- if( DeviceType::DeviceType == tnlHostDevice )
+ if( DeviceType::DeviceType == ( int ) tnlHostDevice )
{
for( int i = 0; i < 10; i++ )
for( int j = 0; j <= i; j++ )
m.setElementFast( i, j, i + j );
}
- if( DeviceType::DeviceType == tnlCudaDevice )
+ if( DeviceType::DeviceType == ( int ) tnlCudaDevice )
{
#ifdef HAVE_CUDA
MatrixType* kernel_matrix = tnlCuda::passToDevice( m );
@@ -491,13 +491,13 @@ class tnlSparseMatrixTester : public CppUnit :: TestCase
m.reset();
m.setDimensions( 10, 10 );
m.setRowLengths( rowLengths );
- if( DeviceType::DeviceType == tnlHostDevice )
+ if( DeviceType::DeviceType == ( int ) tnlHostDevice )
{
for( int i = 9; i >= 0; i-- )
for( int j = i; j >= 0; j-- )
m.setElementFast( i, j, i + j );
}
- if( DeviceType::DeviceType == tnlCudaDevice )
+ if( DeviceType::DeviceType == ( int ) tnlCudaDevice )
{
#ifdef HAVE_CUDA
MatrixType* kernel_matrix = tnlCuda::passToDevice( m );
@@ -597,7 +597,7 @@ class tnlSparseMatrixTester : public CppUnit :: TestCase
m.setRowLengths( rowLengths );
- if( DeviceType::DeviceType == tnlHostDevice )
+ if( DeviceType::DeviceType == ( int ) tnlHostDevice )
{
RealType values[ 1 ];
IndexType columnIndexes[ 1 ];
@@ -608,7 +608,7 @@ class tnlSparseMatrixTester : public CppUnit :: TestCase
m.setRowFast( i, columnIndexes, values, 1 );
}
}
- if( DeviceType::DeviceType == tnlCudaDevice )
+ if( DeviceType::DeviceType == ( int ) tnlCudaDevice )
{
#ifdef HAVE_CUDA
MatrixType* kernel_matrix = tnlCuda::passToDevice( m );
@@ -702,7 +702,7 @@ class tnlSparseMatrixTester : public CppUnit :: TestCase
for( int i = 0; i < 10; i++ )
columnIndexes[ i ] = i;
- if( DeviceType::DeviceType == tnlHostDevice )
+ if( DeviceType::DeviceType == ( int ) tnlHostDevice )
{
for( int i = 0; i < 10; i++ )
{
@@ -715,7 +715,7 @@ class tnlSparseMatrixTester : public CppUnit :: TestCase
m.setRowFast( i, columnIndexes, values, 10 );
}
}
- if( DeviceType::DeviceType == tnlCudaDevice )
+ if( DeviceType::DeviceType == ( int ) tnlCudaDevice )
{
#ifdef HAVE_CUDA
MatrixType* kernel_matrix = tnlCuda::passToDevice( m );
@@ -746,7 +746,7 @@ class tnlSparseMatrixTester : public CppUnit :: TestCase
m.setDimensions( 10, 10 );
m.setRowLengths( rowLengths );
- if( DeviceType::DeviceType == tnlHostDevice )
+ if( DeviceType::DeviceType == ( int ) tnlHostDevice )
{
for( int i = 9; i >= 0; i-- )
{
@@ -755,7 +755,7 @@ class tnlSparseMatrixTester : public CppUnit :: TestCase
m.setRowFast( i, columnIndexes, values, 10 );
}
}
- if( DeviceType::DeviceType == tnlCudaDevice )
+ if( DeviceType::DeviceType == ( int ) tnlCudaDevice )
{
#ifdef HAVE_CUDA
MatrixType* kernel_matrix = tnlCuda::passToDevice( m );
@@ -847,7 +847,7 @@ class tnlSparseMatrixTester : public CppUnit :: TestCase
for( int i = 0; i < 10; i++ )
columnIndexes[ i ] = i;
- if( DeviceType::DeviceType == tnlHostDevice )
+ if( DeviceType::DeviceType == ( int ) tnlHostDevice )
{
for( int i = 0; i < 10; i++ )
{
@@ -856,7 +856,7 @@ class tnlSparseMatrixTester : public CppUnit :: TestCase
m.setRowFast( i, columnIndexes, values, i + 1 );
}
}
- if( DeviceType::DeviceType == tnlCudaDevice )
+ if( DeviceType::DeviceType == ( int ) tnlCudaDevice )
{
#ifdef HAVE_CUDA
MatrixType* kernel_matrix = tnlCuda::passToDevice( m );
@@ -887,7 +887,7 @@ class tnlSparseMatrixTester : public CppUnit :: TestCase
m.setDimensions( 10, 10 );
m.setRowLengths( rowLengths );
- if( DeviceType::DeviceType == tnlHostDevice )
+ if( DeviceType::DeviceType == ( int ) tnlHostDevice )
{
for( int i = 9; i >= 0; i-- )
{
@@ -896,7 +896,7 @@ class tnlSparseMatrixTester : public CppUnit :: TestCase
m.setRowFast( i, columnIndexes, values, i + 1 );
}
}
- if( DeviceType::DeviceType == tnlCudaDevice )
+ if( DeviceType::DeviceType == ( int ) tnlCudaDevice )
{
#ifdef HAVE_CUDA
MatrixType* kernel_matrix = tnlCuda::passToDevice( m );
diff --git a/tests/unit-tests/operators/diffusion/tnlLinearDiffusionTest.cpp b/tests/unit-tests/operators/diffusion/tnlLinearDiffusionTest.cpp
index 525725954c..ddffcf8229 100644
--- a/tests/unit-tests/operators/diffusion/tnlLinearDiffusionTest.cpp
+++ b/tests/unit-tests/operators/diffusion/tnlLinearDiffusionTest.cpp
@@ -31,8 +31,10 @@ template< int Dimensions,
typename Real,
typename Device,
typename Index,
- typename TestFunction >
+ typename TestFunction,
+ typename ApproximationMethod >
class tnlPDEOperatorEocTestResult< tnlLinearDiffusion< tnlGrid< Dimensions, Real, Device, Index >, Real, Index >,
+ ApproximationMethod,
TestFunction >
{
public:
@@ -52,19 +54,47 @@ int main( int argc, char* argv[] )
const bool verbose( true );
const int MeshSize( 32 );
#ifdef HAVE_CPPUNIT
+ /****
+ * Explicit approximation
+ */
if( ! tnlUnitTestStarter :: run< tnlPDEOperatorEocTester< tnlLinearDiffusion< tnlGrid< 1, double, tnlHost, int >, double, int >,
tnlExactLinearDiffusion< 1 >,
tnlExpBumpFunction< 1, double >,
+ tnlExplicitApproximation,
MeshSize,
verbose > >() ||
! tnlUnitTestStarter :: run< tnlPDEOperatorEocTester< tnlLinearDiffusion< tnlGrid< 2, double, tnlHost, int >, double, int >,
tnlExactLinearDiffusion< 2 >,
tnlExpBumpFunction< 2, double >,
+ tnlExplicitApproximation,
MeshSize,
verbose > >() ||
! tnlUnitTestStarter :: run< tnlPDEOperatorEocTester< tnlLinearDiffusion< tnlGrid< 3, double, tnlHost, int >, double, int >,
tnlExactLinearDiffusion< 3 >,
tnlExpBumpFunction< 3, double >,
+ tnlExplicitApproximation,
+ MeshSize,
+ verbose > >() )
+ return EXIT_FAILURE;
+ /****
+ * Implicit (matrix) approximation
+ */
+ if( ! tnlUnitTestStarter :: run< tnlPDEOperatorEocTester< tnlLinearDiffusion< tnlGrid< 1, double, tnlHost, int >, double, int >,
+ tnlExactLinearDiffusion< 1 >,
+ tnlExpBumpFunction< 1, double >,
+ tnlImplicitApproximation,
+ MeshSize,
+ verbose > >() ||
+ ! tnlUnitTestStarter :: run< tnlPDEOperatorEocTester< tnlLinearDiffusion< tnlGrid< 2, double, tnlHost, int >, double, int >,
+ tnlExactLinearDiffusion< 2 >,
+ tnlExpBumpFunction< 2, double >,
+ tnlImplicitApproximation,
+ MeshSize,
+ verbose > >() ||
+ ! tnlUnitTestStarter :: run< tnlPDEOperatorEocTester< tnlLinearDiffusion< tnlGrid< 3, double, tnlHost, int >, double, int >,
+ tnlExactLinearDiffusion< 3 >,
+ tnlExpBumpFunction< 3, double >,
+ tnlImplicitApproximation,
MeshSize,
verbose > >()
)
diff --git a/tests/unit-tests/operators/tnlPDEOperatorEocTestResult.h b/tests/unit-tests/operators/tnlPDEOperatorEocTestResult.h
index 0214a086e0..7e1954dcf2 100644
--- a/tests/unit-tests/operators/tnlPDEOperatorEocTestResult.h
+++ b/tests/unit-tests/operators/tnlPDEOperatorEocTestResult.h
@@ -19,6 +19,7 @@
#define TNLPDEOPERATOREOCTESTRESULT_H_
template< typename ApproximateOperator,
+ typename ApproximationMethod,
typename TestFunction >
class tnlPDEOperatorEocTestResult
{
diff --git a/tests/unit-tests/operators/tnlPDEOperatorEocTestSetter.h b/tests/unit-tests/operators/tnlPDEOperatorEocTestSetter.h
index 2ed0ce0a0e..16542c201d 100644
--- a/tests/unit-tests/operators/tnlPDEOperatorEocTestSetter.h
+++ b/tests/unit-tests/operators/tnlPDEOperatorEocTestSetter.h
@@ -23,6 +23,7 @@
template< typename ApproximateOperator,
typename ExactOperator,
+ typename ApproximationMethod,
typename Mesh,
typename TestFunction >
class tnlPDEOperatorEocTestSetter
@@ -31,11 +32,13 @@ class tnlPDEOperatorEocTestSetter
template< typename ApproximateOperator,
typename ExactOperator,
+ typename ApproximationMethod,
typename Real,
typename Device,
typename Index >
class tnlPDEOperatorEocTestSetter< ApproximateOperator,
ExactOperator,
+ ApproximationMethod,
tnlGrid< 1, Real, Device, Index >,
tnlExpBumpFunction< 1, Real > >
{
@@ -72,11 +75,13 @@ class tnlPDEOperatorEocTestSetter< ApproximateOperator,
template< typename ApproximateOperator,
typename ExactOperator,
+ typename ApproximationMethod,
typename Real,
typename Device,
typename Index >
class tnlPDEOperatorEocTestSetter< ApproximateOperator,
ExactOperator,
+ ApproximationMethod,
tnlGrid< 2, Real, Device, Index >,
tnlExpBumpFunction< 2, Real > >
{
@@ -116,11 +121,13 @@ class tnlPDEOperatorEocTestSetter< ApproximateOperator,
template< typename ApproximateOperator,
typename ExactOperator,
+ typename ApproximationMethod,
typename Real,
typename Device,
typename Index >
class tnlPDEOperatorEocTestSetter< ApproximateOperator,
ExactOperator,
+ ApproximationMethod,
tnlGrid< 3, Real, Device, Index >,
tnlExpBumpFunction< 3, Real > >
{
diff --git a/tests/unit-tests/operators/tnlPDEOperatorEocTester.h b/tests/unit-tests/operators/tnlPDEOperatorEocTester.h
index 9e1b50ff54..97381cf3ec 100644
--- a/tests/unit-tests/operators/tnlPDEOperatorEocTester.h
+++ b/tests/unit-tests/operators/tnlPDEOperatorEocTester.h
@@ -31,18 +31,19 @@
template< typename ApproximateOperator,
typename ExactOperator,
typename TestFunction,
+ typename ApproximationMethod,
int MeshSize,
bool verbose = false >
class tnlPDEOperatorEocTester : public CppUnit :: TestCase
{
public:
- typedef tnlPDEOperatorEocTester< ApproximateOperator, ExactOperator, TestFunction, MeshSize, verbose > TesterType;
+ typedef tnlPDEOperatorEocTester< ApproximateOperator, ExactOperator, TestFunction, ApproximationMethod, MeshSize, verbose > TesterType;
typedef typename CppUnit::TestCaller< TesterType > TestCallerType;
typedef typename ApproximateOperator::MeshType MeshType;
typedef typename ApproximateOperator::RealType RealType;
typedef typename ApproximateOperator::IndexType IndexType;
- typedef tnlPDEOperatorEocTestSetter< ApproximateOperator, ExactOperator, MeshType, TestFunction > TestSetter;
- typedef tnlPDEOperatorEocTestResult< ApproximateOperator, TestFunction > TestResult;
+ typedef tnlPDEOperatorEocTestSetter< ApproximateOperator, ExactOperator, ApproximationMethod, MeshType, TestFunction > TestSetter;
+ typedef tnlPDEOperatorEocTestResult< ApproximateOperator, ApproximationMethod, TestFunction > TestResult;
tnlPDEOperatorEocTester(){};
@@ -55,6 +56,7 @@ class tnlPDEOperatorEocTester : public CppUnit :: TestCase
ApproximateOperator::getType() + ", " +
ExactOperator::getType() + ", " +
TestFunction::getType() + ", " +
+ ApproximationMethod::getType() + ", " +
tnlString( MeshSize ) + ", " +
tnlString( verbose ) + " >";
CppUnit :: TestSuite* suiteOfTests = new CppUnit :: TestSuite( testName.getString() );
@@ -80,7 +82,8 @@ class tnlPDEOperatorEocTester : public CppUnit :: TestCase
tnlApproximationError< MeshType,
ExactOperator,
ApproximateOperator,
- TestFunction >
+ TestFunction,
+ ApproximationMethod >
::getError( mesh,
exactOperator,
approximateOperator,
@@ -94,7 +97,8 @@ class tnlPDEOperatorEocTester : public CppUnit :: TestCase
tnlApproximationError< MeshType,
ExactOperator,
ApproximateOperator,
- TestFunction >
+ TestFunction,
+ ApproximationMethod >
::getError( mesh,
exactOperator,
approximateOperator,
diff --git a/tests/unit-tests/tnlApproximationError.h b/tests/unit-tests/tnlApproximationError.h
index 47ad750ceb..fed513b22c 100644
--- a/tests/unit-tests/tnlApproximationError.h
+++ b/tests/unit-tests/tnlApproximationError.h
@@ -19,12 +19,43 @@
#define TNLAPPROXIMATIONERROR_H_
#include <mesh/tnlGrid.h>
+#include <functions/tnlConstantFunction.h>
+#include <operators/tnlDirichletBoundaryConditions.h>
+
+class tnlExplicitApproximation
+{
+ public:
+
+ static tnlString getType()
+ {
+ return tnlString( "tnlExplicitApproximation" );
+ };
+};
+
+class tnlImplicitApproximation
+{
+ public:
+
+ static tnlString getType()
+ {
+ return tnlString( "tnlImplicitApproximation" );
+ };
+};
template< typename Mesh,
typename ExactOperator,
typename ApproximateOperator,
- typename Function >
+ typename Function,
+ typename ApproximationMethod >
class tnlApproximationError
+{
+};
+
+template< typename Mesh,
+ typename ExactOperator,
+ typename ApproximateOperator,
+ typename Function >
+class tnlApproximationError< Mesh, ExactOperator, ApproximateOperator, Function, tnlExplicitApproximation >
{
public:
@@ -33,6 +64,8 @@ class tnlApproximationError
typedef typename MeshType::DeviceType DeviceType;
typedef typename MeshType::IndexType IndexType;
typedef typename MeshType::VertexType VertexType;
+ typedef tnlConstantFunction< MeshType::Dimensions, RealType > ConstantFunctionType;
+ typedef tnlDirichletBoundaryConditions< MeshType, ConstantFunctionType > BoundaryConditionsType;
static void getError( const Mesh& mesh,
const ExactOperator& exactOperator,
@@ -50,7 +83,7 @@ template< int Dimensions,
typename ExactOperator,
typename ApproximateOperator,
typename Function >
-class tnlApproximationError< tnlGrid< Dimensions, Real, Device, Index >, ExactOperator, ApproximateOperator, Function >
+class tnlApproximationError< tnlGrid< Dimensions, Real, Device, Index >, ExactOperator, ApproximateOperator, Function, tnlExplicitApproximation >
{
public:
@@ -70,7 +103,30 @@ class tnlApproximationError< tnlGrid< Dimensions, Real, Device, Index >, ExactOp
RealType& maxErr );
};
+template< typename Mesh,
+ typename ExactOperator,
+ typename ApproximateOperator,
+ typename Function >
+class tnlApproximationError< Mesh, ExactOperator, ApproximateOperator, Function, tnlImplicitApproximation >
+{
+ public:
+
+ typedef typename ApproximateOperator::RealType RealType;
+ typedef Mesh MeshType;
+ typedef typename MeshType::DeviceType DeviceType;
+ typedef typename MeshType::IndexType IndexType;
+ typedef typename MeshType::VertexType VertexType;
+ typedef tnlConstantFunction< MeshType::Dimensions, RealType > ConstantFunctionType;
+ typedef tnlDirichletBoundaryConditions< MeshType, ConstantFunctionType > BoundaryConditionsType;
+ static void getError( const Mesh& mesh,
+ const ExactOperator& exactOperator,
+ const ApproximateOperator& approximateOperator,
+ const Function& function,
+ RealType& l1Err,
+ RealType& l2Err,
+ RealType& maxErr );
+};
#include "tnlApproximationError_impl.h"
diff --git a/tests/unit-tests/tnlApproximationError_impl.h b/tests/unit-tests/tnlApproximationError_impl.h
index f1c8fc1ab8..4b41d873a3 100644
--- a/tests/unit-tests/tnlApproximationError_impl.h
+++ b/tests/unit-tests/tnlApproximationError_impl.h
@@ -21,13 +21,16 @@
#include <mesh/tnlTraversal.h>
#include <core/vectors/tnlVector.h>
#include <functions/tnlFunctionDiscretizer.h>
+#include <matrices/tnlCSRMatrix.h>
+#include <matrices/tnlMatrixSetter.h>
+#include <solvers/pde/tnlLinearSystemAssembler.h>
template< typename Mesh,
typename ExactOperator,
typename ApproximateOperator,
typename Function >
void
-tnlApproximationError< Mesh, ExactOperator, ApproximateOperator, Function >::
+tnlApproximationError< Mesh, ExactOperator, ApproximateOperator, Function, tnlExplicitApproximation >::
getError( const Mesh& mesh,
const ExactOperator& exactOperator,
const ApproximateOperator& approximateOperator,
@@ -77,7 +80,7 @@ template< int Dimensions,
typename ApproximateOperator,
typename Function >
void
-tnlApproximationError< tnlGrid< Dimensions, Real, Device, Index >, ExactOperator, ApproximateOperator, Function >::
+tnlApproximationError< tnlGrid< Dimensions, Real, Device, Index >, ExactOperator, ApproximateOperator, Function, tnlExplicitApproximation >::
getError( const MeshType& mesh,
const ExactOperator& exactOperator,
const ApproximateOperator& approximateOperator,
@@ -120,5 +123,72 @@ getError( const MeshType& mesh,
maxErr = mesh.getDifferenceAbsMax( exactData, approximateData );
}
+/****
+ * Implicit (matrix) approximation
+ */
+
+template< typename Mesh,
+ typename ExactOperator,
+ typename ApproximateOperator,
+ typename Function >
+void
+tnlApproximationError< Mesh, ExactOperator, ApproximateOperator, Function, tnlImplicitApproximation >::
+getError( const Mesh& mesh,
+ const ExactOperator& exactOperator,
+ const ApproximateOperator& approximateOperator,
+ const Function& function,
+ RealType& l1Err,
+ RealType& l2Err,
+ RealType& maxErr )
+{
+ typedef tnlVector< RealType, DeviceType, IndexType > Vector;
+ typedef tnlCSRMatrix< RealType, DeviceType, IndexType > MatrixType;
+ typedef typename MatrixType::RowLengthsVector RowLengthsVectorType;
+ Vector functionData, exactData, approximateData;
+ MatrixType matrix;
+ RowLengthsVectorType rowLengths;
+ BoundaryConditionsType boundaryConditions;
+ ConstantFunctionType zeroFunction;
+
+ const IndexType entities = mesh.getNumberOfCells();
+
+ if( ! functionData.setSize( entities ) ||
+ ! exactData.setSize( entities ) ||
+ ! approximateData.setSize( entities ) ||
+ ! rowLengths.setSize( entities ) )
+ return;
+
+ tnlFunctionDiscretizer< Mesh, Function, Vector >::template discretize< 0, 0, 0 >( mesh, function, functionData );
+
+ tnlMatrixSetter< MeshType, ApproximateOperator, BoundaryConditionsType, RowLengthsVectorType > matrixSetter;
+ matrixSetter.template getRowLengths< Mesh::Dimensions >( mesh,
+ approximateOperator,
+ boundaryConditions,
+ rowLengths );
+ matrix.setDimensions( entities, entities );
+ if( ! matrix.setRowLengths( rowLengths ) )
+ return;
+
+ tnlLinearSystemAssembler< Mesh, Vector, ApproximateOperator, BoundaryConditionsType, ConstantFunctionType, MatrixType > systemAssembler;
+ systemAssembler.template assembly< Mesh::Dimensions >( 0.0, // time
+ 0.0, // tau
+ mesh,
+ approximateOperator,
+ boundaryConditions,
+ zeroFunction,
+ functionData,
+ matrix,
+ approximateData // this has no meaning here
+ );
+
+ matrix.vectorProduct( functionData, approximateData );
+ for( IndexType i = 0; i < entities; i++ )
+ if( mesh.isBoundaryCell( i ) )
+ approximateData.setElement( i, exactData.getElement( i ) );
+
+ l1Err = mesh.getDifferenceLpNorm( exactData, approximateData, ( RealType ) 1.0 );
+ l2Err = mesh.getDifferenceLpNorm( exactData, approximateData, ( RealType ) 2.0 );
+ maxErr = mesh.getDifferenceAbsMax( exactData, approximateData );
+}
#endif /* TNLAPPROXIMATIONERROR_IMPL_H_ */
--
GitLab