Merge branch 'diffusion' of geraldine.fjfi.cvut.cz:/local/projects/tnl/tnl into diffusion (ef8c7f8c) · Commits · TNL / tnl-dev

examples/heat-equation/heatEquationSolver.h

+13 −3

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"


		@@ -54,6 +56,9 @@ class heatEquationSolver
		bool setInitialCondition( const tnlParameterContainer& parameters,
		const MeshType& mesh );

		bool setupLinearSystem( const MeshType& mesh,
		MatrixType& matrix );

		bool makeSnapshot( const RealType& time,
		const IndexType& step,
		const MeshType& mesh );
		@@ -68,20 +73,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

+108 −41

Original line number	Diff line number	Diff line
		@@ -19,34 +19,51 @@
		#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 >
		::setup( const tnlParameterContainer& parameters )
		template< typename Mesh,
		typename DifferentialOperator,
		typename BoundaryCondition,
		typename RightHandSide >
		bool
		heatEquationSolver< Mesh, DifferentialOperator, BoundaryCondition, RightHandSide >::
		setup( const tnlParameterContainer& parameters )
		{
		if( ! boundaryCondition.setup( parameters ) \|\|
		! rightHandSide.setup( parameters ) )
		@@ -55,11 +72,12 @@ bool heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide >
		}

		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,11 +86,12 @@ 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
		@@ -80,11 +99,11 @@ typename heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide >::I
		}

		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,19 +112,19 @@ 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 )
		{
		@@ -119,13 +138,29 @@ bool heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide >
		}

		template< typename Mesh,
		typename Diffusion,
		typename DifferentialOperator,
		typename BoundaryCondition,
		typename RightHandSide >
		bool heatEquationSolver< Mesh,
		Diffusion,
		BoundaryCondition,
		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, DifferentialOperator, BoundaryCondition, RightHandSide >::
		makeSnapshot( const RealType& time,
		const IndexType& step,
		const MeshType& mesh )
		@@ -139,9 +174,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 +196,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,
		@@ -166,11 +206,38 @@ void heatEquationSolver< Mesh,Diffusion,BoundaryCondition,RightHandSide >
		_fu );
		}

		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 > systemAssembler;
		systemAssembler.template assembly< Mesh::Dimensions >( time,
		mesh,
		differentialOperator,
		boundaryConditions,
		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;
		}

src/implementation/matrices/CMakeLists.txt

+2 −1

Original line number	Diff line number	Diff line
		@@ -8,7 +8,8 @@ SET( headers tnlMatrix_impl.h
		tnlChunkedEllpackMatrix_impl.h
		tnlCSRMatrix_impl.h
		tnlMatrixReader_impl.h
		tnlMatrixWriter_impl.h )
		tnlMatrixWriter_impl.h
		tnlMatrixSetter_impl.h )

		SET( CURRENT_DIR ${CMAKE_SOURCE_DIR}/src/implementation/matrices )
		set( common_SOURCES

src/implementation/matrices/tnlMatrixSetter_impl.h

0 → 100644

+75 −0

Original line number	Diff line number	Diff line
		/***************************************************************************
		tnlMatrixSetter_impl.h - description
		-------------------
		begin : Oct 11, 2014
		copyright : (C) 2014 by Tomas Oberhuber
		email : tomas.oberhuber@fjfi.cvut.cz
		***************************************************************************/

		/***************************************************************************
		* *
		* This program is free software; you can redistribute it and/or modify *
		* it under the terms of the GNU General Public License as published by *
		* the Free Software Foundation; either version 2 of the License, or *
		* (at your option) any later version. *
		* *
		***************************************************************************/

		#ifndef TNLMATRIXSETTER_IMPL_H_
		#define TNLMATRIXSETTER_IMPL_H_

		template< typename Mesh,
		typename DifferentialOperator,
		typename BoundaryConditions,
		typename RowLengthsVector >
		template< int EntityDimensions >
		void
		tnlMatrixSetter< Mesh, DifferentialOperator, BoundaryConditions, RowLengthsVector >::
		getRowLengths( const Mesh& mesh,
		DifferentialOperator& differentialOperator,
		BoundaryConditions& boundaryConditions,
		RowLengthsVector& rowLengths ) const
		{
		TraversalUserData userData( differentialOperator, boundaryConditions, rowLengths );
		TraversalBoundaryEntitiesProcessor boundaryEntitiesProcessor;
		TraversalInteriorEntitiesProcessor interiorEntitiesProcessor;
		tnlTraversal< MeshType, EntityDimensions > meshTraversal;
		meshTraversal.template processEntities< TraversalUserData,
		TraversalBoundaryEntitiesProcessor,
		TraversalInteriorEntitiesProcessor >
		( mesh,
		userData,
		boundaryEntitiesProcessor,
		interiorEntitiesProcessor );
		}

		template< int Dimensions,
		typename Real,
		typename Device,
		typename Index,
		typename DifferentialOperator,
		typename BoundaryConditions,
		typename RowLengthsVector >
		template< int EntityDimensions >
		void
		tnlMatrixSetter< tnlGrid< Dimensions, Real, Device, Index >, DifferentialOperator, BoundaryConditions, RowLengthsVector >::
		getRowLengths( const MeshType& mesh,
		DifferentialOperator& differentialOperator,
		BoundaryConditions& boundaryConditions,
		RowLengthsVector& rowLengths ) const
		{
		TraversalUserData userData( differentialOperator, boundaryConditions, rowLengths );
		TraversalBoundaryEntitiesProcessor boundaryEntitiesProcessor;
		TraversalInteriorEntitiesProcessor interiorEntitiesProcessor;
		tnlTraversal< MeshType, EntityDimensions > meshTraversal;
		meshTraversal.template processEntities< TraversalUserData,
		TraversalBoundaryEntitiesProcessor,
		TraversalInteriorEntitiesProcessor >
		( mesh,
		userData,
		boundaryEntitiesProcessor,
		interiorEntitiesProcessor );
		}


		#endif /* TNLMATRIXSETTER_IMPL_H_ */

src/implementation/operators/diffusion/tnlLinearDiffusion_impl.h

+43 −72

Original line number	Diff line number	Diff line
		@@ -20,7 +20,7 @@ getType()
		::getType< Index >() + " >";
		}

		/*template< typename MeshReal,
		template< typename MeshReal,
		typename Device,
		typename MeshIndex,
		typename Real,
		@@ -29,41 +29,34 @@ getType()
		#ifdef HAVE_CUDA
		__device__ __host__
		#endif
		void
		Real
		tnlLinearDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
		explicitUpdate( const RealType& time,
		const RealType& tau,
		const MeshType& mesh,
		getValue( const MeshType& mesh,
		const IndexType cellIndex,
		const CoordinatesType& coordinates,
		Vector& u,
		Vector& fu ) const
		const Vector& u,
		const Real& time ) const
		{
		fu[ cellIndex ] = ( u[ mesh.getCellXPredecessor( cellIndex ) ]
		return ( u[ mesh.getCellXPredecessor( cellIndex ) ]
		- 2.0 * u[ cellIndex ]
		+ u[ mesh.getCellXSuccessor( cellIndex ) ] ) * mesh.getHxSquareInverse();
		}*/
		}

		template< typename MeshReal,
		typename Device,
		typename MeshIndex,
		typename Real,
		typename Index >
		template< typename Vector >
		#ifdef HAVE_CUDA
		__device__ __host__
		#endif
		Real
		Index
		tnlLinearDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
		getValue( const MeshType& mesh,
		const IndexType cellIndex,
		const CoordinatesType& coordinates,
		const Vector& u,
		const Real& time ) const
		getLinearSystemRowLength( const MeshType& mesh,
		const IndexType& index,
		const CoordinatesType& coordinates ) const
		{
		return ( u[ mesh.getCellXPredecessor( cellIndex ) ]
		- 2.0 * u[ cellIndex ]
		+ u[ mesh.getCellXSuccessor( cellIndex ) ] ) * mesh.getHxSquareInverse();
		return 3;
		}

		template< typename MeshReal,
		@@ -81,32 +74,23 @@ getType()
		::getType< Index >() + " >";
		}

		/*template< typename MeshReal,
		template< typename MeshReal,
		typename Device,
		typename MeshIndex,
		typename Real,
		typename Index >
		template< typename Vector >
		#ifdef HAVE_CUDA
		__device__ __host__
		#endif
		void
		Index
		tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
		explicitUpdate( const RealType& time,
		const RealType& tau,
		const MeshType& mesh,
		const IndexType cellIndex,
		const CoordinatesType& coordinates,
		Vector& u,
		Vector& fu ) const
		getLinearSystemRowLength( const MeshType& mesh,
		const IndexType& index,
		const CoordinatesType& coordinates ) const
		{
		fu[ cellIndex ] = ( u[ mesh.getCellXPredecessor( cellIndex ) ]
		- 2.0 * u[ cellIndex ]
		+ u[ mesh.getCellXSuccessor( cellIndex ) ] ) * mesh.getHxSquareInverse() +
		( u[ mesh.getCellYPredecessor( cellIndex ) ]
		- 2.0 * u[ cellIndex ]
		+ u[ mesh.getCellYSuccessor( cellIndex ) ] ) * mesh.getHySquareInverse();
		}*/
		return 5;
		}


		template< typename MeshReal,
		typename Device,
		@@ -148,7 +132,7 @@ getType()
		::getType< Index >() + " >";
		}

		/*template< typename MeshReal,
		template< typename MeshReal,
		typename Device,
		typename MeshIndex,
		typename Real,
		@@ -157,17 +141,15 @@ getType()
		#ifdef HAVE_CUDA
		__device__ __host__
		#endif
		void
		Real
		tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
		explicitUpdate( const RealType& time,
		const RealType& tau,
		const MeshType& mesh,
		getValue( const MeshType& mesh,
		const IndexType cellIndex,
		const CoordinatesType& coordinates,
		Vector& u,
		Vector& fu ) const
		const Vector& u,
		const Real& time ) const
		{
		fu[ cellIndex ] = ( u[ mesh.getCellXPredecessor( cellIndex ) ]
		return ( u[ mesh.getCellXPredecessor( cellIndex ) ]
		- 2.0 * u[ cellIndex ]
		+ u[ mesh.getCellXSuccessor( cellIndex ) ] ) * mesh.getHxSquareInverse() +
		( u[ mesh.getCellYPredecessor( cellIndex ) ]
		@@ -176,34 +158,23 @@ explicitUpdate( const RealType& time,
		( u[ mesh.getCellZPredecessor( cellIndex ) ]
		- 2.0 * u[ cellIndex ]
		+ u[ mesh.getCellZSuccessor( cellIndex ) ] ) * mesh.getHzSquareInverse();
		}*/
		}

		template< typename MeshReal,
		typename Device,
		typename MeshIndex,
		typename Real,
		typename Index >
		template< typename Vector >
		#ifdef HAVE_CUDA
		__device__ __host__
		#endif
		Real
		Index
		tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
		getValue( const MeshType& mesh,
		const IndexType cellIndex,
		const CoordinatesType& coordinates,
		const Vector& u,
		const Real& time ) const
		getLinearSystemRowLength( const MeshType& mesh,
		const IndexType& index,
		const CoordinatesType& coordinates ) const
		{
		return ( u[ mesh.getCellXPredecessor( cellIndex ) ]
		- 2.0 * u[ cellIndex ]
		+ u[ mesh.getCellXSuccessor( cellIndex ) ] ) * mesh.getHxSquareInverse() +
		( u[ mesh.getCellYPredecessor( cellIndex ) ]
		- 2.0 * u[ cellIndex ]
		+ u[ mesh.getCellYSuccessor( cellIndex ) ] ) * mesh.getHySquareInverse() +
		( u[ mesh.getCellZPredecessor( cellIndex ) ]
		- 2.0 * u[ cellIndex ]
		+ u[ mesh.getCellZSuccessor( cellIndex ) ] ) * mesh.getHzSquareInverse();
		return 7;
		}