Commit e3943328 authored by Ondrej Szekely's avatar Ondrej Szekely
Browse files

Adding Finite Volume method for 2D mesh

parent f53c548b

examples/mean-curvature-flow/tnl-mean-curvature-flow-eoc.h

+5 −3
Original line number Diff line number Diff line
@@ -29,8 +29,10 @@ 
#include <operators/diffusion/tnlExactNonlinearDiffusion.h>

#include <operators/diffusion/tnlNonlinearDiffusion.h>

#include <operators/operator-Q/tnlOneSideDiffOperatorQ.h>

#include <operators/operator-Q/tnlFiniteVolumeOperatorQ.h>

#include <operators/diffusion/tnlExactNonlinearDiffusion.h>

#include <operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator.h>

#include <operators/diffusion/nonlinear-diffusion-operators/tnlFiniteVolumeNonlinearOperator.h>

#include <operators/operator-Q/tnlExactOperatorQ.h>



//typedef tnlDefaultConfigTag BuildConfig;
@@ -45,7 +47,7 @@ class meanCurvatureFlowEocConfig 
         config.addDelimiter( "Mean Curvature Flow EOC settings:" );

         config.addEntry< double >( "eps", "This sets a eps in operator Q.", 1.0 );

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

         tnlTestFunction< 2, double >::configSetup( config );

         tnlTestFunction< 3, double >::configSetup( config );

      }

};
@@ -68,8 +70,8 @@ class meanCurvatureFlowEocSetter 
   static bool run( const tnlParameterContainer& parameters )

   {

      enum { Dimensions = MeshType::Dimensions };

      typedef tnlOneSideDiffOperatorQ<MeshType, Real, Index, 0> OperatorQ;

      typedef tnlOneSideDiffNonlinearOperator<MeshType, OperatorQ, Real, Index > NonlinearOperator;

      typedef tnlFiniteVolumeOperatorQ<MeshType, Real, Index, 0> OperatorQ;

      typedef tnlFiniteVolumeNonlinearOperator<MeshType, OperatorQ, Real, Index > NonlinearOperator;

      typedef tnlNonlinearDiffusion< MeshType, NonlinearOperator, Real, Index > ApproximateOperator;

      typedef tnlExactNonlinearDiffusion< tnlExactOperatorQ<Dimensions>, Dimensions > ExactOperator;

      typedef tnlTestFunction< MeshType::Dimensions, Real, Device > TestFunction;

src/operators/diffusion/nonlinear-diffusion-operators/CMakeLists.txt

+3 −1
Original line number Diff line number Diff line 
SET( headers tnlOneSideDiffNonlinearOperator.h

             tnlOneSideDiffNonlinearOperator_impl.h )

             tnlOneSideDiffNonlinearOperator_impl.h

	     tnlFiniteVolumeNonlinearOperator.h

             tnlFiniteVolumeNonlinearOperator_impl.h )



INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/operators/diffusion/nonlinear-diffusion-operators )

src/operators/diffusion/nonlinear-diffusion-operators/tnlFiniteVolumeNonlinearOperator.h

0 → 100644
+187 −0
Original line number Diff line number Diff line 
#ifndef TNLFINITEVOLUMENONLINEAROPERATOR_H

#define	TNLFINITEVOLUMENONLINEAROPERATOR_H



#include <core/vectors/tnlVector.h>

#include <mesh/tnlGrid.h>



template< typename Mesh,

          typename NonlinearDiffusionOperator,

	  typename OperatorQ,

          typename Real = typename Mesh::RealType,

          typename Index = typename Mesh::IndexType >

class tnlFiniteVolumeNonlinearOperator

{

 

};





template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index,

          typename OperatorQ >

class tnlFiniteVolumeNonlinearOperator< tnlGrid< 1,MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >

{

   public: 

   

   typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;

   typedef typename MeshType::CoordinatesType CoordinatesType;

   typedef Real RealType;

   typedef Device DeviceType;

   typedef Index IndexType;

   typedef OperatorQ OperatorQType;



   static tnlString getType();

   

   template< typename Vector >

#ifdef HAVE_CUDA

   __device__ __host__

#endif

   Real getValue( const MeshType& mesh,

                  const IndexType cellIndex,

                  const CoordinatesType& coordinates,

                  const Vector& u,

                  const RealType& time) const;

   

#ifdef HAVE_CUDA

   __device__ __host__

#endif

   Index getLinearSystemRowLength( const MeshType& mesh,

                                   const IndexType& index,

                                   const CoordinatesType& coordinates ) const;



   template< typename Vector, typename MatrixRow >

#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,

                               MatrixRow& matrixRow ) const;

   

   public:

   

   OperatorQ operatorQ;

};





template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index,

          typename OperatorQ >

class tnlFiniteVolumeNonlinearOperator< tnlGrid< 2, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >

{

   public: 

   

   typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;

   typedef typename MeshType::CoordinatesType CoordinatesType;

   typedef Real RealType;

   typedef Device DeviceType;

   typedef Index IndexType;

   typedef OperatorQ OperatorQType;

   



   static tnlString getType();

   

   template< typename Vector >

#ifdef HAVE_CUDA

   __device__ __host__

#endif

   Real getValue( const MeshType& mesh,

                  const IndexType cellIndex,

                  const CoordinatesType& coordinates,

                  const Vector& u,

                  const RealType& time) const;

   

#ifdef HAVE_CUDA

   __device__ __host__

#endif

   Index getLinearSystemRowLength( const MeshType& mesh,

                                   const IndexType& index,

                                   const CoordinatesType& coordinates ) const;



   template< typename Vector, typename MatrixRow >

#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,

                               MatrixRow& matrixRow ) const;

   

   public:

   

   OperatorQ operatorQ;

};





template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index,

          typename OperatorQ >

class tnlFiniteVolumeNonlinearOperator< tnlGrid< 3, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >

{

   public: 

   

   typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;

   typedef typename MeshType::CoordinatesType CoordinatesType;

   typedef Real RealType;

   typedef Device DeviceType;

   typedef Index IndexType;

   typedef OperatorQ OperatorQType;



   static tnlString getType();

   

   template< typename Vector >

#ifdef HAVE_CUDA

   __device__ __host__

#endif

   Real getValue( const MeshType& mesh,

                  const IndexType cellIndex,

                  const CoordinatesType& coordinates,

                  const Vector& u,

                  const RealType& time) const;

   

#ifdef HAVE_CUDA

   __device__ __host__

#endif

   Index getLinearSystemRowLength( const MeshType& mesh,

                                   const IndexType& index,

                                   const CoordinatesType& coordinates ) const;



   template< typename Vector, typename MatrixRow >

#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,

                               MatrixRow& matrixRow ) const;

   

   public:

   

   OperatorQ operatorQ;

};





#include "tnlFiniteVolumeNonlinearOperator_impl.h"





#endif	/* TNLFINITEVOLUMENONLINEAROPERATOR_H */

src/operators/diffusion/nonlinear-diffusion-operators/tnlFiniteVolumeNonlinearOperator_impl.h

0 → 100644
+307 −0
Original line number Diff line number Diff line 


#ifndef TNLFINITEVOLUMENONLINEAROPERATOR__IMPL_H

#define	TNLFINITEVOLUMENONLINEAROPERATOR__IMPL_H



#include "tnlFiniteVolumeNonlinearOperator.h"



#include <mesh/tnlGrid.h>



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index,

          typename OperatorQ >

tnlString

tnlFiniteVolumeNonlinearOperator< tnlGrid< 1, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::

getType()

{

   return tnlString( "tnlFiniteVolumeNonlinearOperator< " ) +

          MeshType::getType() + ", " +

          ::getType< Real >() + ", " +

          ::getType< Index >() + ", " +

	  OperatorQ::getType() + " >";

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index,

          typename OperatorQ >

template< typename Vector >

#ifdef HAVE_CUDA

__device__ __host__

#endif

Real

tnlFiniteVolumeNonlinearOperator< tnlGrid< 1, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::

getValue( const MeshType& mesh,

          const IndexType cellIndex,

          const CoordinatesType& coordinates,

          const Vector& u,

          const Real& time ) const

{

   return 0.0;

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index,

          typename OperatorQ >

#ifdef HAVE_CUDA

   __device__ __host__

#endif

Index

tnlFiniteVolumeNonlinearOperator< tnlGrid< 1, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::

getLinearSystemRowLength( const MeshType& mesh,

                          const IndexType& index,

                          const CoordinatesType& coordinates ) const

{

   return 1;

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index,

          typename OperatorQ >

template< typename Vector, typename MatrixRow >

#ifdef HAVE_CUDA

__device__ __host__

#endif

void

tnlFiniteVolumeNonlinearOperator< tnlGrid< 1, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::

updateLinearSystem( const RealType& time,

                    const RealType& tau,

                    const MeshType& mesh,

                    const IndexType& index,

                    const CoordinatesType& coordinates,

                    Vector& u,

                    Vector& b,

                    MatrixRow& matrixRow ) const

{

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index,

	  typename OperatorQ >

tnlString

tnlFiniteVolumeNonlinearOperator< tnlGrid< 2, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::

getType()

{

   return tnlString( "tnlFiniteVolumeNonlinearOperator< " ) +

          MeshType::getType() + ", " +

          ::getType< Real >() + ", " +

          ::getType< Index >() + ", " +

	  OperatorQ::getType() + " >";

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index,

          typename OperatorQ >

template< typename Vector >

#ifdef HAVE_CUDA

__device__ __host__

#endif

Real

tnlFiniteVolumeNonlinearOperator< tnlGrid< 2, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::

getValue( const MeshType& mesh,

          const IndexType cellIndex,

          const CoordinatesType& coordinates,

          const Vector& u,

          const Real& time ) const

{

   return operatorQ.getValue(mesh,cellIndex, coordinates, u, time)*( (u[ mesh.template getCellNextToCell< 1,0 >( cellIndex ) ] - u[ cellIndex ]) 

          * mesh.getHxSquareInverse() / operatorQ.getValue(mesh,cellIndex, coordinates, u, time, 1 )

          + ( u[ mesh.template getCellNextToCell< 0,1 >( cellIndex ) ] - u[ cellIndex ]) * mesh.getHySquareInverse()/

          operatorQ.getValue(mesh,cellIndex, coordinates, u, time, 0, 1 ) 

          - ( - u[ mesh.template getCellNextToCell< -1,0 >( cellIndex ) ]  + u[ cellIndex ]) 

          * mesh.getHxSquareInverse()/operatorQ.getValue(mesh,cellIndex, coordinates, u, time, -1)

          -( - u[ mesh.template getCellNextToCell< 0,-1 >( cellIndex ) ] + u[ cellIndex ]) * mesh.getHySquareInverse()

          /operatorQ.getValue(mesh,cellIndex,coordinates, u, time, 0, -1) );

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index,

          typename OperatorQ >

#ifdef HAVE_CUDA

   __device__ __host__

#endif

Index

tnlFiniteVolumeNonlinearOperator< tnlGrid< 2, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::

getLinearSystemRowLength( const MeshType& mesh,

                          const IndexType& index,

                          const CoordinatesType& coordinates ) const

{

   return 5;

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index,

          typename OperatorQ >

template< typename Vector, typename MatrixRow >

#ifdef HAVE_CUDA

__device__ __host__

#endif

void

tnlFiniteVolumeNonlinearOperator< tnlGrid< 2, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::

updateLinearSystem( const RealType& time,

                    const RealType& tau,

                    const MeshType& mesh,

                    const IndexType& index,

                    const CoordinatesType& coordinates,

                    Vector& u,

                    Vector& b,

                    MatrixRow& matrixRow ) const

{

   const RealType aCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time ) * mesh.getHySquareInverse() / 

                       operatorQ.getValue(mesh, index, coordinates, u, time, 0, -1 );

   const RealType bCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time ) * mesh.getHxSquareInverse() / 

                       operatorQ.getValue(mesh, index, coordinates, u, time, -1 );

   const RealType cCoef = tau * operatorQ.getValue(mesh, index, coordinates, u, time ) * ( mesh.getHxSquareInverse() / 

                       operatorQ.getValue(mesh, index, coordinates, u, time, 1 ) + mesh.getHySquareInverse() / 

                       operatorQ.getValue(mesh, index, coordinates, u, time, 0, 1 )

                       + mesh.getHxSquareInverse() / operatorQ.getValue(mesh, index, coordinates, u, time, -1 ) + 

                       mesh.getHySquareInverse() / operatorQ.getValue(mesh, index, coordinates, u, time, 0, -1 ) );

   const RealType dCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time ) * mesh.getHxSquareInverse() / 

                       operatorQ.getValue(mesh, index, coordinates, u, time, 1 );

   const RealType eCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time ) * mesh.getHySquareInverse() / 

                       operatorQ.getValue(mesh, index, coordinates, u, time, 0, 1 );

   matrixRow.setElement( 0, mesh.template getCellNextToCell< 0,-1 >( index ),     aCoef );

   matrixRow.setElement( 1, mesh.template getCellNextToCell< -1,0 >( index ),     bCoef );

   matrixRow.setElement( 2, index                                           ,     cCoef );

   matrixRow.setElement( 3, mesh.template getCellNextToCell< 1,0 >( index ),      dCoef );

   matrixRow.setElement( 4, mesh.template getCellNextToCell< 0,1 >( index ),      eCoef );

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index,

 	  typename OperatorQ >

tnlString

tnlFiniteVolumeNonlinearOperator< tnlGrid< 3, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::

getType()

{

   return tnlString( "tnlFiniteVolumeNonlinearOperator< " ) +

          MeshType::getType() + ", " +

          ::getType< Real >() + ", " +

          ::getType< Index >() + ", " +

	  OperatorQ::getType() + " >";

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index,

          typename OperatorQ >

template< typename Vector >

#ifdef HAVE_CUDA

__device__ __host__

#endif

Real

tnlFiniteVolumeNonlinearOperator< tnlGrid< 3, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::

getValue( const MeshType& mesh,

          const IndexType cellIndex,

          const CoordinatesType& coordinates,

          const Vector& u,

          const Real& time ) const

{

   return operatorQ.getValue(mesh,cellIndex, coordinates, u, time)*( (u[ mesh.template getCellNextToCell< 1,0,0 >( cellIndex ) ] - u[ cellIndex ]) 

          * mesh.getHxSquareInverse() / operatorQ.getValue(mesh,cellIndex, coordinates, u, time, 1 )

          + ( u[ mesh.template getCellNextToCell< 0,1,0 >( cellIndex ) ] - u[ cellIndex ]) * mesh.getHySquareInverse()/

          operatorQ.getValue(mesh,cellIndex, coordinates, u, time, 0, 1 ) 

          + ( u[ mesh.template getCellNextToCell< 0,0,1 >( cellIndex ) ] - u[ cellIndex ]) * mesh.getHzSquareInverse()/

          operatorQ.getValue(mesh,cellIndex, coordinates, u, time, 0, 0, 1 ) 

          - ( - u[ mesh.template getCellNextToCell< -1,0,0 >( cellIndex ) ]  + u[ cellIndex ]) 

          * mesh.getHxSquareInverse()/operatorQ.getValue(mesh,cellIndex, coordinates, u, time, -1)

          -( - u[ mesh.template getCellNextToCell< 0,-1,0 >( cellIndex ) ] + u[ cellIndex ]) * mesh.getHySquareInverse()

          /operatorQ.getValue(mesh,cellIndex,coordinates, u, time, 0, -1) 

          -( - u[ mesh.template getCellNextToCell< 0,0,-1 >( cellIndex ) ] + u[ cellIndex ]) * mesh.getHzSquareInverse()

          /operatorQ.getValue(mesh,cellIndex,coordinates, u, time, 0, 0, -1) );

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index,

          typename OperatorQ >

#ifdef HAVE_CUDA

   __device__ __host__

#endif

Index

tnlFiniteVolumeNonlinearOperator< tnlGrid< 3, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::

getLinearSystemRowLength( const MeshType& mesh,

                          const IndexType& index,

                          const CoordinatesType& coordinates ) const

{

   return 7;

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index,

          typename OperatorQ >

template< typename Vector, typename MatrixRow >

#ifdef HAVE_CUDA

__device__ __host__

#endif

void

tnlFiniteVolumeNonlinearOperator< tnlGrid< 3, MeshReal, Device, MeshIndex >, OperatorQ, Real, Index >::

updateLinearSystem( const RealType& time,

                    const RealType& tau,

                    const MeshType& mesh,

                    const IndexType& index,

                    const CoordinatesType& coordinates,

                    Vector& u,

                    Vector& b,

                    MatrixRow& matrixRow ) const

{

   const RealType aCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time ) * mesh.getHzSquareInverse() / 

                       operatorQ.getValue(mesh, index, coordinates, u, time, 0, 0, -1 );

   const RealType bCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time ) * mesh.getHySquareInverse() / 

                       operatorQ.getValue(mesh, index, coordinates, u, time, 0, -1, 0 );

   const RealType cCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time ) * mesh.getHxSquareInverse() / 

                       operatorQ.getValue(mesh, index, coordinates, u, time, -1, 0, 0 );

   const RealType dCoef = tau * operatorQ.getValue(mesh, index, coordinates, u, time ) * ( mesh.getHxSquareInverse() / 

                       operatorQ.getValue(mesh, index, coordinates, u, time, 1, 0, 0 ) + mesh.getHySquareInverse() / 

                       operatorQ.getValue(mesh, index, coordinates, u, time, 0, 1, 0 )

                       + mesh.getHzSquareInverse() / operatorQ.getValue(mesh, index, coordinates, u, time, 0, 0, 1 ) + 

                       mesh.getHxSquareInverse() / operatorQ.getValue(mesh, index, coordinates, u, time, -1, 0, 0 )

                       + mesh.getHySquareInverse() / operatorQ.getValue(mesh, index, coordinates, u, time, 0, -1, 0 ) + 

                       mesh.getHzSquareInverse() / operatorQ.getValue(mesh, index, coordinates, u, time, 0, 0, -1 ) );

   const RealType eCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time, 1, 0, 0 ) * mesh.getHxSquareInverse() / 

                       operatorQ.getValue(mesh, index, coordinates, u, time );

   const RealType fCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time, 0, 1, 0 ) * mesh.getHySquareInverse() / 

                       operatorQ.getValue(mesh, index, coordinates, u, time );

   const RealType gCoef = - tau * operatorQ.getValue(mesh, index, coordinates, u, time, 0, 0, 1 ) * mesh.getHzSquareInverse() / 

                       operatorQ.getValue(mesh, index, coordinates, u, time );

   matrixRow.setElement( 0, mesh.template getCellNextToCell< 0,0,-1 >( index ),     aCoef );

   matrixRow.setElement( 1, mesh.template getCellNextToCell< 0,-1,0 >( index ),     bCoef );

   matrixRow.setElement( 2, mesh.template getCellNextToCell< -1,0,0 >( index ),     cCoef );

   matrixRow.setElement( 3, index,                                                  dCoef );

   matrixRow.setElement( 4, mesh.template getCellNextToCell< 1,0,0 >( index ),      eCoef );

   matrixRow.setElement( 5, mesh.template getCellNextToCell< 0,1,0 >( index ),      fCoef );

   matrixRow.setElement( 6, mesh.template getCellNextToCell< 0,0,1 >( index ),      gCoef );

}

#endif	/* TNLFINITEVOLUMENONLINEAROPERATOR__IMPL_H */

src/operators/operator-Q/CMakeLists.txt

+3 −1
Original line number Diff line number Diff line 
SET( headers tnlExactOperatorQ.h

             tnlExactOperatorQ_impl.h

             tnlOneSideDiffOperatorQ.h             

             tnlOneSideDiffOperatorQ_impl.h )

             tnlOneSideDiffOperatorQ_impl.h

	     tnlFiniteVolumeOperatorQ.h             

             tnlFiniteVolumeOperatorQ_impl.h )



INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/operators/operator-Q )