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

Modification for experimental faster explicit heat equation.

parent 8775d239

src/Examples/heat-equation/tnl-heat-equation.h

+81 −3
Original line number Diff line number Diff line
@@ -16,6 +16,7 @@ 
#include <TNL/Solvers/FastBuildConfigTag.h>

#include <TNL/Solvers/BuildConfigTags.h>

#include <TNL/Operators/diffusion/LinearDiffusion.h>

#include <TNL/Operators/diffusion/DirichletLinearDiffusion.h>

#include <TNL/Operators/DirichletBoundaryConditions.h>

#include <TNL/Operators/NeumannBoundaryConditions.h>

#include <TNL/Functions/Analytic/Constant.h>
@@ -24,12 +25,88 @@ 
#include <TNL/Meshes/Grid.h>

#include "HeatEquationBuildConfigTag.h"



#include <TNL/Operators/Operator.h>



using namespace TNL;

using namespace TNL::Problems;



//typedef Solvers::DefaultBuildConfigTag BuildConfig;

typedef Solvers::FastBuildConfigTag BuildConfig;



template< typename Mesh,

          int MeshEntitiesDimension = Mesh::getMeshDimension(),

          typename Real = typename Mesh::RealType,

          typename Index = typename Mesh::GlobalIndexType >

class DummyBC

: public Operators::Operator< Mesh,

                   Functions::MeshBoundaryDomain,

                   MeshEntitiesDimension,

                   MeshEntitiesDimension,

                   Real,

                   Index >

{

   public:



      typedef Mesh MeshType;

      typedef Real RealType;

      typedef typename MeshType::DeviceType DeviceType;

      typedef Index IndexType;



      typedef Pointers::SharedPointer<  Mesh > MeshPointer;

      typedef Containers::Vector< RealType, DeviceType, IndexType> DofVectorType;

      typedef typename MeshType::PointType PointType;



      static constexpr int getMeshDimension() { return MeshType::getMeshDimension(); }



      static void configSetup( Config::ConfigDescription& config,

                               const String& prefix = "" ){};



      bool setup( const MeshPointer& meshPointer,

                  const Config::ParameterContainer& parameters,

                  const String& prefix = "" )

      {

         return true;

      }



      template< typename EntityType,

                typename MeshFunction >

      __cuda_callable__

      const RealType operator()( const MeshFunction& u,

                                 const EntityType& entity,

                                 const RealType& time = 0 ) const

      {

         //static_assert( EntityType::getDimension() == MeshEntitiesDimension, "Wrong mesh entity dimension." );

         return 0.0;

      }



      template< typename EntityType >

      __cuda_callable__

      IndexType getLinearSystemRowLength( const MeshType& mesh,

                                          const IndexType& index,

                                          const EntityType& entity ) const

      {

         return 1;

      }



      template< typename PreimageFunction,

                typename MeshEntity,

                typename Matrix,

                typename Vector >

      __cuda_callable__

      void setMatrixElements( const PreimageFunction& u,

                              const MeshEntity& entity,

                              const RealType& time,

                              const RealType& tau,

                              Matrix& matrix,

                              Vector& b ) const

      {

         typename Matrix::MatrixRow matrixRow = matrix.getRow( entity.getIndex() );

         const IndexType& index = entity.getIndex();

         matrixRow.setElement( 0, index, 1.0 );

         b[ index ] = 0.0;

      }

};



template< typename MeshConfig >

class heatEquationConfig

{
@@ -73,8 +150,9 @@ class heatEquationSetter 
   static bool run( const Config::ParameterContainer& parameters )

   {

      enum { Dimension = MeshType::getMeshDimension() };

      typedef Operators::LinearDiffusion< MeshType, Real, Index > ApproximateOperator;

      typedef Functions::Analytic::Constant< Dimension, Real > RightHandSide;

      using ApproximateOperator = Operators::LinearDiffusion< MeshType, Real, Index >;

      using DirichletApproximateOperator = Operators::DirichletLinearDiffusion< MeshType, Real, Index >;

      using RightHandSide = Functions::Analytic::Constant< Dimension, Real >;



      String boundaryConditionsType = parameters.getParameter< String >( "boundary-conditions-type" );

      if( parameters.checkParameter( "boundary-conditions-constant" ) )
@@ -83,7 +161,7 @@ class heatEquationSetter 
         if( boundaryConditionsType == "dirichlet" )

         {

            typedef Operators::DirichletBoundaryConditions< MeshType, Constant > BoundaryConditions;

            typedef HeatEquationProblem< MeshType, BoundaryConditions, RightHandSide,CommunicatorType, ApproximateOperator > Problem;

            typedef HeatEquationProblem< MeshType, DummyBC< MeshType >, RightHandSide,CommunicatorType, DirichletApproximateOperator > Problem;

            SolverStarter solverStarter;

            return solverStarter.template run< Problem >( parameters );

         }

src/TNL/Operators/diffusion/DirichletLinearDiffusion.h

0 → 100644
+185 −0
Original line number Diff line number Diff line 
/***************************************************************************

                          DirichletLinearDiffusion.h  -  description

                             -------------------

    begin                : Aug 8, 2014

    copyright            : (C) 2014 by Tomas Oberhuber

    email                : tomas.oberhuber@fjfi.cvut.cz

 ***************************************************************************/



/* See Copyright Notice in tnl/Copyright */



/***

 * Authors:

 * Oberhuber Tomas, tomas.oberhuber@fjfi.cvut.cz

 * Szekely Ondrej, ondra.szekely@gmail.com

 */



#pragma once



#include <TNL/Containers/Vector.h>

#include <TNL/Functions/MeshFunction.h>

#include <TNL/Meshes/Grid.h>

#include <TNL/Operators/Operator.h>

#include <TNL/Operators/diffusion/ExactLinearDiffusion.h>



namespace TNL {

namespace Operators {



template< typename Mesh,

          typename Real = typename Mesh::RealType,

          typename Index = typename Mesh::GlobalIndexType >

class DirichletLinearDiffusion

{



};





template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

class DirichletLinearDiffusion< Meshes::Grid< 1,MeshReal, Device, MeshIndex >, Real, Index >

: public Operator< Meshes::Grid< 1, MeshReal, Device, MeshIndex >,

                      Functions::MeshInteriorDomain, 1, 1, Real, Index >

{

   public:



      typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType;

      typedef typename MeshType::CoordinatesType CoordinatesType;

      typedef Real RealType;

      typedef Device DeviceType;

      typedef Index IndexType;

      typedef ExactLinearDiffusion< 1 > ExactOperatorType;



      static const int Dimension = MeshType::getMeshDimension();



      static constexpr int getMeshDimension() { return Dimension; }



      template< typename PreimageFunction,

                typename MeshEntity >

      __cuda_callable__

      inline Real operator()( const PreimageFunction& u,

                              const MeshEntity& entity,

                              const RealType& time = 0.0 ) const;



      template< typename MeshEntity >

      __cuda_callable__

      inline Index getLinearSystemRowLength( const MeshType& mesh,

                                             const IndexType& index,

                                             const MeshEntity& entity ) const;



      template< typename PreimageFunction,

                typename MeshEntity,

                typename Matrix,

                typename Vector >

      __cuda_callable__

      inline void setMatrixElements( const PreimageFunction& u,

                                     const MeshEntity& entity,

                                     const RealType& time,

                                     const RealType& tau,

                                     Matrix& matrix,

                                     Vector& b ) const;

};





template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

class DirichletLinearDiffusion< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >

: public Operator< Meshes::Grid< 2, MeshReal, Device, MeshIndex >,

                      Functions::MeshInteriorDomain, 2, 2, Real, Index >

{

   public:



      typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType;

      typedef typename MeshType::CoordinatesType CoordinatesType;

      typedef Real RealType;

      typedef Device DeviceType;

      typedef Index IndexType;

      typedef ExactLinearDiffusion< 2 > ExactOperatorType;



      static const int Dimension = MeshType::getMeshDimension();



      static constexpr int getMeshDimension() { return Dimension; }



      template< typename PreimageFunction, typename EntityType >

      __cuda_callable__

      inline Real operator()( const PreimageFunction& u,

                              const EntityType& entity,

                              const Real& time = 0.0 ) const;



      template< typename EntityType >

      __cuda_callable__

      inline Index getLinearSystemRowLength( const MeshType& mesh,

                                             const IndexType& index,

                                             const EntityType& entity ) const;



      template< typename PreimageFunction,

                typename MeshEntity,

                typename Matrix,

                typename Vector >

      __cuda_callable__

      inline void setMatrixElements( const PreimageFunction& u,

                                     const MeshEntity& entity,

                                     const RealType& time,

                                     const RealType& tau,

                                     Matrix& matrix,

                                     Vector& b ) const;

};





template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

class DirichletLinearDiffusion< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >

: public Operator< Meshes::Grid< 3, MeshReal, Device, MeshIndex >,

                      Functions::MeshInteriorDomain, 3, 3, Real, Index >

{

   public:



      typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;

      typedef typename MeshType::CoordinatesType CoordinatesType;

      typedef Real RealType;

      typedef Device DeviceType;

      typedef Index IndexType;

      typedef ExactLinearDiffusion< 3 > ExactOperatorType;



      static const int Dimension = MeshType::getMeshDimension();



      static constexpr int getMeshDimension() { return Dimension; }



      template< typename PreimageFunction,

                typename EntityType >

      __cuda_callable__

      inline Real operator()( const PreimageFunction& u,

                              const EntityType& entity,

                              const Real& time = 0.0 ) const;



      template< typename EntityType >

      __cuda_callable__

      inline Index getLinearSystemRowLength( const MeshType& mesh,

                                             const IndexType& index,

                                             const EntityType& entity ) const;



      template< typename PreimageFunction,

                typename MeshEntity,

                typename Matrix,

                typename Vector >

      __cuda_callable__

      inline void setMatrixElements( const PreimageFunction& u,

                                     const MeshEntity& entity,

                                     const RealType& time,

                                     const RealType& tau,

                                     Matrix& matrix,

                                     Vector& b ) const;

};



} // namespace Operators

} // namespace TNL



#include <TNL/Operators/diffusion/DirichletLinearDiffusion.hpp>

src/TNL/Operators/diffusion/DirichletLinearDiffusion.hpp

0 → 100644
+261 −0
Original line number Diff line number Diff line 
/***************************************************************************

                          DirichletLinearDiffusion_impl.h  -  description

                             -------------------

    begin                : Aug 8, 2014

    copyright            : (C) 2014 by Tomas Oberhuber

    email                : tomas.oberhuber@fjfi.cvut.cz

 ***************************************************************************/



/* See Copyright Notice in tnl/Copyright */



/***

 * Authors:

 * Oberhuber Tomas, tomas.oberhuber@fjfi.cvut.cz

 * Szekely Ondrej, ondra.szekely@gmail.com

 */



#pragma once



#include <TNL/Operators/diffusion/DirichletLinearDiffusion.h>

#include <TNL/Meshes/Grid.h>



namespace TNL {

namespace Operators {



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

template< typename PreimageFunction,

          typename MeshEntity >

__cuda_callable__

inline

Real

DirichletLinearDiffusion< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >::

operator()( const PreimageFunction& u,

            const MeshEntity& entity,

            const Real& time ) const

{

   static_assert( MeshEntity::getEntityDimension() == 1, "Wrong mesh entity dimensions." );

   static_assert( PreimageFunction::getEntitiesDimension() == 1, "Wrong preimage function" );

   const typename MeshEntity::template NeighborEntities< 1 >& neighborEntities = entity.getNeighborEntities();

   const RealType& hxSquareInverse = entity.getMesh().template getSpaceStepsProducts< - 2 >();

   return ( u[ neighborEntities.template getEntityIndex< -1 >() ]

            - 2.0 * u[ entity.getIndex() ]

            + u[ neighborEntities.template getEntityIndex< 1 >() ] ) * hxSquareInverse;

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

   template< typename MeshEntity >

__cuda_callable__

inline

Index

DirichletLinearDiffusion< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >::

getLinearSystemRowLength( const MeshType& mesh,

                          const IndexType& index,

                          const MeshEntity& entity ) const

{

   return 3;

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

   template< typename PreimageFunction,

             typename MeshEntity,

             typename Matrix,

             typename Vector >

__cuda_callable__

inline

void

DirichletLinearDiffusion< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >::

setMatrixElements( const PreimageFunction& u,

                   const MeshEntity& entity,

                   const RealType& time,

                   const RealType& tau,

                   Matrix& matrix,

                   Vector& b ) const

{

   static_assert( MeshEntity::getEntityDimension() == 1, "Wrong mesh entity dimensions." );

   static_assert( PreimageFunction::getEntitiesDimension() == 1, "Wrong preimage function" );

   const typename MeshEntity::template NeighborEntities< 1 >& neighborEntities = entity.getNeighborEntities();

   const IndexType& index = entity.getIndex();

   typename Matrix::MatrixRow matrixRow = matrix.getRow( index );

   const RealType lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2 >();

   matrixRow.setElement( 0, neighborEntities.template getEntityIndex< -1 >(),      - lambdaX );

   matrixRow.setElement( 1, index,                                              2.0 * lambdaX );

   matrixRow.setElement( 2, neighborEntities.template getEntityIndex< 1 >(),       - lambdaX );

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

   template< typename EntityType >

__cuda_callable__

inline

Index

DirichletLinearDiffusion< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >::

getLinearSystemRowLength( const MeshType& mesh,

                          const IndexType& index,

                          const EntityType& entity ) const

{

   return 5;

}





template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

template< typename PreimageFunction,

          typename EntityType >

__cuda_callable__

inline

Real

DirichletLinearDiffusion< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >::

operator()( const PreimageFunction& u,

            const EntityType& entity,

            const Real& time ) const

{

   static_assert( EntityType::getEntityDimension() == 2, "Wrong mesh entity dimensions." );

   static_assert( PreimageFunction::getEntitiesDimension() == 2, "Wrong preimage function" );

   const typename EntityType::template NeighborEntities< 2 >& neighborEntities = entity.getNeighborEntities();

   const RealType& hxSquareInverse = entity.getMesh().template getSpaceStepsProducts< -2, 0 >();

   const RealType& hySquareInverse = entity.getMesh().template getSpaceStepsProducts< 0, -2 >();

   return ( u[ neighborEntities.template getEntityIndex< -1,  0 >() ]

          + u[ neighborEntities.template getEntityIndex<  1,  0 >() ] ) * hxSquareInverse +

          ( u[ neighborEntities.template getEntityIndex<  0, -1 >() ]

          + u[ neighborEntities.template getEntityIndex<  0,  1 >() ] ) * hySquareInverse

          - 2.0 * u[ entity.getIndex() ] * ( hxSquareInverse + hySquareInverse );

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

   template< typename PreimageFunction,

             typename MeshEntity,

             typename Matrix,

             typename Vector >

__cuda_callable__

inline

void

DirichletLinearDiffusion< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >::

setMatrixElements( const PreimageFunction& u,

                   const MeshEntity& entity,

                   const RealType& time,

                   const RealType& tau,

                   Matrix& matrix,

                   Vector& b ) const

{

   static_assert( MeshEntity::getEntityDimension() == 2, "Wrong mesh entity dimensions." );

   static_assert( PreimageFunction::getEntitiesDimension() == 2, "Wrong preimage function" );

   const IndexType& index = entity.getIndex();

   typename Matrix::MatrixRow matrixRow = matrix.getRow( index );

   const RealType lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2, 0 >();

   const RealType lambdaY = tau * entity.getMesh().template getSpaceStepsProducts< 0, -2 >();

   const typename MeshEntity::template NeighborEntities< 2 >& neighborEntities = entity.getNeighborEntities();

   matrixRow.setElement( 0, neighborEntities.template getEntityIndex< 0, -1 >(), -lambdaY );

   matrixRow.setElement( 1, neighborEntities.template getEntityIndex< -1, 0 >(), -lambdaX );

   matrixRow.setElement( 2, index,                                                        2.0 * ( lambdaX + lambdaY ) );

   matrixRow.setElement( 3, neighborEntities.template getEntityIndex< 1, 0 >(),   -lambdaX );

   matrixRow.setElement( 4, neighborEntities.template getEntityIndex< 0, 1 >(),   -lambdaY );

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

template< typename PreimageFunction,

          typename EntityType >

__cuda_callable__

inline

Real

DirichletLinearDiffusion< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >::

operator()( const PreimageFunction& u,

            const EntityType& entity,

            const Real& time ) const

{

   static_assert( EntityType::getEntityDimension() == 3, "Wrong mesh entity dimensions." );

   static_assert( PreimageFunction::getEntitiesDimension() == 3, "Wrong preimage function" );

   const typename EntityType::template NeighborEntities< 3 >& neighborEntities = entity.getNeighborEntities();

   const RealType& hxSquareInverse = entity.getMesh().template getSpaceStepsProducts< -2,  0,  0 >();

   const RealType& hySquareInverse = entity.getMesh().template getSpaceStepsProducts<  0, -2,  0 >();

   const RealType& hzSquareInverse = entity.getMesh().template getSpaceStepsProducts<  0,  0, -2 >();

   return (   u[ neighborEntities.template getEntityIndex< -1,  0,  0 >() ]

            + u[ neighborEntities.template getEntityIndex<  1,  0,  0 >() ] ) * hxSquareInverse +

          (   u[ neighborEntities.template getEntityIndex<  0, -1,  0 >() ]

            + u[ neighborEntities.template getEntityIndex<  0,  1,  0 >() ] ) * hySquareInverse +

          (   u[ neighborEntities.template getEntityIndex<  0,  0, -1 >() ]

            + u[ neighborEntities.template getEntityIndex<  0,  0,  1 >() ] ) * hzSquareInverse

         - 2.0 * u[ entity.getIndex() ] * ( hxSquareInverse + hySquareInverse + hzSquareInverse );

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

   template< typename EntityType >

__cuda_callable__

inline

Index

DirichletLinearDiffusion< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >::

getLinearSystemRowLength( const MeshType& mesh,

                          const IndexType& index,

                          const EntityType& entity ) const

{

   return 7;

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

   template< typename PreimageFunction,

             typename MeshEntity,

             typename Matrix,

             typename Vector >

__cuda_callable__

inline

void

DirichletLinearDiffusion< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >::

setMatrixElements( const PreimageFunction& u,

                   const MeshEntity& entity,

                   const RealType& time,

                   const RealType& tau,

                   Matrix& matrix,

                   Vector& b ) const

{

   static_assert( MeshEntity::getEntityDimension() == 3, "Wrong mesh entity dimensions." );

   static_assert( PreimageFunction::getEntitiesDimension() == 3, "Wrong preimage function" );

   const typename MeshEntity::template NeighborEntities< 3 >& neighborEntities = entity.getNeighborEntities();

   const IndexType& index = entity.getIndex();

   typename Matrix::MatrixRow matrixRow = matrix.getRow( index );

   const RealType lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2, 0, 0 >();

   const RealType lambdaY = tau * entity.getMesh().template getSpaceStepsProducts< 0, -2, 0 >();

   const RealType lambdaZ = tau * entity.getMesh().template getSpaceStepsProducts< 0, 0, -2 >();

   matrixRow.setElement( 0, neighborEntities.template getEntityIndex< 0, 0, -1 >(), -lambdaZ );

   matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 0, -1, 0 >(), -lambdaY );

   matrixRow.setElement( 2, neighborEntities.template getEntityIndex< -1, 0, 0 >(), -lambdaX );

   matrixRow.setElement( 3, index,                             2.0 * ( lambdaX + lambdaY + lambdaZ ) );

   matrixRow.setElement( 4, neighborEntities.template getEntityIndex< 1, 0, 0 >(),   -lambdaX );

   matrixRow.setElement( 5, neighborEntities.template getEntityIndex< 0, 1, 0 >(),   -lambdaY );

   matrixRow.setElement( 6, neighborEntities.template getEntityIndex< 0, 0, 1 >(),   -lambdaZ );

}



} // namespace Operators

} // namespace TNL

src/TNL/Problems/HeatEquationProblem_impl.h

+7 −7
Original line number Diff line number Diff line
@@ -265,8 +265,8 @@ HeatEquationProblem< Mesh, BoundaryCondition, RightHandSide, Communicator, Diffe 
applyBoundaryConditions( const RealType& time,

                         DofVectorPointer& uDofs )

{

   this->bindDofs( uDofs );

   this->explicitUpdater.template applyBoundaryConditions< typename Mesh::Cell >( this->getMesh(), time, this->uPointer );

   //this->bindDofs( uDofs );

   //this->explicitUpdater.template applyBoundaryConditions< typename Mesh::Cell >( this->getMesh(), time, this->uPointer );

}



template< typename Mesh,

src/TNL/Solvers/FastBuildConfigTag.h

+1 −1
Original line number Diff line number Diff line
@@ -26,7 +26,7 @@ class FastBuildConfigTag 
/****

 * Turn off support for float and long double.

 */

template<> struct ConfigTagReal< FastBuildConfigTag, float > { enum { enabled = false }; };

template<> struct ConfigTagReal< FastBuildConfigTag, float > { enum { enabled = true }; };

template<> struct ConfigTagReal< FastBuildConfigTag, long double > { enum { enabled = false }; };



/****

src/TNL/Problems/HeatEquationProblem.h

+2 −2

File changed.

Contains only whitespace changes.

src/TNL/Solvers/PDE/ExplicitUpdater.h

+17 −17

File changed.

Contains only whitespace changes.