Commit 5634a143 authored by Tomáš Oberhuber's avatar Tomáš Oberhuber
Browse files

Merge branch 'diffusion'

parents c6a62119 15f0a73c

examples/heat-equation/heatEquationSolver.h

+18 −5
Original line number Diff line number Diff line
@@ -25,6 +25,8 @@ 
#include <core/vectors/tnlVector.h>

#include <core/vectors/tnlSharedVector.h>

#include <solvers/pde/tnlExplicitUpdater.h>

#include <solvers/pde/tnlLinearSystemAssembler.h>

#include <matrices/tnlCSRMatrix.h>

#include "heatEquationSolver.h"
@@ -41,6 +43,8 @@ class heatEquationSolver 
   typedef typename DifferentialOperator::IndexType IndexType;

   typedef Mesh MeshType;

   typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;

   typedef tnlCSRMatrix< RealType, DeviceType, IndexType > MatrixType;

   typedef typename MatrixType::RowLengthsVector RowLengthsVectorType;



   static tnlString getTypeStatic();
@@ -49,10 +53,14 @@ class heatEquationSolver 
   void writeProlog( tnlLogger& logger,

                     const tnlParameterContainer& parameters ) const;



   bool init( const tnlParameterContainer& parameters );

   bool setup( const tnlParameterContainer& parameters );



   bool setInitialCondition( const tnlParameterContainer& parameters,

                             const MeshType& mesh );

                             const MeshType& mesh,

                             DofVectorType& dofs );



   bool setupLinearSystem( const MeshType& mesh,

                           MatrixType& matrix );



   bool makeSnapshot( const RealType& time,

                      const IndexType& step,
@@ -68,20 +76,25 @@ class heatEquationSolver 
   void bindAuxiliaryDofs( const MeshType& mesh,

                           DofVectorType& auxiliaryDofs );



   void GetExplicitRHS( const RealType& time,

   void getExplicitRHS( const RealType& time,

                        const RealType& tau,

                        const MeshType& mesh,

                        DofVectorType& _u,

                        DofVectorType& _fu );



   void assemblyLinearSystem( const RealType& time,

                              const RealType& tau,

                              const MeshType& mesh,

                              DofVectorType& u,

                              MatrixType& matrix,

                              DofVectorType& rightHandSide );



   tnlSolverMonitor< RealType, IndexType >* getSolverMonitor();

   

   protected:



   tnlSharedVector< RealType, DeviceType, IndexType > solution;



   tnlExplicitUpdater< Mesh, DofVectorType, DifferentialOperator, BoundaryCondition, RightHandSide  > explicitUpdater;



   DifferentialOperator differentialOperator;



   BoundaryCondition boundaryCondition;

examples/heat-equation/heatEquationSolver_impl.h

+121 −46
Original line number Diff line number Diff line
@@ -19,47 +19,65 @@ 
#define HEATEQUATIONSOLVER_IMPL_H_



#include <core/mfilename.h>

#include <matrices/tnlMatrixSetter.h>

#include "heatEquationSolver.h"





template< typename Mesh, typename Diffusion, typename BoundaryCondition, typename RightHandSide >

tnlString heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide > 

::getTypeStatic()

template< typename Mesh,

          typename DifferentialOperator,

          typename BoundaryCondition,

          typename RightHandSide >

tnlString

heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide >::

getTypeStatic()

{

   return tnlString( "heatEquationSolver< " ) + Mesh :: getTypeStatic() + " >";

}



template< typename Mesh, typename Diffusion, typename BoundaryCondition, typename RightHandSide >

tnlString heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide >

:: getPrologHeader() const

template< typename Mesh,

          typename DifferentialOperator,

          typename BoundaryCondition,

          typename RightHandSide >

tnlString

heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide >::

getPrologHeader() const

{

   return tnlString( "Heat equation" );

}



template< typename Mesh, typename Diffusion, typename BoundaryCondition, typename RightHandSide >

void heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide >

:: writeProlog( tnlLogger& logger, const tnlParameterContainer& parameters ) const

template< typename Mesh,

          typename DifferentialOperator,

          typename BoundaryCondition,

          typename RightHandSide >

void

heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide >::

writeProlog( tnlLogger& logger, const tnlParameterContainer& parameters ) const

{

   //logger. WriteParameter< tnlString >( "Problem name:", "problem-name", parameters );

   //logger. WriteParameter< int >( "Simple parameter:", 1 );

}



template< typename Mesh, typename Diffusion, typename BoundaryCondition, typename RightHandSide >

bool heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide >

:: init( const tnlParameterContainer& parameters )

template< typename Mesh,

          typename DifferentialOperator,

          typename BoundaryCondition,

          typename RightHandSide >

bool

heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide >::

setup( const tnlParameterContainer& parameters )

{

   if( ! boundaryCondition.init( parameters ) ||

       ! rightHandSide.init( parameters ) )

   if( ! boundaryCondition.setup( parameters ) ||

       ! rightHandSide.setup( parameters ) )

      return false;

   return true;

}



template< typename Mesh,

          typename Diffusion,

          typename DifferentialOperator,

          typename BoundaryCondition,

          typename RightHandSide >

typename heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide >::IndexType 

   heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide >::getDofs( const Mesh& mesh ) const

typename heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide >::IndexType

heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide >::

getDofs( const Mesh& mesh ) const

{

   /****

    * Set-up DOFs and supporting grid functions
@@ -68,23 +86,25 @@ typename heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide >::I 
}



template< typename Mesh,

          typename Diffusion,

          typename DifferentialOperator,

          typename BoundaryCondition,

          typename RightHandSide >

typename heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide >::IndexType

   heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide >::getAuxiliaryDofs( const Mesh& mesh ) const

typename heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide >::IndexType

heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide >::

getAuxiliaryDofs( const Mesh& mesh ) const

{

   /****

    * Set-up DOFs and supporting grid functions which will not appear in the discrete solver

    */

   return 0;

}



template< typename Mesh,

          typename Diffusion,

          typename DifferentialOperator,

          typename BoundaryCondition,

          typename RightHandSide >

void

heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide >::

heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide >::

bindDofs( const MeshType& mesh,

          DofVectorType& dofVector )

{
@@ -93,39 +113,57 @@ bindDofs( const MeshType& mesh, 
}



template< typename Mesh,

          typename Diffusion,

          typename DifferentialOperator,

          typename BoundaryCondition,

          typename RightHandSide >

void

heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide >::

heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide >::

bindAuxiliaryDofs( const MeshType& mesh,

                   DofVectorType& auxiliaryDofVector )

{

}





template< typename Mesh, typename Diffusion, typename BoundaryCondition, typename RightHandSide >

bool heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide >

template< typename Mesh, typename DifferentialOperator, typename BoundaryCondition, typename RightHandSide >

bool heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide >

:: setInitialCondition( const tnlParameterContainer& parameters,

                        const MeshType& mesh )

                        const MeshType& mesh,

                        DofVectorType& dofs )

{

   /*const tnlString& initialConditionFile = parameters.GetParameter< tnlString >( "initial-condition" );

   this->bindDofs( mesh, dofs );

   const tnlString& initialConditionFile = parameters.GetParameter< tnlString >( "initial-condition" );

   if( ! this->solution.load( initialConditionFile ) )

   {

      cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." << endl;

      return false;

   }*/

   }

   return true;

}



template< typename Mesh,

          typename Diffusion,

          typename DifferentialOperator,

          typename BoundaryCondition,

          typename RightHandSide >

bool

heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide >::

setupLinearSystem( const MeshType& mesh,

                   MatrixType& matrix )

{

   RowLengthsVectorType rowLengths;

   tnlMatrixSetter< MeshType, DifferentialOperator, BoundaryCondition, RowLengthsVectorType > matrixSetter;

   matrixSetter.template getRowLengths< Mesh::Dimensions >( mesh,

                                                            differentialOperator,

                                                            boundaryCondition,

                                                            rowLengths );

   matrix.setRowLengths( rowLengths );

}



template< typename Mesh,

          typename DifferentialOperator,

          typename BoundaryCondition,

          typename RightHandSide >

bool heatEquationSolver< Mesh,

                         Diffusion,

                         BoundaryCondition,

                         RightHandSide >::

bool

heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide >::

makeSnapshot( const RealType& time,

              const IndexType& step,

              const MeshType& mesh )
@@ -139,9 +177,13 @@ makeSnapshot( const RealType& time, 
   return true;

}



template< typename Mesh, typename Diffusion, typename BoundaryCondition, typename RightHandSide >

void heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide > 

:: GetExplicitRHS( const RealType& time,

template< typename Mesh,

          typename DifferentialOperator,

          typename BoundaryCondition,

          typename RightHandSide >

void

heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide >::

getExplicitRHS( const RealType& time,

                const RealType& tau,

                const Mesh& mesh,

                DofVectorType& _u,
@@ -157,6 +199,7 @@ void heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide > 
    */



   this->bindDofs( mesh, _u );

   tnlExplicitUpdater< Mesh, DofVectorType, DifferentialOperator, BoundaryCondition, RightHandSide > explicitUpdater;

   explicitUpdater.template update< Mesh::Dimensions >( time,

                                                        mesh,

                                                        this->differentialOperator,
@@ -164,13 +207,45 @@ void heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide > 
                                                        this->rightHandSide,

                                                        _u,

                                                        _fu );

   //_u.save( "u.tnl" );

   //_fu.save( "fu.tnl" );

   //getchar();

}



template< typename Mesh, typename Diffusion, typename BoundaryCondition, typename RightHandSide >

tnlSolverMonitor< typename heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide >::RealType,

                  typename heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide >::IndexType >*

                  heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide > 

::  getSolverMonitor()

template< typename Mesh,

          typename DifferentialOperator,

          typename BoundaryCondition,

          typename RightHandSide >

void

heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide >::

assemblyLinearSystem( const RealType& time,

                      const RealType& tau,

                      const MeshType& mesh,

                      DofVectorType& u,

                      MatrixType& matrix,

                      DofVectorType& b )

{

   tnlLinearSystemAssembler< Mesh, DofVectorType, DifferentialOperator, BoundaryCondition, RightHandSide, MatrixType > systemAssembler;

   /*systemAssembler.template assembly< Mesh::Dimensions >( time,

                                                          tau,

                                                          mesh,

                                                          this->differentialOperator,

                                                          this->boundaryCondition,

                                                          this->rightHandSide,

                                                          u,

                                                          matrix,

                                                          b );

    */

}



template< typename Mesh,

          typename DifferentialOperator,

          typename BoundaryCondition,

          typename RightHandSide >

tnlSolverMonitor< typename heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide >::RealType,

                  typename heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide >::IndexType >*

heatEquationSolver< Mesh,DifferentialOperator,BoundaryCondition,RightHandSide >::

getSolverMonitor()

{

   return 0;

}

examples/heat-equation/tnl-heat-equation.cpp

+1 −2
Original line number Diff line number Diff line
@@ -32,8 +32,7 @@ class heatEquationConfig 
   public:

      static void configSetup( tnlConfigDescription& config )

      {

         config.addDelimiter( "Heat equation settings:" );

         config.addDelimiter( "Tests setting::" );

         //config.addDelimiter( "Heat equation settings:" );

      }

};

examples/heat-equation/tnl-run-heat-equation-eoc-test

+11 −8
Original line number Diff line number Diff line
@@ -5,17 +5,17 @@ dimensions="1D 2D 3D" 
sizes1D="16 32 64 128 256 512"

sizes2D="16 32 64 128 256 512"

sizes3D="16 32 64 128"

testFunctions="exp-bump"

testFunctions="sin-wave"

snapshotPeriod=0.1

finalTime=1

timeDependence="linear"

finalTime=1.5

timeDependence="cosine"

solverName="tnl-heat-equation-eoc-test"



setupTestFunction()

{

   testFunction=$1

   if test x${testFunction} = "xexp-bump";

   then

#   if test x${testFunction} = "xexp-bump";

#   then

      origin=-1.0

      proportions=2.0

      amplitude=1.0
@@ -32,7 +32,7 @@ setupTestFunction() 
      phaseY=0.0

      phaseZ=0.0

      sigma=0.5

   fi

#   fi

}



setupGrid()
@@ -54,7 +54,7 @@ setupGrid() 
setInitialCondition()

{

   testFunction=$1

   tnl-init --function ${testFunction} \

   tnl-init --test-function ${testFunction} \

            --output-file exact-u.tnl \

            --amplitude ${amplitude} \

            --wave-length ${waveLength} \
@@ -82,7 +82,9 @@ solve() 
   ${solverName} --mesh mesh.tnl \

                 --initial-condition exact-u-00000.tnl \

                 --time-discretisation ${timeDiscretisation} \

                 --time-step 0.001 \

                 --discrete-solver ${discreteSolver} \

                 --merson-adaptivity 1.0e-5 \

                 --test-function ${testFunction}\

                 --amplitude ${amplitude} \

                 --wave-length ${waveLength} \
@@ -110,7 +112,8 @@ computeError() 
            --input-files exact-u-*.tnl u-*.tnl \

            --mode halves \

            --snapshot-period ${snapshotPeriod} \

            --output-file errors.txt 

            --output-file errors.txt \

            --write-difference yes

}



runTest()

examples/heat-equation/tnlHeatEquationEocRhs.h

+4 −3
Original line number Diff line number Diff line
@@ -27,9 +27,9 @@ class tnlHeatEquationEocRhs 
      typedef ExactOperator ExactOperatorType;

      typedef TestFunction TestFunctionType;



      bool init( const tnlParameterContainer& parameters )

      bool setup( const tnlParameterContainer& parameters )

      {

         if( ! testFunction.init( parameters ) )

         if( ! testFunction.setup( parameters ) )

            return false;

         return true;

      };
@@ -39,7 +39,8 @@ class tnlHeatEquationEocRhs 
      Real getValue( const Vertex& vertex,

                     const Real& time )

      {

         return testFunction.getTimeDerivative( vertex, time ) - exactOperator.getValue( testFunction, vertex, time );

         return testFunction.getTimeDerivative( vertex, time )

                - exactOperator.getValue( testFunction, vertex, time );

      };



   protected: