From 13bde21f3de637614468f2c1e15286561c2e3602 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= <oberhuber.tomas@gmail.com> Date: Sat, 18 Feb 2017 12:42:06 +0100 Subject: [PATCH] Refactoring the inviscid flow solver. --- examples/inviscid-flow/1d/CMakeLists.txt | 20 - .../inviscid-flow/1d/EulerPressureGetter.h | 63 --- examples/inviscid-flow/1d/EulerVelGetter.h | 56 --- examples/inviscid-flow/1d/LaxFridrichs.h | 36 -- .../inviscid-flow/1d/LaxFridrichsContinuity.h | 184 --------- .../1d/LaxFridrichsContinuity_impl.h | 345 ---------------- .../inviscid-flow/1d/LaxFridrichsEnergy.h | 200 --------- .../1d/LaxFridrichsEnergy_impl.h | 352 ---------------- .../inviscid-flow/1d/LaxFridrichsMomentum.h | 201 --------- .../1d/LaxFridrichsMomentum_impl.h | 347 ---------------- examples/inviscid-flow/1d/euler.cpp | 1 - examples/inviscid-flow/1d/euler.cu | 1 - examples/inviscid-flow/1d/euler.h | 119 ------ .../inviscid-flow/1d/eulerBuildConfigTag.h | 51 --- examples/inviscid-flow/1d/eulerProblem.h | 120 ------ examples/inviscid-flow/1d/eulerProblem_impl.h | 383 ------------------ examples/inviscid-flow/1d/eulerRhs.h | 35 -- examples/inviscid-flow/1d/tnl-run-euler-1d | 25 -- examples/inviscid-flow/CMakeLists.txt | 7 +- examples/inviscid-flow/EulerPressureGetter.h | 222 ---------- .../inviscid-flow/EulerPressureGetter_impl.h | 338 ---------------- examples/inviscid-flow/EulerVelGetter.h | 186 --------- examples/inviscid-flow/EulerVelGetter_impl.h | 340 ---------------- examples/inviscid-flow/EulerVelXGetter.h | 183 --------- examples/inviscid-flow/EulerVelXGetter_impl.h | 339 ---------------- examples/inviscid-flow/EulerVelYGetter.h | 183 --------- examples/inviscid-flow/EulerVelYGetter_impl.h | 332 --------------- examples/inviscid-flow/LaxFridrichs.h | 26 +- .../inviscid-flow/LaxFridrichsContinuity.h | 265 +++++++----- .../LaxFridrichsContinuity_impl .h | 349 ---------------- examples/inviscid-flow/LaxFridrichsEnergy.h | 317 +++++++++------ .../inviscid-flow/LaxFridrichsEnergy_impl.h | 350 ---------------- .../inviscid-flow/LaxFridrichsMomentumBase.h | 59 +++ .../inviscid-flow/LaxFridrichsMomentumX.h | 281 +++++++------ .../LaxFridrichsMomentumX_impl.h | 351 ---------------- .../inviscid-flow/LaxFridrichsMomentumY.h | 263 +++++++----- .../LaxFridrichsMomentumY_impl.h | 351 ---------------- .../inviscid-flow/LaxFridrichsMomentumZ.h | 239 +++++++++++ examples/inviscid-flow/eulerBuildConfigTag.h | 4 +- examples/inviscid-flow/eulerProblem.h | 5 +- examples/inviscid-flow/eulerProblem_impl.h | 67 +-- examples/inviscid-flow/run-euler | 19 +- src/TNL/Functions/VectorField.h | 2 + 43 files changed, 1060 insertions(+), 6557 deletions(-) delete mode 100644 examples/inviscid-flow/1d/CMakeLists.txt delete mode 100644 examples/inviscid-flow/1d/EulerPressureGetter.h delete mode 100644 examples/inviscid-flow/1d/EulerVelGetter.h delete mode 100644 examples/inviscid-flow/1d/LaxFridrichs.h delete mode 100644 examples/inviscid-flow/1d/LaxFridrichsContinuity.h delete mode 100644 examples/inviscid-flow/1d/LaxFridrichsContinuity_impl.h delete mode 100644 examples/inviscid-flow/1d/LaxFridrichsEnergy.h delete mode 100644 examples/inviscid-flow/1d/LaxFridrichsEnergy_impl.h delete mode 100644 examples/inviscid-flow/1d/LaxFridrichsMomentum.h delete mode 100644 examples/inviscid-flow/1d/LaxFridrichsMomentum_impl.h delete mode 100644 examples/inviscid-flow/1d/euler.cpp delete mode 100644 examples/inviscid-flow/1d/euler.cu delete mode 100644 examples/inviscid-flow/1d/euler.h delete mode 100644 examples/inviscid-flow/1d/eulerBuildConfigTag.h delete mode 100644 examples/inviscid-flow/1d/eulerProblem.h delete mode 100644 examples/inviscid-flow/1d/eulerProblem_impl.h delete mode 100644 examples/inviscid-flow/1d/eulerRhs.h delete mode 100644 examples/inviscid-flow/1d/tnl-run-euler-1d delete mode 100644 examples/inviscid-flow/EulerPressureGetter.h delete mode 100644 examples/inviscid-flow/EulerPressureGetter_impl.h delete mode 100644 examples/inviscid-flow/EulerVelGetter.h delete mode 100644 examples/inviscid-flow/EulerVelGetter_impl.h delete mode 100644 examples/inviscid-flow/EulerVelXGetter.h delete mode 100644 examples/inviscid-flow/EulerVelXGetter_impl.h delete mode 100644 examples/inviscid-flow/EulerVelYGetter.h delete mode 100644 examples/inviscid-flow/EulerVelYGetter_impl.h delete mode 100644 examples/inviscid-flow/LaxFridrichsContinuity_impl .h delete mode 100644 examples/inviscid-flow/LaxFridrichsEnergy_impl.h create mode 100644 examples/inviscid-flow/LaxFridrichsMomentumBase.h delete mode 100644 examples/inviscid-flow/LaxFridrichsMomentumX_impl.h delete mode 100644 examples/inviscid-flow/LaxFridrichsMomentumY_impl.h create mode 100644 examples/inviscid-flow/LaxFridrichsMomentumZ.h diff --git a/examples/inviscid-flow/1d/CMakeLists.txt b/examples/inviscid-flow/1d/CMakeLists.txt deleted file mode 100644 index 6fe0766eec..0000000000 --- a/examples/inviscid-flow/1d/CMakeLists.txt +++ /dev/null @@ -1,20 +0,0 @@ -set( tnl_heat_equation_SOURCES - euler.cpp - euler.cu ) - -IF( BUILD_CUDA ) - CUDA_ADD_EXECUTABLE(tnl-euler-1d${debugExt} euler.cu) - target_link_libraries (tnl-euler-1d${debugExt} tnl${debugExt}-${tnlVersion} ${CUSPARSE_LIBRARY} ) -ELSE( BUILD_CUDA ) - ADD_EXECUTABLE(tnl-euler-1d${debugExt} euler.cpp) - target_link_libraries (tnl-euler-1d${debugExt} tnl${debugExt}-${tnlVersion} ) -ENDIF( BUILD_CUDA ) - - -INSTALL( TARGETS tnl-euler-1d${debugExt} - RUNTIME DESTINATION bin - PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE ) - -INSTALL( FILES tnl-run-euler-1d - ${tnl_heat_equation_SOURCES} - DESTINATION share/tnl-${tnlVersion}/examples/inviscid-flow-1d ) diff --git a/examples/inviscid-flow/1d/EulerPressureGetter.h b/examples/inviscid-flow/1d/EulerPressureGetter.h deleted file mode 100644 index fc98847a27..0000000000 --- a/examples/inviscid-flow/1d/EulerPressureGetter.h +++ /dev/null @@ -1,63 +0,0 @@ -#ifndef EulerPressureGetter_H -#define EulerPressureGetter_H - -#include <TNL/Containers/Vector.h> -#include <TNL/Meshes/Grid.h> -#include <TNL/Functions/Domain.h> - -namespace TNL { - -template< typename Mesh, - typename Real = typename Mesh::RealType, - typename Index = typename Mesh::IndexType > -class EulerPressureGetter -: public Functions::Domain< Mesh::getMeshDimensions(), Functions::MeshDomain > -{ - public: - - typedef Mesh MeshType; - typedef typename MeshType::DeviceType DeviceType; - typedef Real RealType; - typedef Index IndexType; - typedef Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - - EulerPressureGetter( const MeshFunctionType& velocity, - const MeshFunctionType& rhoVel, - const MeshFunctionType& energy, - const RealType& gamma ) - : rho( rho ), rhoVel( rhoVel ), energy( energy ), gamma( gamma ) - {} - - template< typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshEntity& entity, - const RealType& time = 0.0 ) const - { - return this->operator[]( entity.getIndex() ); - } - - __cuda_callable__ - Real operator[]( const IndexType& idx ) const - { - return ( this->gamma - 1.0 ) * ( this->energy[ idx ] - 0.5 * this->rhoVel[ idx ] * this->rhoVel[ idx ] / this->rho[ idx ]); - } - - - protected: - - Real gamma; - - const MeshFunctionType& rho; - - const MeshFunctionType& rhoVel; - - const MeshFunctionType& energy; - -}; - -} //namespace TNL - -#endif /* EulerPressureGetter_H */ diff --git a/examples/inviscid-flow/1d/EulerVelGetter.h b/examples/inviscid-flow/1d/EulerVelGetter.h deleted file mode 100644 index 7e58ae77d6..0000000000 --- a/examples/inviscid-flow/1d/EulerVelGetter.h +++ /dev/null @@ -1,56 +0,0 @@ -#ifndef EulerVelGetter_H -#define EulerVelGetter_H - -#include <TNL/Containers/Vector.h> -#include <TNL/Meshes/Grid.h> - -namespace TNL { - -template< typename Mesh, - typename Real = typename Mesh::RealType, - typename Index = typename Mesh::IndexType > -class EulerVelGetter -: public Functions::Domain< Mesh::getMeshDimensions(), Functions::MeshDomain > -{ - public: - - typedef Mesh MeshType; - typedef typename MeshType::DeviceType DeviceType; - typedef Real RealType; - typedef Index IndexType; - typedef Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - - EulerVelGetter( const MeshFunctionType& rho, - const MeshFunctionType& rhoVel) - : rho( rho ), rhoVel( rhoVel ) - {} - - template< typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshEntity& entity, - const RealType& time = 0.0 ) const - { - return this->operator[]( entity.getIndex() ); - } - - __cuda_callable__ - Real operator[]( const IndexType& idx ) const - { - return this->rho[ idx ] / this->rhoVel[ idx ]; - } - - - protected: - - const MeshFunctionType& rho; - - const MeshFunctionType& rhoVel; - -}; - -} // namespace TNL - -#endif /* EulerVelGetter_H */ diff --git a/examples/inviscid-flow/1d/LaxFridrichs.h b/examples/inviscid-flow/1d/LaxFridrichs.h deleted file mode 100644 index a5cb7eea5c..0000000000 --- a/examples/inviscid-flow/1d/LaxFridrichs.h +++ /dev/null @@ -1,36 +0,0 @@ -#ifndef LaxFridrichs_H -#define LaxFridrichs_H - -#include <TNL/Containers/Vector.h> -#include <TNL/Meshes/Grid.h> - -#include "LaxFridrichsContinuity.h" -#include "LaxFridrichsMomentum.h" -#include "LaxFridrichsEnergy.h" -#include "EulerVelGetter.h" -#include "EulerPressureGetter.h" - -namespace TNL { - -template< typename Mesh, - typename Real = typename Mesh::RealType, - typename Index = typename Mesh::IndexType > -class LaxFridrichs -{ - public: - typedef Real RealType; - typedef typename Mesh::DeviceType DeviceType; - typedef Index IndexType; - typedef Functions::MeshFunction< Mesh > MeshFunctionType; - - typedef LaxFridrichsContinuity< Mesh, Real, Index > Continuity; - typedef LaxFridrichsMomentum< Mesh, Real, Index > Momentum; - typedef LaxFridrichsEnergy< Mesh, Real, Index > Energy; - typedef EulerVelGetter< Mesh, Real, Index > Velocity; - typedef EulerPressureGetter< Mesh, Real, Index > Pressure; - -}; - -} // namespace TNL - -#endif /* LaxFridrichs_H */ diff --git a/examples/inviscid-flow/1d/LaxFridrichsContinuity.h b/examples/inviscid-flow/1d/LaxFridrichsContinuity.h deleted file mode 100644 index 8053895eda..0000000000 --- a/examples/inviscid-flow/1d/LaxFridrichsContinuity.h +++ /dev/null @@ -1,184 +0,0 @@ -#ifndef LaxFridrichsContinuity_H -#define LaxFridrichsContinuity_H - -#include <TNL/Containers/Vector.h> -#include <TNL/Meshes/Grid.h> - -namespace TNL { - -template< typename Mesh, - typename Real = typename Mesh::RealType, - typename Index = typename Mesh::IndexType > -class LaxFridrichsContinuity -{ -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class LaxFridrichsContinuity< Meshes::Grid< 1,MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocity; - - void setTau(const Real& tau) - { - this->tau = tau; - }; - - void setVelocity(const MeshFunctionType velocity) - { - //this->velocity = velocity; - }; - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class LaxFridrichsContinuity< Meshes::Grid< 2,MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocity; - - void setTau(const Real& tau) - { - this->tau = tau; - }; - - void setVelocity(const MeshFunctionType& velocity) - { - this->velocity = velocity; - }; - - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class LaxFridrichsContinuity< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocity; - - void setTau(const Real& tau) - { - this->tau = tau; - }; - - void setVelocity(const MeshFunctionType& velocity) - { - this->velocity = velocity; - }; - - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -} // namespace TNL - -#include "LaxFridrichsContinuity_impl.h" - -#endif /* LaxFridrichsContinuity_H */ diff --git a/examples/inviscid-flow/1d/LaxFridrichsContinuity_impl.h b/examples/inviscid-flow/1d/LaxFridrichsContinuity_impl.h deleted file mode 100644 index 3be1a5b6e4..0000000000 --- a/examples/inviscid-flow/1d/LaxFridrichsContinuity_impl.h +++ /dev/null @@ -1,345 +0,0 @@ -#ifndef LaxFridrichsContinuity_IMPL_H -#define LaxFridrichsContinuity_IMPL_H - -namespace TNL { - -/**** - * 1D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsContinuity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsContinuity< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsContinuity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - static_assert( MeshEntity::entityDimensions == 1, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 1, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - - //rho - const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - return (0.5 * this->tau) * ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) - - 0.5 * hxInverse * ( u[ west ] * this->velocity[ west ] - u[ east ] * this->velocity[ east ] ); -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsContinuity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsContinuity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - matrixRow.setElement( 0, west, - lambdaX ); - matrixRow.setElement( 1, center, 2.0 * lambdaX ); - matrixRow.setElement( 2, east, - lambdaX ); -} - -/**** - * 2D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsContinuity< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsContinuity< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsContinuity< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 2, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 2, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - - const RealType& hxSquareInverse = entity.getMesh().template getSpaceStepsProducts< -2, 0 >(); - const RealType& hySquareInverse = entity.getMesh().template getSpaceStepsProducts< 0, -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse + - ( u[ south ] - 2.0 * u[ center ] + u[ north ] ) * hySquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsContinuity< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsContinuity< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2, 0 >(); - const RealType& lambdaY = tau * entity.getMesh().template getSpaceStepsProducts< 0, -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); - matrixRow.setElement( 0, south, -lambdaY ); - matrixRow.setElement( 1, west, -lambdaX ); - matrixRow.setElement( 2, center, 2.0 * ( lambdaX + lambdaY ) ); - matrixRow.setElement( 3, east, -lambdaX ); - matrixRow.setElement( 4, north, -lambdaY ); -} - -/**** - * 3D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsContinuity< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsContinuity< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsContinuity< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 3, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 3, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse + - ( u[ south ] - 2.0 * u[ center ] + u[ north ] ) * hySquareInverse + - ( u[ up ] - 2.0 * u[ center ] + u[ down ] ) * hzSquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsContinuity< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsContinuity< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - matrixRow.setElement( 0, down, -lambdaZ ); - matrixRow.setElement( 1, south, -lambdaY ); - matrixRow.setElement( 2, west, -lambdaX ); - matrixRow.setElement( 3, center, 2.0 * ( lambdaX + lambdaY + lambdaZ ) ); - matrixRow.setElement( 4, east, -lambdaX ); - matrixRow.setElement( 5, north, -lambdaY ); - matrixRow.setElement( 6, up, -lambdaZ ); -} - -} // namespace TNL - -#endif /* LaxFridrichsContinuityIMPL_H */ - diff --git a/examples/inviscid-flow/1d/LaxFridrichsEnergy.h b/examples/inviscid-flow/1d/LaxFridrichsEnergy.h deleted file mode 100644 index fa0012eb8d..0000000000 --- a/examples/inviscid-flow/1d/LaxFridrichsEnergy.h +++ /dev/null @@ -1,200 +0,0 @@ -#ifndef LaxFridrichsEnergy_H -#define LaxFridrichsEnergy_H - -#include <TNL/Containers/Vector.h> -#include <TNL/Meshes/Grid.h> - -namespace TNL { - -template< typename Mesh, - typename Real = typename Mesh::RealType, - typename Index = typename Mesh::IndexType > -class LaxFridrichsEnergy -{ -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class LaxFridrichsEnergy< Meshes::Grid< 1,MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocity; - MeshFunctionType pressure; - - void setTau(const Real& tau) - { - this->tau = tau; - }; - - void setVelocity(const MeshFunctionType& velocity) - { - this->velocity = velocity; - }; - - void setPressure(const MeshFunctionType& pressure) - { - this->pressure = pressure; - }; - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class LaxFridrichsEnergy< Meshes::Grid< 2,MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocity; - MeshFunctionType pressure; - - void setTau(const Real& tau) - { - this->tau = tau; - }; - - void setVelocity(const MeshFunctionType& velocity) - { - this->velocity = velocity; - }; - - void setPressure(const MeshFunctionType& pressure) - { - this->pressure = pressure; - }; - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class LaxFridrichsEnergy< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocity; - MeshFunctionType pressure; - - void setTau(const Real& tau) - { - this->tau = tau; - }; - - void setVelocity(const MeshFunctionType& velocity) - { - this->velocity = velocity; - }; - - void setPressure(const MeshFunctionType& pressure) - { - this->pressure = pressure; - }; - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -} // namespace TNL - -#include "LaxFridrichsEnergy_impl.h" - -#endif /* LaxFridrichsEnergy_H */ diff --git a/examples/inviscid-flow/1d/LaxFridrichsEnergy_impl.h b/examples/inviscid-flow/1d/LaxFridrichsEnergy_impl.h deleted file mode 100644 index ba8092e807..0000000000 --- a/examples/inviscid-flow/1d/LaxFridrichsEnergy_impl.h +++ /dev/null @@ -1,352 +0,0 @@ -#ifndef LaxFridrichsEnergy_IMPL_H -#define LaxFridrichsEnergy_IMPL_H - -namespace TNL { - -/**** - * 1D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsEnergy< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsEnergy< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsEnergy< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 1, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 1, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - - //energy - const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - return (0.5 * this->tau) * ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) - - 0.5 * hxInverse * - (( u[ west ] + this->pressure[ west ] ) * velocity[ west ] - - (u[ east ] + this->pressure[ east ] ) * velocity[ east ] ); -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsEnergy< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsEnergy< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - matrixRow.setElement( 0, west, - lambdaX ); - matrixRow.setElement( 1, center, 2.0 * lambdaX ); - matrixRow.setElement( 2, east, - lambdaX ); -} - -/**** - * 2D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsEnergy< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsEnergy< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsEnergy< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 2, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 2, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - - const RealType& hxSquareInverse = entity.getMesh().template getSpaceStepsProducts< -2, 0 >(); - const RealType& hySquareInverse = entity.getMesh().template getSpaceStepsProducts< 0, -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse + - ( u[ south ] - 2.0 * u[ center ] + u[ north ] ) * hySquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsEnergy< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsEnergy< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2, 0 >(); - const RealType& lambdaY = tau * entity.getMesh().template getSpaceStepsProducts< 0, -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); - matrixRow.setElement( 0, south, -lambdaY ); - matrixRow.setElement( 1, west, -lambdaX ); - matrixRow.setElement( 2, center, 2.0 * ( lambdaX + lambdaY ) ); - matrixRow.setElement( 3, east, -lambdaX ); - matrixRow.setElement( 4, north, -lambdaY ); -} - -/**** - * 3D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsEnergy< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsEnergy< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsEnergy< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 3, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 3, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse + - ( u[ south ] - 2.0 * u[ center ] + u[ north ] ) * hySquareInverse + - ( u[ up ] - 2.0 * u[ center ] + u[ down ] ) * hzSquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsEnergy< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsEnergy< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - matrixRow.setElement( 0, down, -lambdaZ ); - matrixRow.setElement( 1, south, -lambdaY ); - matrixRow.setElement( 2, west, -lambdaX ); - matrixRow.setElement( 3, center, 2.0 * ( lambdaX + lambdaY + lambdaZ ) ); - matrixRow.setElement( 4, east, -lambdaX ); - matrixRow.setElement( 5, north, -lambdaY ); - matrixRow.setElement( 6, up, -lambdaZ ); -} - -} //namespace TNL - -#endif /* LaxFridrichsEnergyIMPL_H */ - diff --git a/examples/inviscid-flow/1d/LaxFridrichsMomentum.h b/examples/inviscid-flow/1d/LaxFridrichsMomentum.h deleted file mode 100644 index 5bb8818833..0000000000 --- a/examples/inviscid-flow/1d/LaxFridrichsMomentum.h +++ /dev/null @@ -1,201 +0,0 @@ -#ifndef LaxFridrichsMomentum_H -#define LaxFridrichsMomentum_H - -#include <TNL/Containers/Vector.h> -#include <TNL/Meshes/Grid.h> - -namespace TNL { - -template< typename Mesh, - typename Real = typename Mesh::RealType, - typename Index = typename Mesh::IndexType > -class LaxFridrichsMomentum -{ -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class LaxFridrichsMomentum< Meshes::Grid< 1,MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocity; - MeshFunctionType pressure; - - void setTau(const Real& tau) - { - this->tau = tau; - }; - - void setVelocity(const MeshFunctionType& velocity) - { - this->velocity = velocity; - }; - - void setPressure(const MeshFunctionType& pressure) - { - this->pressure = pressure; - }; - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class LaxFridrichsMomentum< Meshes::Grid< 2,MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocity; - MeshFunctionType pressure; - - void setTau(const Real& tau) - { - this->tau = tau; - }; - - void setVelocity(const MeshFunctionType& velocity) - { - this->velocity = velocity; - }; - - void setPressure(const MeshFunctionType& pressure) - { - this->pressure = pressure; - }; - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class LaxFridrichsMomentum< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocity; - MeshFunctionType pressure; - - void setTau(const Real& tau) - { - this->tau = tau; - }; - - void setVelocity(const MeshFunctionType& velocity) - { - this->velocity = velocity; - }; - - void setPressure(const MeshFunctionType& pressure) - { - this->pressure = pressure; - }; - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -} // namespace TNL - - -#include "LaxFridrichsMomentum_impl.h" - -#endif /* LaxFridrichsMomentum_H */ diff --git a/examples/inviscid-flow/1d/LaxFridrichsMomentum_impl.h b/examples/inviscid-flow/1d/LaxFridrichsMomentum_impl.h deleted file mode 100644 index 75c7a7fea0..0000000000 --- a/examples/inviscid-flow/1d/LaxFridrichsMomentum_impl.h +++ /dev/null @@ -1,347 +0,0 @@ -#ifndef LaxFridrichsMomentum_IMPL_H -#define LaxFridrichsMomentum_IMPL_H - -namespace TNL { - -/**** - * 1D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsMomentum< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsMomentum< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsMomentum< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - static_assert( MeshEntity::entityDimensions == 1, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 1, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - - //rhoVelocity - const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - return (0.5 * this->tau) * ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) - - 0.5 * hxInverse * - (( u[ west ] * this -> velocity[ west ] + this -> pressure [ west ] ) - - (u[ east ] * this -> velocity[ east ] + this -> pressure [ east ] )); -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsMomentum< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsMomentum< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - matrixRow.setElement( 0, west, - lambdaX ); - matrixRow.setElement( 1, center, 2.0 * lambdaX ); - matrixRow.setElement( 2, east, - lambdaX ); -} - -/**** - * 2D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsMomentum< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsMomentum< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsMomentum< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 2, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 2, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - - const RealType& hxSquareInverse = entity.getMesh().template getSpaceStepsProducts< -2, 0 >(); - const RealType& hySquareInverse = entity.getMesh().template getSpaceStepsProducts< 0, -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse + - ( u[ south ] - 2.0 * u[ center ] + u[ north ] ) * hySquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsMomentum< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsMomentum< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2, 0 >(); - const RealType& lambdaY = tau * entity.getMesh().template getSpaceStepsProducts< 0, -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); - matrixRow.setElement( 0, south, -lambdaY ); - matrixRow.setElement( 1, west, -lambdaX ); - matrixRow.setElement( 2, center, 2.0 * ( lambdaX + lambdaY ) ); - matrixRow.setElement( 3, east, -lambdaX ); - matrixRow.setElement( 4, north, -lambdaY ); -} - -/**** - * 3D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsMomentum< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsMomentum< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsMomentum< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 3, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 3, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse + - ( u[ south ] - 2.0 * u[ center ] + u[ north ] ) * hySquareInverse + - ( u[ up ] - 2.0 * u[ center ] + u[ down ] ) * hzSquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsMomentum< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsMomentum< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - matrixRow.setElement( 0, down, -lambdaZ ); - matrixRow.setElement( 1, south, -lambdaY ); - matrixRow.setElement( 2, west, -lambdaX ); - matrixRow.setElement( 3, center, 2.0 * ( lambdaX + lambdaY + lambdaZ ) ); - matrixRow.setElement( 4, east, -lambdaX ); - matrixRow.setElement( 5, north, -lambdaY ); - matrixRow.setElement( 6, up, -lambdaZ ); -} - -} // namespace TNL - -#endif /* LaxFridrichsMomentumIMPL_H */ - diff --git a/examples/inviscid-flow/1d/euler.cpp b/examples/inviscid-flow/1d/euler.cpp deleted file mode 100644 index 4d76005cb1..0000000000 --- a/examples/inviscid-flow/1d/euler.cpp +++ /dev/null @@ -1 +0,0 @@ -#include "euler.h" diff --git a/examples/inviscid-flow/1d/euler.cu b/examples/inviscid-flow/1d/euler.cu deleted file mode 100644 index 4d76005cb1..0000000000 --- a/examples/inviscid-flow/1d/euler.cu +++ /dev/null @@ -1 +0,0 @@ -#include "euler.h" diff --git a/examples/inviscid-flow/1d/euler.h b/examples/inviscid-flow/1d/euler.h deleted file mode 100644 index 65bf558206..0000000000 --- a/examples/inviscid-flow/1d/euler.h +++ /dev/null @@ -1,119 +0,0 @@ -#include <TNL/tnlConfig.h> -#include <TNL/Solvers/Solver.h> -#include <TNL/Solvers/BuildConfigTags.h> -#include <TNL/Operators/DirichletBoundaryConditions.h> -#include <TNL/Operators/NeumannBoundaryConditions.h> -#include <TNL/Functions/Analytic/Constant.h> -#include "eulerProblem.h" -#include "LaxFridrichs.h" - -#include "eulerRhs.h" -#include "eulerBuildConfigTag.h" - -using namespace TNL; - -typedef eulerBuildConfigTag BuildConfig; - -/**** - * Uncoment the following (and comment the previous line) for the complete build. - * This will include support for all floating point precisions, all indexing types - * and more solvers. You may then choose between them from the command line. - * The compile time may, however, take tens of minutes or even several hours, - * esppecially if CUDA is enabled. Use this, if you want, only for the final build, - * not in the development phase. - */ -//typedef tnlDefaultConfigTag BuildConfig; - -template< typename ConfigTag >class eulerConfig -{ - public: - static void configSetup( Config::ConfigDescription & config ) - { - config.addDelimiter( "euler settings:" ); - config.addEntry< String >( "boundary-conditions-type", "Choose the boundary conditions type.", "dirichlet"); - config.addEntryEnum< String >( "dirichlet" ); - config.addEntryEnum< String >( "neumann" ); - config.addEntry< double >( "boundary-conditions-constant", "This sets a value in case of the constant boundary conditions." ); - config.addEntry< double >( "left-density", "This sets a value of left density." ); - config.addEntry< double >( "left-velocity", "This sets a value of left velocity." ); - config.addEntry< double >( "left-pressure", "This sets a value of left pressure." ); - config.addEntry< double >( "riemann-border", "This sets a position of discontinuity." ); - config.addEntry< double >( "right-density", "This sets a value of right density." ); - config.addEntry< double >( "right-velocity", "This sets a value of right velocity." ); - config.addEntry< double >( "right-pressure", "This sets a value of right pressure." ); - config.addEntry< double >( "gamma", "This sets a value of gamma constant." ); - - /**** - * Add definition of your solver command line arguments. - */ - - } -}; - -template< typename Real, - typename Device, - typename Index, - typename MeshType, - typename ConfigTag, - typename SolverStarter > -class eulerSetter -{ - public: - - typedef Real RealType; - typedef Device DeviceType; - typedef Index IndexType; - - static bool run( const Config::ParameterContainer & parameters ) - { - enum { Dimensions = MeshType::getMeshDimensions() }; - typedef LaxFridrichs< MeshType, Real, Index > ApproximateOperator; - typedef eulerRhs< MeshType, Real > RightHandSide; - typedef Containers::StaticVector < MeshType::getMeshDimensions(), Real > Vertex; - - /**** - * Resolve the template arguments of your solver here. - * The following code is for the Dirichlet and the Neumann boundary conditions. - * Both can be constant or defined as descrete values of Vector. - */ - String boundaryConditionsType = parameters.getParameter< String >( "boundary-conditions-type" ); - if( parameters.checkParameter( "boundary-conditions-constant" ) ) - { - typedef Functions::Analytic::Constant< Dimensions, Real > Constant; - if( boundaryConditionsType == "dirichlet" ) - { - typedef Operators::DirichletBoundaryConditions< MeshType, Constant, MeshType::getMeshDimensions(), Real, Index > BoundaryConditions; - typedef eulerProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem; - SolverStarter solverStarter; - return solverStarter.template run< Problem >( parameters ); - } - typedef Operators::NeumannBoundaryConditions< MeshType, Constant, Real, Index > BoundaryConditions; - typedef eulerProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem; - SolverStarter solverStarter; - return solverStarter.template run< Problem >( parameters ); - } - typedef Functions::MeshFunction< MeshType > MeshFunction; - if( boundaryConditionsType == "dirichlet" ) - { - typedef Operators::DirichletBoundaryConditions< MeshType, MeshFunction, MeshType::getMeshDimensions(), Real, Index > BoundaryConditions; - typedef eulerProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem; - SolverStarter solverStarter; - return solverStarter.template run< Problem >( parameters ); - } - typedef Operators::NeumannBoundaryConditions< MeshType, MeshFunction, Real, Index > BoundaryConditions; - typedef eulerProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem; - SolverStarter solverStarter; - return solverStarter.template run< Problem >( parameters ); - } - -}; - -int main( int argc, char* argv[] ) -{ - Solvers::Solver< eulerSetter, eulerConfig, BuildConfig > solver; - if( ! solver. run( argc, argv ) ) - return EXIT_FAILURE; - return EXIT_SUCCESS; -} - - diff --git a/examples/inviscid-flow/1d/eulerBuildConfigTag.h b/examples/inviscid-flow/1d/eulerBuildConfigTag.h deleted file mode 100644 index abede16343..0000000000 --- a/examples/inviscid-flow/1d/eulerBuildConfigTag.h +++ /dev/null @@ -1,51 +0,0 @@ -#ifndef eulerBUILDCONFIGTAG_H_ -#define eulerBUILDCONFIGTAG_H_ - -#include <TNL/Solvers/BuildConfigTags.h> - -namespace TNL { - -class eulerBuildConfigTag{}; - -namespace Solvers { - -/**** - * Turn off support for float and long double. - */ -template<> struct ConfigTagReal< eulerBuildConfigTag, float > { enum { enabled = false }; }; -template<> struct ConfigTagReal< eulerBuildConfigTag, long double > { enum { enabled = false }; }; - -/**** - * Turn off support for short int and long int indexing. - */ -template<> struct ConfigTagIndex< eulerBuildConfigTag, short int >{ enum { enabled = false }; }; -template<> struct ConfigTagIndex< eulerBuildConfigTag, long int >{ enum { enabled = false }; }; - -template< int Dimensions > struct ConfigTagDimensions< eulerBuildConfigTag, Dimensions >{ enum { enabled = ( Dimensions == 1 ) }; }; - -/**** - * Use of Grid is enabled for allowed dimensions and Real, Device and Index types. - */ -template< int Dimensions, typename Real, typename Device, typename Index > - struct ConfigTagMesh< eulerBuildConfigTag, Meshes::Grid< Dimensions, Real, Device, Index > > - { enum { enabled = ConfigTagDimensions< eulerBuildConfigTag, Dimensions >::enabled && - ConfigTagReal< eulerBuildConfigTag, Real >::enabled && - ConfigTagDevice< eulerBuildConfigTag, Device >::enabled && - ConfigTagIndex< eulerBuildConfigTag, Index >::enabled }; }; - -/**** - * Please, chose your preferred time discretisation here. - */ -template<> struct ConfigTagTimeDiscretisation< eulerBuildConfigTag, ExplicitTimeDiscretisationTag >{ enum { enabled = true }; }; -template<> struct ConfigTagTimeDiscretisation< eulerBuildConfigTag, SemiImplicitTimeDiscretisationTag >{ enum { enabled = false }; }; -template<> struct ConfigTagTimeDiscretisation< eulerBuildConfigTag, ImplicitTimeDiscretisationTag >{ enum { enabled = false }; }; - -/**** - * Only the Runge-Kutta-Merson solver is enabled by default. - */ -template<> struct ConfigTagExplicitSolver< eulerBuildConfigTag, ExplicitEulerSolverTag >{ enum { enabled = true }; }; - -} // namespace Solvers -} // namespace TNL - -#endif /* eulerBUILDCONFIGTAG_H_ */ diff --git a/examples/inviscid-flow/1d/eulerProblem.h b/examples/inviscid-flow/1d/eulerProblem.h deleted file mode 100644 index f687ee5460..0000000000 --- a/examples/inviscid-flow/1d/eulerProblem.h +++ /dev/null @@ -1,120 +0,0 @@ -#ifndef eulerPROBLEM_H_ -#define eulerPROBLEM_H_ - -#include <TNL/Problems/PDEProblem.h> -#include <TNL/Functions/MeshFunction.h> - -using namespace TNL::Problems; -namespace TNL { - -template< typename Mesh, - typename BoundaryCondition, - typename RightHandSide, - typename DifferentialOperator > -class eulerProblem: - public PDEProblem< Mesh, - typename DifferentialOperator::RealType, - typename Mesh::DeviceType, - typename DifferentialOperator::IndexType > -{ - public: - - typedef typename DifferentialOperator::RealType RealType; - typedef typename Mesh::DeviceType DeviceType; - typedef typename DifferentialOperator::IndexType IndexType; - typedef Functions::MeshFunction< Mesh > MeshFunctionType; - typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer; - typedef SharedPointer< DifferentialOperator > DifferentialOperatorPointer; - typedef SharedPointer< BoundaryCondition > BoundaryConditionPointer; - typedef SharedPointer< RightHandSide, DeviceType > RightHandSidePointer; - typedef PDEProblem< Mesh, RealType, DeviceType, IndexType > BaseType; - - using typename BaseType::MeshType; - using typename BaseType::MeshPointer; - using typename BaseType::DofVectorType; - using typename BaseType::DofVectorPointer; - using typename BaseType::MeshDependentDataType; - using typename BaseType::MeshDependentDataPointer; - - typedef typename DifferentialOperator::Continuity Continuity; - typedef typename DifferentialOperator::Momentum Momentum; - typedef typename DifferentialOperator::Energy Energy; - typedef typename DifferentialOperator::Velocity Velocity; - typedef typename DifferentialOperator::Pressure Pressure; - - - - static String getTypeStatic(); - - String getPrologHeader() const; - - void writeProlog( Logger& logger, - const Config::ParameterContainer& parameters ) const; - - bool setup( const MeshPointer& meshPointer, - const Config::ParameterContainer& parameters, - const String& prefix ); - - bool setInitialCondition( const Config::ParameterContainer& parameters, - const MeshPointer& mesh, - DofVectorPointer& dofs, - MeshDependentDataPointer& meshDependentData ); - - template< typename Matrix > - bool setupLinearSystem( const MeshPointer& mesh, - Matrix& matrix ); - - bool makeSnapshot( const RealType& time, - const IndexType& step, - const MeshPointer& mesh, - DofVectorPointer& dofs, - MeshDependentDataPointer& meshDependentData ); - - IndexType getDofs( const MeshPointer& mesh ) const; - - void bindDofs( const MeshPointer& mesh, - DofVectorPointer& dofs ); - - void getExplicitRHS( const RealType& time, - const RealType& tau, - const MeshPointer& mesh, - DofVectorPointer& _u, - DofVectorPointer& _fu, - MeshDependentDataPointer& meshDependentData ); - - template< typename Matrix > - void assemblyLinearSystem( const RealType& time, - const RealType& tau, - const MeshPointer& mesh, - DofVectorPointer& dofs, - Matrix& matrix, - DofVectorPointer& rightHandSide, - MeshDependentDataPointer& meshDependentData ); - - bool postIterate( const RealType& time, - const RealType& tau, - const MeshPointer& mesh, - DofVectorPointer& dofs, - MeshDependentDataPointer& meshDependentData ); - - protected: - - DifferentialOperatorPointer differentialOperatorPointer; - BoundaryConditionPointer boundaryConditionsPointer; - RightHandSidePointer rightHandSidePointer; - - MeshFunctionPointer uRho, uRhoVelocity, uEnergy; - MeshFunctionPointer fuRho, fuRhoVelocity, fuEnergy; - - MeshFunctionPointer pressure, velocity, rho, rhoVel, energy; - - RealType gamma; - - -}; - -} // namepsace TNL - -#include "eulerProblem_impl.h" - -#endif /* eulerPROBLEM_H_ */ diff --git a/examples/inviscid-flow/1d/eulerProblem_impl.h b/examples/inviscid-flow/1d/eulerProblem_impl.h deleted file mode 100644 index 0ed9d4f1a5..0000000000 --- a/examples/inviscid-flow/1d/eulerProblem_impl.h +++ /dev/null @@ -1,383 +0,0 @@ -#ifndef eulerPROBLEM_IMPL_H_ -#define eulerPROBLEM_IMPL_H_ - -#include <TNL/FileName.h> -#include <TNL/Matrices/MatrixSetter.h> -#include <TNL/Solvers/PDE/ExplicitUpdater.h> -#include <TNL/Solvers/PDE/LinearSystemAssembler.h> -#include <TNL/Solvers/PDE/BackwardTimeDiscretisation.h> -#include "LaxFridrichsContinuity.h" -#include "LaxFridrichsMomentum.h" -#include "LaxFridrichsEnergy.h" -#include "EulerVelGetter.h" -#include "EulerPressureGetter.h" - -namespace TNL { - -template< typename Mesh, - typename BoundaryCondition, - typename RightHandSide, - typename DifferentialOperator > -String -eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >:: -getTypeStatic() -{ - return String( "eulerProblem< " ) + Mesh :: getTypeStatic() + " >"; -} - -template< typename Mesh, - typename BoundaryCondition, - typename RightHandSide, - typename DifferentialOperator > -String -eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >:: -getPrologHeader() const -{ - return String( "euler" ); -} - -template< typename Mesh, - typename BoundaryCondition, - typename RightHandSide, - typename DifferentialOperator > -void -eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >:: -writeProlog( Logger& logger, const Config::ParameterContainer& parameters ) const -{ - /**** - * Add data you want to have in the computation report (log) as follows: - * logger.writeParameter< double >( "Parameter description", parameter ); - */ -} - -template< typename Mesh, - typename BoundaryCondition, - typename RightHandSide, - typename DifferentialOperator > -bool -eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >:: -setup( const MeshPointer& meshPointer, - const Config::ParameterContainer& parameters, - const String& prefix ) -{ - if( ! this->boundaryConditionsPointer->setup( meshPointer, parameters, prefix + "boundary-conditions-" ) || - ! this->rightHandSidePointer->setup( parameters, prefix + "right-hand-side-" ) ) - return false; - return true; -} - -template< typename Mesh, - typename BoundaryCondition, - typename RightHandSide, - typename DifferentialOperator > -typename eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::IndexType -eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >:: -getDofs( const MeshPointer& mesh ) const -{ - /**** - * Return number of DOFs (degrees of freedom) i.e. number - * of unknowns to be resolved by the main solver. - */ - return 3*mesh->template getEntitiesCount< typename MeshType::Cell >(); -} - -template< typename Mesh, - typename BoundaryCondition, - typename RightHandSide, - typename DifferentialOperator > -void -eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >:: -bindDofs( const MeshPointer& mesh, - DofVectorPointer& dofVector ) -{ -} - -template< typename Mesh, - typename BoundaryCondition, - typename RightHandSide, - typename DifferentialOperator > -bool -eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >:: -setInitialCondition( const Config::ParameterContainer& parameters, - const MeshPointer& mesh, - DofVectorPointer& dofs, - MeshDependentDataPointer& meshDependentData ) -{ - std::cout << std::endl << "get conditions from CML"; - typedef typename MeshType::Cell Cell; - this->gamma = parameters.getParameter< RealType >( "gamma" ); - RealType rhoL = parameters.getParameter< RealType >( "left-density" ); - RealType velL = parameters.getParameter< RealType >( "left-velocity" ); - RealType preL = parameters.getParameter< RealType >( "left-pressure" ); - RealType eL = ( preL / (gamma - 1) ) + 0.5 * rhoL * velL * velL; - RealType rhoR = parameters.getParameter< RealType >( "right-density" ); - RealType velR = parameters.getParameter< RealType >( "right-velocity" ); - RealType preR = parameters.getParameter< RealType >( "right-pressure" ); - RealType eR = ( preR / (gamma - 1) ) + 0.5 * rhoR * velR * velR; - RealType x0 = parameters.getParameter< RealType >( "riemann-border" ); - std::cout << std::endl << gamma << " " << rhoL << " " << velL << " " << preL << " " << eL << " " << rhoR << " " << velR << " " << preR << " " << eR << " " << x0 << " " << gamma << std::endl; - int count = mesh->template getEntitiesCount< Cell >(); - std::cout << count << std::endl; - uRho->bind( mesh, *dofs, 0); - uRhoVelocity->bind( mesh, *dofs, count); - uEnergy->bind( mesh, *dofs, 2 * count); - Containers::Vector < RealType, DeviceType, IndexType > data; - data.setSize(2*count); - velocity->bind( mesh, data, 0); - pressure->bind( mesh, data, count ); - std::cout << std::endl << "set conditions from CML"<< std::endl; - for(IndexType i = 0; i < count; i++) - if (i < x0 * count ) - { - ( *uRho )[i] = rhoL; - ( *uRhoVelocity )[i] = rhoL * velL; - ( *uEnergy )[i] = eL; - ( *velocity )[i] = velL; - ( *pressure )[i] = preL; - } - else - { - ( *uRho )[i] = rhoR; - ( *uRhoVelocity )[i] = rhoR * velR; - ( *uEnergy )[i] = eR; - ( *velocity )[i] = velR; - ( *pressure )[i] = preR; - }; - std::cout << "dofs = " << *dofs << std::endl; - getchar(); - - - /* - const String& initialConditionFile = parameters.getParameter< String >( "initial-condition" ); - if( ! dofs.load( initialConditionFile ) ) - { - std::cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." << std::endl; - return false; - } - */ - return true; -} -template< typename Mesh, - typename BoundaryCondition, - typename RightHandSide, - typename DifferentialOperator > - template< typename Matrix > -bool -eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >:: -setupLinearSystem( const MeshPointer& mesh, - Matrix& matrix ) -{ -/* const IndexType dofs = this->getDofs( mesh ); - typedef typename Matrix::CompressedRowsLengthsVector CompressedRowsLengthsVectorType; - CompressedRowsLengthsVectorType rowLengths; - if( ! rowLengths.setSize( dofs ) ) - return false; - MatrixSetter< MeshType, DifferentialOperator, BoundaryCondition, CompressedRowsLengthsVectorType > matrixSetter; - matrixSetter.template getCompressedRowsLengths< typename Mesh::Cell >( mesh, - differentialOperator, - boundaryCondition, - rowLengths ); - matrix.setDimensions( dofs, dofs ); - if( ! matrix.setCompressedRowsLengths( rowLengths ) ) - return false;*/ - return true; -} - -template< typename Mesh, - typename BoundaryCondition, - typename RightHandSide, - typename DifferentialOperator > -bool -eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >:: -makeSnapshot( const RealType& time, - const IndexType& step, - const MeshPointer& mesh, - DofVectorPointer& dofs, - MeshDependentDataPointer& meshDependentData ) -{ - std::cout << std::endl << "Writing output at time " << time << " step " << step << "." << std::endl; - this->bindDofs( mesh, dofs ); - typedef typename MeshType::Cell Cell; - int count = mesh->template getEntitiesCount< Cell >(); - std::ofstream vysledek; -/* std::cout << "pressure:" << std::endl; - for (IndexType i = 0; i<count; i++)std::cout << this->pressure[i] << " " << i ; - vysledek.open("pressure" + to_string(step) + ".txt"); - for (IndexType i = 0; i<count; i++) - vysledek << 0.01*i << " " << pressure[i] << std::endl; - vysledek.close(); - std::cout << " " << std::endl; - std::cout << "velocity:" << std::endl; - for (IndexType i = 0; i<count; i++)std::cout << this->velocity[i] << " " ; - vysledek.open("velocity" + to_string(step) + ".txt"); - for (IndexType i = 0; i<count; i++) - vysledek << 0.01*i << " " << pressure[i] << std::endl; - vysledek.close(); - std::cout << "energy:" << std::endl; - for (IndexType i = 0; i<count; i++)std::cout << this->uEnergy[i] << " " ; - vysledek.open("energy" + to_string(step) + ".txt"); - for (IndexType i = 0; i<count; i++) - vysledek << 0.01*i << " " << uEnergy[i] << std::endl; - vysledek.close(); - std::cout << " " << std::endl; - std::cout << "density:" << std::endl; - for (IndexType i = 0; i<count; i++)std::cout << this->uRho[i] << " " ; - vysledek.open("density" + to_string(step) + ".txt"); - for (IndexType i = 0; i<count; i++) - vysledek << 0.01*i << " " << uRho[i] << std::endl; - vysledek.close(); -*/ getchar(); - - FileName fileName; - fileName.setExtension( "tnl" ); - fileName.setIndex( step ); - fileName.setFileNameBase( "rho-" ); - - if( ! uRho->save( fileName.getFileName() ) ) - return false; - fileName.setFileNameBase( "rhoVel-" ); - if( ! uRhoVelocity->save( fileName.getFileName() ) ) - return false; - fileName.setFileNameBase( "energy-" ); - if( ! uEnergy->save( fileName.getFileName() ) ) - return false; - return true; -} - -template< typename Mesh, - typename BoundaryCondition, - typename RightHandSide, - typename DifferentialOperator > -void -eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >:: -getExplicitRHS( const RealType& time, - const RealType& tau, - const MeshPointer& mesh, - DofVectorPointer& _u, - DofVectorPointer& _fu, - MeshDependentDataPointer& meshDependentData ) -{ - std::cout << "explicitRHS" << std::endl; - typedef typename MeshType::Cell Cell; - int count = mesh->template getEntitiesCount< Cell >(); - //bind _u - this->uRho->bind(mesh, _u, 0); - this->uRhoVelocity->bind(mesh, _u ,count); - this->uEnergy->bind(mesh, _u, 2 * count); - - //bind _fu - this->fuRho->bind(mesh, _u, 0); - this->fuRhoVelocity->bind(mesh, _u, count); - this->fuEnergy->bind(mesh, _u, 2 * count); - - //generating Differential operator object - SharedPointer< Continuity > lF1DContinuity; - SharedPointer< Momentum > lF1DMomentum; - SharedPointer< Energy > lF1DEnergy; - - - - std::cout << "explicitRHSrho" << std::endl; - //rho - this->bindDofs( mesh, _u ); - lF1DContinuity->setTau(tau); - lF1DContinuity->setVelocity( *velocity); - /*ExplicitUpdater< Mesh, MeshFunctionType, Continuity, BoundaryCondition, RightHandSide > explicitUpdaterContinuity; - explicitUpdaterContinuity.template update< typename Mesh::Cell >( time, - mesh, - lF1DContinuity, - this->boundaryCondition, - this->rightHandSide, - uRho, - fuRho );*/ - - std::cout << "explicitRHSrhovel" << std::endl; - //rhoVelocity - lF1DMomentum->setTau(tau); - lF1DMomentum->setVelocity( *velocity); - lF1DMomentum->setPressure( *pressure); - Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, Momentum, BoundaryCondition, RightHandSide > explicitUpdaterMomentum; - explicitUpdaterMomentum.template update< typename Mesh::Cell >( time, - mesh, - lF1DMomentum, - this->boundaryConditionsPointer, - this->rightHandSidePointer, - uRhoVelocity, - fuRhoVelocity ); - - std::cout << "explicitRHSenergy" << std::endl; - //energy - lF1DEnergy->setTau(tau); - lF1DEnergy->setPressure( *pressure); - lF1DEnergy->setVelocity( *velocity); - Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, Energy, BoundaryCondition, RightHandSide > explicitUpdaterEnergy; - explicitUpdaterEnergy.template update< typename Mesh::Cell >( time, - mesh, - lF1DEnergy, - this->boundaryConditionsPointer, - this->rightHandSidePointer, - uEnergy, - fuEnergy ); - - } - -template< typename Mesh, - typename BoundaryCondition, - typename RightHandSide, - typename DifferentialOperator > - template< typename Matrix > -void -eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >:: -assemblyLinearSystem( const RealType& time, - const RealType& tau, - const MeshPointer& mesh, - DofVectorPointer& _u, - Matrix& matrix, - DofVectorPointer& b, - MeshDependentDataPointer& meshDependentData ) -{ -/* LinearSystemAssembler< Mesh, - MeshFunctionType, - DifferentialOperator, - BoundaryCondition, - RightHandSide, - BackwardTimeDiscretisation, - Matrix, - DofVectorType > systemAssembler; - - MeshFunction< Mesh > u( mesh, _u ); - systemAssembler.template assembly< typename Mesh::Cell >( time, - tau, - mesh, - this->differentialOperator, - this->boundaryCondition, - this->rightHandSide, - u, - matrix, - b );*/ -} -template< typename Mesh, - typename BoundaryCondition, - typename RightHandSide, - typename DifferentialOperator > -bool -eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >:: -postIterate( const RealType& time, - const RealType& tau, - const MeshPointer& mesh, - DofVectorPointer& dofs, - MeshDependentDataPointer& meshDependentData ) -{ - //velocity - this->velocity->setMesh( mesh ); - Velocity velocityGetter( *uRho, *uRhoVelocity ); - *this->velocity = velocityGetter; - //pressure - this->pressure->setMesh( mesh ); - Pressure pressureGetter( *uRho, *uRhoVelocity, *uEnergy, gamma ); - *this->pressure = pressureGetter; -} - -} // namespace TNL - -#endif /* eulerPROBLEM_IMPL_H_ */ diff --git a/examples/inviscid-flow/1d/eulerRhs.h b/examples/inviscid-flow/1d/eulerRhs.h deleted file mode 100644 index 1b46dc831f..0000000000 --- a/examples/inviscid-flow/1d/eulerRhs.h +++ /dev/null @@ -1,35 +0,0 @@ -#ifndef eulerRHS_H_ -#define eulerRHS_H_ - -#include <TNL/Functions/Domain.h> - -namespace TNL { - -template< typename Mesh, typename Real >class eulerRhs - : public Functions::Domain< Mesh::meshDimensions, Functions::MeshDomain > - { - public: - - typedef Mesh MeshType; - typedef Real RealType; - - bool setup( const Config::ParameterContainer& parameters, - const String& prefix = "" ) - { - return true; - } - - template< typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshEntity& entity, - const Real& time = 0.0 ) const - { - typedef typename MeshEntity::MeshType::VertexType VertexType; - VertexType v = entity.getCenter(); - return 0.0; - } -}; - -} //namespace TNL - -#endif /* eulerRHS_H_ */ diff --git a/examples/inviscid-flow/1d/tnl-run-euler-1d b/examples/inviscid-flow/1d/tnl-run-euler-1d deleted file mode 100644 index a45f664042..0000000000 --- a/examples/inviscid-flow/1d/tnl-run-euler-1d +++ /dev/null @@ -1,25 +0,0 @@ -#!/usr/bin/env bash - -tnl-grid-setup --dimensions 1 \ - --origin-x 0.0 \ - --proportions-x 1.0 \ - --size-x 100 - -#tnl-init --test-function sin-wave \ -# --output-file init.tnl -tnl-euler-1d-dbg --time-discretisation explicit \ - --time-step 1.0e-3 \ - --boundary-conditions-constant 0 \ - --discrete-solver euler \ - --snapshot-period 0.1 \ - --final-time 1.0 \ - --left-density 1.0 \ - --left-velocity 0.75 \ - --left-pressure 1.0 \ - --right-density 0.125 \ - --right-velocity 0 \ - --right-pressure 0.1 \ - --gamma 1.4 \ - --riemann-border 0.3 \ - -tnl-view --mesh mesh.tnl --input-files *tnl diff --git a/examples/inviscid-flow/CMakeLists.txt b/examples/inviscid-flow/CMakeLists.txt index 98bb2139da..9eb17cce74 100644 --- a/examples/inviscid-flow/CMakeLists.txt +++ b/examples/inviscid-flow/CMakeLists.txt @@ -1,4 +1,7 @@ -set( tnl_heat_equation_SOURCES +set( tnl_inviscid_flow_HEADERS + CompressibleConservativeVariables.h ) + +set( tnl_inviscid_flow_SOURCES euler.cpp euler.cu ) @@ -16,6 +19,6 @@ INSTALL( TARGETS tnl-euler-2d${debugExt} PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE ) INSTALL( FILES run-euler - ${tnl_heat_equation_SOURCES} + ${tnl_inviscid_flow_SOURCES} DESTINATION share/tnl-${tnlVersion}/examples/inviscid-flow-2d ) diff --git a/examples/inviscid-flow/EulerPressureGetter.h b/examples/inviscid-flow/EulerPressureGetter.h deleted file mode 100644 index 49e50d7257..0000000000 --- a/examples/inviscid-flow/EulerPressureGetter.h +++ /dev/null @@ -1,222 +0,0 @@ -#ifndef EulerPressureGetter_H -#define EulerPressureGetter_H - -#include <TNL/Containers/Vector.h> -#include <TNL/Meshes/Grid.h> - -namespace TNL { - -template< typename Mesh, - typename Real = typename Mesh::RealType, - typename Index = typename Mesh::IndexType > -class EulerPressureGetter -{ -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerPressureGetter< Meshes::Grid< 1,MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real gamma; - MeshFunctionType velocity; - MeshFunctionType rho; - MeshFunctionType energy; - - void setGamma(const Real& gamma) - { - this->gamma = gamma; - }; - - void setVelocity(const MeshFunctionType& velocity) - { - this->velocity = velocity; - }; - - void setRho(const MeshFunctionType& rho) - { - this->rho = rho; - }; - - void setEnergy(const MeshFunctionType& energy) - { - this->energy = energy; - }; - - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerPressureGetter< Meshes::Grid< 2,MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real gamma; - MeshFunctionType velocity; - MeshFunctionType rho; - MeshFunctionType energy; - - void setGamma(const Real& gamma) - { - this->gamma = gamma; - }; - - void setVelocity(const MeshFunctionType& velocity) - { - this->velocity = velocity; - }; - - void setRho(const MeshFunctionType& rho) - { - this->rho = rho; - }; - - void setEnergy(const MeshFunctionType& energy) - { - this->energy = energy; - }; - - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerPressureGetter< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real gamma; - MeshFunctionType velocity; - MeshFunctionType rho; - MeshFunctionType energy; - - void setGamma(const Real& gamma) - { - this->gamma = gamma; - }; - - void setVelocity(const MeshFunctionType& velocity) - { - this->velocity = velocity; - }; - - void setRho(const MeshFunctionType& rho) - { - this->rho = rho; - }; - - void setEnergy(const MeshFunctionType& energy) - { - this->energy = energy; - }; - - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - - -} //namespace TNL - -#include "EulerPressureGetter_impl.h" - -#endif /* EulerPressureGetter_H */ diff --git a/examples/inviscid-flow/EulerPressureGetter_impl.h b/examples/inviscid-flow/EulerPressureGetter_impl.h deleted file mode 100644 index 7d6b9310df..0000000000 --- a/examples/inviscid-flow/EulerPressureGetter_impl.h +++ /dev/null @@ -1,338 +0,0 @@ -#ifndef EulerPressureGetter_IMPL_H -#define EulerPressureGetter_IMPL_H - -namespace TNL { - -/**** - * 1D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -EulerPressureGetter< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "EulerPressureGetter< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -EulerPressureGetter< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 1, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 1, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& hxSquareInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -EulerPressureGetter< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -EulerPressureGetter< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - matrixRow.setElement( 0, west, - lambdaX ); - matrixRow.setElement( 1, center, 2.0 * lambdaX ); - matrixRow.setElement( 2, east, - lambdaX ); -} - -/**** - * 2D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -EulerPressureGetter< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "EulerPressureGetter< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -EulerPressureGetter< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - //pressure - const IndexType& center = entity.getIndex(); - - return ( this->gamma - 1 ) * ( energy[ center ] - 0.5 * rho[ center ] * ::pow( velocity[ center ],2 )); -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -EulerPressureGetter< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -EulerPressureGetter< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2, 0 >(); - const RealType& lambdaY = tau * entity.getMesh().template getSpaceStepsProducts< 0, -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); - matrixRow.setElement( 0, south, -lambdaY ); - matrixRow.setElement( 1, west, -lambdaX ); - matrixRow.setElement( 2, center, 2.0 * ( lambdaX + lambdaY ) ); - matrixRow.setElement( 3, east, -lambdaX ); - matrixRow.setElement( 4, north, -lambdaY ); -} - -/**** - * 3D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -EulerPressureGetter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "EulerPressureGetter< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -EulerPressureGetter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 3, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 3, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse + - ( u[ south ] - 2.0 * u[ center ] + u[ north ] ) * hySquareInverse + - ( u[ up ] - 2.0 * u[ center ] + u[ down ] ) * hzSquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -EulerPressureGetter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -EulerPressureGetter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - matrixRow.setElement( 0, down, -lambdaZ ); - matrixRow.setElement( 1, south, -lambdaY ); - matrixRow.setElement( 2, west, -lambdaX ); - matrixRow.setElement( 3, center, 2.0 * ( lambdaX + lambdaY + lambdaZ ) ); - matrixRow.setElement( 4, east, -lambdaX ); - matrixRow.setElement( 5, north, -lambdaY ); - matrixRow.setElement( 6, up, -lambdaZ ); -} - -} // namespace TNL - -#endif /* EulerPressureGetterIMPL_H */ - diff --git a/examples/inviscid-flow/EulerVelGetter.h b/examples/inviscid-flow/EulerVelGetter.h deleted file mode 100644 index 2da0052b87..0000000000 --- a/examples/inviscid-flow/EulerVelGetter.h +++ /dev/null @@ -1,186 +0,0 @@ -#ifndef EulerVelGetter_H -#define EulerVelGetter_H - -#include <TNL/Containers/Vector.h> -#include <TNL/Meshes/Grid.h> -#include <TNL/Functions/Domain.h> - -namespace TNL { - -template< typename Mesh, - typename Real = typename Mesh::RealType, - typename Index = typename Mesh::IndexType > -class EulerVelGetter -{ -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerVelGetter< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index > - : public TNL::Functions::Domain< 1, TNL::Functions::MeshDomain > -{ - public: - typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType; - typedef typename MeshType::CoordinatesType CoordinatesType; - typedef Real RealType; - typedef Device DeviceType; - typedef Index IndexType; - typedef Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - MeshFunctionType velX; - MeshFunctionType velY; - - void setVelX(const MeshFunctionType& velX) - { - this->velX = velX; - }; - - void setVelY(const MeshFunctionType& velY) - { - this->velY = velY; - }; - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerVelGetter< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index > - : public TNL::Functions::Domain< 2, TNL::Functions::MeshDomain > -{ - public: - typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType; - typedef typename MeshType::CoordinatesType CoordinatesType; - typedef Real RealType; - typedef Device DeviceType; - typedef Index IndexType; - typedef Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - MeshFunctionType velX; - MeshFunctionType velY; - - void setVelX(const MeshFunctionType& velX) - { - this->velX = velX; - }; - - void setVelY(const MeshFunctionType& velY) - { - this->velY = velY; - }; - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerVelGetter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index > - : public TNL::Functions::Domain< 3, TNL::Functions::MeshDomain > -{ - public: - typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType; - typedef typename MeshType::CoordinatesType CoordinatesType; - typedef Real RealType; - typedef Device DeviceType; - typedef Index IndexType; - typedef Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - MeshFunctionType velX; - MeshFunctionType velY; - - void setVelX(const MeshFunctionType& velX) - { - this->velX = velX; - }; - - void setVelY(const MeshFunctionType& velY) - { - this->velY = velY; - }; - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -} //namespace TNL - -#include "EulerVelGetter_impl.h" - -#endif /* EulerVelGetter_H */ diff --git a/examples/inviscid-flow/EulerVelGetter_impl.h b/examples/inviscid-flow/EulerVelGetter_impl.h deleted file mode 100644 index eeac07e570..0000000000 --- a/examples/inviscid-flow/EulerVelGetter_impl.h +++ /dev/null @@ -1,340 +0,0 @@ -#ifndef EulerVelGetter_IMPL_H -#define EulerVelGetter_IMPL_H - -namespace TNL { - -/**** - * 1D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -EulerVelGetter< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "EulerVelGetter< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -EulerVelGetter< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 1, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 1, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& hxSquareInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -EulerVelGetter< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -EulerVelGetter< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - matrixRow.setElement( 0, west, - lambdaX ); - matrixRow.setElement( 1, center, 2.0 * lambdaX ); - matrixRow.setElement( 2, east, - lambdaX ); -} - -/**** - * 2D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -EulerVelGetter< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "EulerVelGetter< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -EulerVelGetter< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 2, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 2, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - //vel - const IndexType& center = entity.getIndex(); - return ::sqrt(pow( velX[ center ],2)+pow( velY[ center ],2)); -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -EulerVelGetter< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -EulerVelGetter< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2, 0 >(); - const RealType& lambdaY = tau * entity.getMesh().template getSpaceStepsProducts< 0, -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); - matrixRow.setElement( 0, south, -lambdaY ); - matrixRow.setElement( 1, west, -lambdaX ); - matrixRow.setElement( 2, center, 2.0 * ( lambdaX + lambdaY ) ); - matrixRow.setElement( 3, east, -lambdaX ); - matrixRow.setElement( 4, north, -lambdaY ); -} - -/**** - * 3D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -EulerVelGetter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "EulerVelGetter< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -EulerVelGetter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 3, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 3, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse + - ( u[ south ] - 2.0 * u[ center ] + u[ north ] ) * hySquareInverse + - ( u[ up ] - 2.0 * u[ center ] + u[ down ] ) * hzSquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -EulerVelGetter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -EulerVelGetter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - matrixRow.setElement( 0, down, -lambdaZ ); - matrixRow.setElement( 1, south, -lambdaY ); - matrixRow.setElement( 2, west, -lambdaX ); - matrixRow.setElement( 3, center, 2.0 * ( lambdaX + lambdaY + lambdaZ ) ); - matrixRow.setElement( 4, east, -lambdaX ); - matrixRow.setElement( 5, north, -lambdaY ); - matrixRow.setElement( 6, up, -lambdaZ ); -} - -} //namespace TNL - -#endif /* EulerVelGetterIMPL_H */ - diff --git a/examples/inviscid-flow/EulerVelXGetter.h b/examples/inviscid-flow/EulerVelXGetter.h deleted file mode 100644 index 93e3363bfa..0000000000 --- a/examples/inviscid-flow/EulerVelXGetter.h +++ /dev/null @@ -1,183 +0,0 @@ -#ifndef EulerVelXGetter_H -#define EulerVelXGetter_H - -#include <TNL/Containers/Vector.h> -#include <TNL/Meshes/Grid.h> - -namespace TNL { - -template< typename Mesh, - typename Real = typename Mesh::RealType, - typename Index = typename Mesh::IndexType > -class EulerVelXGetter -{ -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerVelXGetter< Meshes::Grid< 1,MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - MeshFunctionType rhoVelX; - MeshFunctionType rho; - - void setRhoVelX(const MeshFunctionType& rhoVelX) - { - this->rhoVelX = rhoVelX; - }; - - void setRho(const MeshFunctionType& rho) - { - this->rho = rho; - }; - - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerVelXGetter< Meshes::Grid< 2,MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - MeshFunctionType rhoVelX; - MeshFunctionType rho; - - void setRhoVelX(const MeshFunctionType& rhoVelX) - { - this->rhoVelX = rhoVelX; - }; - - void setRho(const MeshFunctionType& rho) - { - this->rho = rho; - }; - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerVelXGetter< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - MeshFunctionType rhoVelX; - MeshFunctionType rho; - - void setRhoVelX(const MeshFunctionType& rhoVelX) - { - this->rhoVelX = rhoVelX; - }; - - void setRho(const MeshFunctionType& rho) - { - this->rho = rho; - }; - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -} // namespace TNL - -#include "EulerVelXGetter_impl.h" - -#endif /* EulerVelXGetter_H */ diff --git a/examples/inviscid-flow/EulerVelXGetter_impl.h b/examples/inviscid-flow/EulerVelXGetter_impl.h deleted file mode 100644 index 6008509056..0000000000 --- a/examples/inviscid-flow/EulerVelXGetter_impl.h +++ /dev/null @@ -1,339 +0,0 @@ -#pragma once - -namespace TNL { - -/**** - * 1D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -EulerVelXGetter< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "EulerVelXGetter< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -EulerVelXGetter< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 1, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 1, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& hxSquareInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -EulerVelXGetter< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -EulerVelXGetter< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - matrixRow.setElement( 0, west, - lambdaX ); - matrixRow.setElement( 1, center, 2.0 * lambdaX ); - matrixRow.setElement( 2, east, - lambdaX ); -} - -/**** - * 2D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -EulerVelXGetter< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "EulerVelXGetter< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -EulerVelXGetter< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 2, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 2, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - //velX - const IndexType& center = entity.getIndex(); - return ( rhoVelX[ center ] / rho[ center ]); -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -EulerVelXGetter< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -EulerVelXGetter< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2, 0 >(); - const RealType& lambdaY = tau * entity.getMesh().template getSpaceStepsProducts< 0, -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); - matrixRow.setElement( 0, south, -lambdaY ); - matrixRow.setElement( 1, west, -lambdaX ); - matrixRow.setElement( 2, center, 2.0 * ( lambdaX + lambdaY ) ); - matrixRow.setElement( 3, east, -lambdaX ); - matrixRow.setElement( 4, north, -lambdaY ); -} - -/**** - * 3D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -EulerVelXGetter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "EulerVelXGetter< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -EulerVelXGetter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 3, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 3, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse + - ( u[ south ] - 2.0 * u[ center ] + u[ north ] ) * hySquareInverse + - ( u[ up ] - 2.0 * u[ center ] + u[ down ] ) * hzSquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -EulerVelXGetter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -EulerVelXGetter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - matrixRow.setElement( 0, down, -lambdaZ ); - matrixRow.setElement( 1, south, -lambdaY ); - matrixRow.setElement( 2, west, -lambdaX ); - matrixRow.setElement( 3, center, 2.0 * ( lambdaX + lambdaY + lambdaZ ) ); - matrixRow.setElement( 4, east, -lambdaX ); - matrixRow.setElement( 5, north, -lambdaY ); - matrixRow.setElement( 6, up, -lambdaZ ); -} - -}// namespace TNL - - - diff --git a/examples/inviscid-flow/EulerVelYGetter.h b/examples/inviscid-flow/EulerVelYGetter.h deleted file mode 100644 index 49474808d3..0000000000 --- a/examples/inviscid-flow/EulerVelYGetter.h +++ /dev/null @@ -1,183 +0,0 @@ -#ifndef EulerVelYGetter_H -#define EulerVelYGetter_H - -#include <TNL/Containers/Vector.h> -#include <TNL/Meshes/Grid.h> - -namespace TNL { - -template< typename Mesh, - typename Real = typename Mesh::RealType, - typename Index = typename Mesh::IndexType > -class EulerVelYGetter -{ -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerVelYGetter< Meshes::Grid< 1,MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - MeshFunctionType rhoVelY; - MeshFunctionType rho; - - void setRhoVelY(const MeshFunctionType& rhoVelY) - { - this->rhoVelY = rhoVelY; - }; - - void setRho(const MeshFunctionType& rho) - { - this->rho = rho; - }; - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerVelYGetter< Meshes::Grid< 2,MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - MeshFunctionType rhoVelY; - MeshFunctionType rho; - - void setRhoVelY(const MeshFunctionType& rhoVelY) - { - this->rhoVelY = rhoVelY; - }; - - void setRho(const MeshFunctionType& rho) - { - this->rho = rho; - }; - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class EulerVelYGetter< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - MeshFunctionType rhoVelY; - MeshFunctionType rho; - - void setRhoVelY(const MeshFunctionType& rhoVelY) - { - this->rhoVelY = rhoVelY; - }; - - void setRho(const MeshFunctionType& rho) - { - this->rho = rho; - }; - - template< typename MeshFunction, typename MeshEntity > - __cuda_callable__ - Real operator()( const MeshFunction& u, - const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > - __cuda_callable__ - Index getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const; - - template< typename MeshEntity, typename Vector, typename MatrixRow > - __cuda_callable__ - void updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const; -}; - -} //namespace TNL - - -#include "EulerVelYGetter_impl.h" - -#endif /* EulerVelYGetter_H */ diff --git a/examples/inviscid-flow/EulerVelYGetter_impl.h b/examples/inviscid-flow/EulerVelYGetter_impl.h deleted file mode 100644 index 337f4208c3..0000000000 --- a/examples/inviscid-flow/EulerVelYGetter_impl.h +++ /dev/null @@ -1,332 +0,0 @@ -#ifndef EulerVelYGetter_IMPL_H -#define EulerVelYGetter_IMPL_H - -namespace TNL { - -/**** - * 1D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -EulerVelYGetter< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "EulerVelYGetter< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -EulerVelYGetter< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 1, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 1, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& hxSquareInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -EulerVelYGetter< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -EulerVelYGetter< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - matrixRow.setElement( 0, west, - lambdaX ); - matrixRow.setElement( 1, center, 2.0 * lambdaX ); - matrixRow.setElement( 2, east, - lambdaX ); -} - -/**** - * 2D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -EulerVelYGetter< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "EulerVelYGetter< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -EulerVelYGetter< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - //velY - const IndexType& center = entity.getIndex(); - return ( rhoVelY[ center ] / rho[ center ]); -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -EulerVelYGetter< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -EulerVelYGetter< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2, 0 >(); - const RealType& lambdaY = tau * entity.getMesh().template getSpaceStepsProducts< 0, -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); - matrixRow.setElement( 0, south, -lambdaY ); - matrixRow.setElement( 1, west, -lambdaX ); - matrixRow.setElement( 2, center, 2.0 * ( lambdaX + lambdaY ) ); - matrixRow.setElement( 3, east, -lambdaX ); - matrixRow.setElement( 4, north, -lambdaY ); -} - -/**** - * 3D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -EulerVelYGetter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "EulerVelYGetter< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -EulerVelYGetter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 3, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 3, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse + - ( u[ south ] - 2.0 * u[ center ] + u[ north ] ) * hySquareInverse + - ( u[ up ] - 2.0 * u[ center ] + u[ down ] ) * hzSquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -EulerVelYGetter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -EulerVelYGetter< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - matrixRow.setElement( 0, down, -lambdaZ ); - matrixRow.setElement( 1, south, -lambdaY ); - matrixRow.setElement( 2, west, -lambdaX ); - matrixRow.setElement( 3, center, 2.0 * ( lambdaX + lambdaY + lambdaZ ) ); - matrixRow.setElement( 4, east, -lambdaX ); - matrixRow.setElement( 5, north, -lambdaY ); - matrixRow.setElement( 6, up, -lambdaZ ); -} - -} // namespace TNL - -#endif /* EulerVelYGetterIMPL_H */ - diff --git a/examples/inviscid-flow/LaxFridrichs.h b/examples/inviscid-flow/LaxFridrichs.h index b689e00212..549acba199 100644 --- a/examples/inviscid-flow/LaxFridrichs.h +++ b/examples/inviscid-flow/LaxFridrichs.h @@ -1,5 +1,15 @@ -#ifndef LaxFridrichs_H -#define LaxFridrichs_H +/*************************************************************************** + LaxFridrichs.h - description + ------------------- + begin : Feb 18, 2017 + copyright : (C) 2017 by Tomas Oberhuber + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/* See Copyright Notice in tnl/Copyright */ + + +#pragma once #include <TNL/Containers/Vector.h> #include <TNL/Meshes/Grid.h> @@ -8,10 +18,11 @@ #include "LaxFridrichsEnergy.h" #include "LaxFridrichsMomentumX.h" #include "LaxFridrichsMomentumY.h" -#include "EulerPressureGetter.h" +#include "LaxFridrichsMomentumZ.h" +/*#include "EulerPressureGetter.h" #include "EulerVelXGetter.h" #include "EulerVelYGetter.h" -#include "EulerVelGetter.h" +#include "EulerVelGetter.h"*/ namespace TNL { @@ -29,14 +40,13 @@ class LaxFridrichs typedef LaxFridrichsContinuity< Mesh, Real, Index > Continuity; typedef LaxFridrichsMomentumX< Mesh, Real, Index > MomentumX; typedef LaxFridrichsMomentumY< Mesh, Real, Index > MomentumY; + typedef LaxFridrichsMomentumZ< Mesh, Real, Index > MomentumZ; typedef LaxFridrichsEnergy< Mesh, Real, Index > Energy; - typedef EulerVelXGetter< Mesh, Real, Index > VelocityX; + /*typedef EulerVelXGetter< Mesh, Real, Index > VelocityX; typedef EulerVelYGetter< Mesh, Real, Index > VelocityY; typedef EulerVelGetter< Mesh, Real, Index > Velocity; - typedef EulerPressureGetter< Mesh, Real, Index > Pressure; + typedef EulerPressureGetter< Mesh, Real, Index > Pressure;*/ }; } //namespace TNL - -#endif /* LaxFridrichs_H */ diff --git a/examples/inviscid-flow/LaxFridrichsContinuity.h b/examples/inviscid-flow/LaxFridrichsContinuity.h index 813cdcd8a3..260bc9c6cd 100644 --- a/examples/inviscid-flow/LaxFridrichsContinuity.h +++ b/examples/inviscid-flow/LaxFridrichsContinuity.h @@ -1,8 +1,20 @@ -#ifndef LaxFridrichsContinuity_H -#define LaxFridrichsContinuity_H +/*************************************************************************** + LaxFridrichsContinuity.h - description + ------------------- + begin : Feb 17, 2017 + copyright : (C) 2017 by Tomas Oberhuber + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/* See Copyright Notice in tnl/Copyright */ + + +#pragma once #include <TNL/Containers/Vector.h> #include <TNL/Meshes/Grid.h> +#include <TNL/Functions/VectorField.h> +#include <TNL/SharedPointer.h> namespace TNL { @@ -13,6 +25,41 @@ template< typename Mesh, class LaxFridrichsContinuityBase { public: + + typedef Real RealType; + typedef Index IndexType; + typedef Mesh MeshType; + typedef typename MeshType::DeviceType DeviceType; + typedef typename MeshType::CoordinatesType CoordinatesType; + typedef Functions::MeshFunction< MeshType > MeshFunctionType; + static const int Dimensions = MeshType::getMeshDimensions(); + typedef Functions::VectorField< Dimensions, MeshFunctionType > VelocityFieldType; + typedef SharedPointer< VelocityFieldType > VelocityFieldPointer; + + static String getType() + { + return String( "LaxFridrichsContinuity< " ) + + MeshType::getType() + ", " + + TNL::getType< Real >() + ", " + + TNL::getType< Index >() + " >"; + } + + void setTau(const Real& tau) + { + this->tau = tau; + }; + + void setVelocity( const VelocityFieldPointer& velocity ) + { + this->velocity = velocity; + }; + + + protected: + + RealType tau; + + VelocityFieldPointer velocity; }; @@ -30,45 +77,43 @@ template< typename MeshReal, typename MeshIndex, typename Real, typename Index > -class LaxFridrichsContinuity< Meshes::Grid< 1,MeshReal, Device, MeshIndex >, Real, Index > +class LaxFridrichsContinuity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index > + : public LaxFridrichsContinuityBase< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocityX; - MeshFunctionType velocityY; - - void setTau(const Real& tau) - { - this->tau = tau; - }; - - void setVelocityX(const MeshFunctionType& velocityX) - { - this->velocityX = velocityX; - }; - - void setVelocityY(const MeshFunctionType& velocityY) - { - this->velocityY = velocityY; - }; - + typedef LaxFridrichsContinuityBase< MeshType, Real, Index > BaseType; + + using typename BaseType::RealType; + using typename BaseType::IndexType; + using typename BaseType::DeviceType; + using typename BaseType::CoordinatesType; + using typename BaseType::MeshFunctionType; + using typename BaseType::VelocityFieldType; + using typename BaseType::VelocityFieldPointer; + using typename BaseType::Dimensions; template< typename MeshFunction, typename MeshEntity > __cuda_callable__ Real operator()( const MeshFunction& u, const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > + const RealType& time = 0.0 ) const + { + static_assert( MeshEntity::entityDimensions == 1, "Wrong mesh entity dimensions." ); + static_assert( MeshFunction::getEntitiesDimensions() == 1, "Wrong preimage function" ); + const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); + + const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); + const IndexType& center = entity.getIndex(); + const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); + const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); + const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ]; + const RealType& velocity_x_east = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ]; + return (0.5 * this->tau) * ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) + - 0.5 * hxInverse * ( u[ west ] * velocity_x_west - u[ east ] * velocity_x_east ); + } + + /*template< typename MeshEntity > __cuda_callable__ Index getLinearSystemRowLength( const MeshType& mesh, const IndexType& index, @@ -83,7 +128,7 @@ class LaxFridrichsContinuity< Meshes::Grid< 1,MeshReal, Device, MeshIndex >, Rea const MeshEntity& entity, const MeshFunctionType& u, Vector& b, - MatrixRow& matrixRow ) const; + MatrixRow& matrixRow ) const;*/ }; template< typename MeshReal, @@ -91,45 +136,51 @@ template< typename MeshReal, typename MeshIndex, typename Real, typename Index > -class LaxFridrichsContinuity< Meshes::Grid< 2,MeshReal, Device, MeshIndex >, Real, Index > +class LaxFridrichsContinuity< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index > + : public LaxFridrichsContinuityBase< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocityX; // TODO: change to pointer - MeshFunctionType velocityY; - - void setTau(const Real& tau) - { - this->tau = tau; - }; - - void setVelocityX(const MeshFunctionType& velocityX) - { - this->velocityX = velocityX; - }; - - void setVelocityY(const MeshFunctionType& velocityY) - { - this->velocityY = velocityY; - }; - + typedef LaxFridrichsContinuityBase< MeshType, Real, Index > BaseType; + + using typename BaseType::RealType; + using typename BaseType::IndexType; + using typename BaseType::DeviceType; + using typename BaseType::CoordinatesType; + using typename BaseType::MeshFunctionType; + using typename BaseType::VelocityFieldType; + using typename BaseType::VelocityFieldPointer; + using BaseType::Dimensions; template< typename MeshFunction, typename MeshEntity > __cuda_callable__ Real operator()( const MeshFunction& u, const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > + const RealType& time = 0.0 ) const + { + static_assert( MeshEntity::entityDimensions == 2, "Wrong mesh entity dimensions." ); + static_assert( MeshFunction::getEntitiesDimensions() == 2, "Wrong preimage function" ); + const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); + + //rho + const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1, 0 >(); + const RealType& hyInverse = entity.getMesh().template getSpaceStepsProducts< 0, -1 >(); + const IndexType& center = entity.getIndex(); + const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); + const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); + const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); + const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); + const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ]; + const RealType& velocity_x_east = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ]; + const RealType& velocity_y_north = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ north ]; + const RealType& velocity_y_south = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ]; + + return 0.5 * ( this->tau * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] - 4.0 * u[ center ] ) + - ( u[ west ] * velocity_x_west - u[ east ] * velocity_x_east ) * hxInverse + - ( u[ north ] * velocity_y_north - u[ south ] * velocity_y_south ) * hyInverse ); + } + + /*template< typename MeshEntity > __cuda_callable__ Index getLinearSystemRowLength( const MeshType& mesh, const IndexType& index, @@ -144,7 +195,7 @@ class LaxFridrichsContinuity< Meshes::Grid< 2,MeshReal, Device, MeshIndex >, Rea const MeshEntity& entity, const MeshFunctionType& u, Vector& b, - MatrixRow& matrixRow ) const; + MatrixRow& matrixRow ) const;*/ }; template< typename MeshReal, @@ -152,45 +203,59 @@ template< typename MeshReal, typename MeshIndex, typename Real, typename Index > -class LaxFridrichsContinuity< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, Index > +class LaxFridrichsContinuity< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index > + : public LaxFridrichsContinuityBase< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocityX; - MeshFunctionType velocityY; - - void setTau(const Real& tau) - { - this->tau = tau; - }; - - void setVelocityX(const MeshFunctionType& velocityX) - { - this->velocityX = velocityX; - }; - - void setVelocityY(const MeshFunctionType& velocityY) - { - this->velocityY = velocityY; - }; - + typedef LaxFridrichsContinuityBase< MeshType, Real, Index > BaseType; + + using typename BaseType::RealType; + using typename BaseType::IndexType; + using typename BaseType::DeviceType; + using typename BaseType::Dimensions; + using typename BaseType::CoordinatesType; + using typename BaseType::MeshFunctionType; + using typename BaseType::VelocityFieldType; + using typename BaseType::VelocityFieldPointer; template< typename MeshFunction, typename MeshEntity > __cuda_callable__ Real operator()( const MeshFunction& u, const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > + const RealType& time = 0.0 ) const + { + static_assert( MeshEntity::entityDimensions == 3, "Wrong mesh entity dimensions." ); + static_assert( MeshFunction::getEntitiesDimensions() == 3, "Wrong preimage function" ); + const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); + + //rho + const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1, 0, 0 >(); + const RealType& hyInverse = entity.getMesh().template getSpaceStepsProducts< 0, -1, 0 >(); + const RealType& hzInverse = entity.getMesh().template getSpaceStepsProducts< 0, 0, -1 >(); + const IndexType& center = entity.getIndex(); + const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); + const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); + const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); + const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); + const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); + const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); + + const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ]; + const RealType& velocity_x_east = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ]; + const RealType& velocity_y_north = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ north ]; + const RealType& velocity_y_south = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ]; + const RealType& velocity_z_up = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ up ]; + const RealType& velocity_z_down = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ down ]; + + return 0.5 * ( this->tau * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] + u[ up ] + u[ down ]- 6.0 * u[ center ] ) + - ( u[ west ] * velocity_x_west - u[ east ] * velocity_x_east ) * hxInverse + - ( u[ north ] * velocity_y_north - u[ south ] * velocity_y_south ) * hyInverse + - ( u[ up ] * velocity_z_up - u[ down ] * velocity_z_down ) * hzInverse ); + + } + + /*template< typename MeshEntity > __cuda_callable__ Index getLinearSystemRowLength( const MeshType& mesh, const IndexType& index, @@ -205,12 +270,8 @@ class LaxFridrichsContinuity< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Rea const MeshEntity& entity, const MeshFunctionType& u, Vector& b, - MatrixRow& matrixRow ) const; + MatrixRow& matrixRow ) const;*/ }; } //namespace TNL - -#include "LaxFridrichsContinuity_impl .h" - -#endif /* LaxFridrichsContinuity_H */ diff --git a/examples/inviscid-flow/LaxFridrichsContinuity_impl .h b/examples/inviscid-flow/LaxFridrichsContinuity_impl .h deleted file mode 100644 index 237b1a5845..0000000000 --- a/examples/inviscid-flow/LaxFridrichsContinuity_impl .h +++ /dev/null @@ -1,349 +0,0 @@ -#ifndef LaxFridrichsContinuity_IMPL_H -#define LaxFridrichsContinuity_IMPL_H - -namespace TNL { - -/**** - * 1D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsContinuity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsContinuity< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsContinuity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 1, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 1, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& hxSquareInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsContinuity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsContinuity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - matrixRow.setElement( 0, west, - lambdaX ); - matrixRow.setElement( 1, center, 2.0 * lambdaX ); - matrixRow.setElement( 2, east, - lambdaX ); -} - -/**** - * 2D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsContinuity< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsContinuity< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsContinuity< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 2, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 2, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - - //rho - const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1, 0 >(); - const RealType& hyInverse = entity.getMesh().template getSpaceStepsProducts< 0, -1 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); - return (0.5 * this->tau) * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] - 4.0 * u[ center ] ) - - 0.5 * hxInverse * ( u[ west ] * this->velocityX[ west ] - u[ east ] * this->velocityX[ east ] ) - - 0.5 * hyInverse * ( u[ north ] * this->velocityY[ north ] - u[ south ] * this->velocityY[ east ] ); -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsContinuity< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsContinuity< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2, 0 >(); - const RealType& lambdaY = tau * entity.getMesh().template getSpaceStepsProducts< 0, -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); - matrixRow.setElement( 0, south, -lambdaY ); - matrixRow.setElement( 1, west, -lambdaX ); - matrixRow.setElement( 2, center, 2.0 * ( lambdaX + lambdaY ) ); - matrixRow.setElement( 3, east, -lambdaX ); - matrixRow.setElement( 4, north, -lambdaY ); -} - -/**** - * 3D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsContinuity< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsContinuity< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsContinuity< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 3, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 3, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse + - ( u[ south ] - 2.0 * u[ center ] + u[ north ] ) * hySquareInverse + - ( u[ up ] - 2.0 * u[ center ] + u[ down ] ) * hzSquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsContinuity< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsContinuity< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - matrixRow.setElement( 0, down, -lambdaZ ); - matrixRow.setElement( 1, south, -lambdaY ); - matrixRow.setElement( 2, west, -lambdaX ); - matrixRow.setElement( 3, center, 2.0 * ( lambdaX + lambdaY + lambdaZ ) ); - matrixRow.setElement( 4, east, -lambdaX ); - matrixRow.setElement( 5, north, -lambdaY ); - matrixRow.setElement( 6, up, -lambdaZ ); -} - -} //namespace TNL - -#endif /* LaxFridrichsContinuityIMPL_H */ - diff --git a/examples/inviscid-flow/LaxFridrichsEnergy.h b/examples/inviscid-flow/LaxFridrichsEnergy.h index 8685697295..7a351d1027 100644 --- a/examples/inviscid-flow/LaxFridrichsEnergy.h +++ b/examples/inviscid-flow/LaxFridrichsEnergy.h @@ -1,67 +1,125 @@ -#ifndef LaxFridrichsEnergy_H -#define LaxFridrichsEnergy_H +/*************************************************************************** + LaxFridrichsEnergy.h - description + ------------------- + begin : Feb 17, 2017 + copyright : (C) 2017 by Tomas Oberhuber + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/* See Copyright Notice in tnl/Copyright */ + +#pragma once #include <TNL/Containers/Vector.h> #include <TNL/Meshes/Grid.h> namespace TNL { - + template< typename Mesh, typename Real = typename Mesh::RealType, typename Index = typename Mesh::IndexType > -class LaxFridrichsEnergy -{ -}; - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -class LaxFridrichsEnergy< Meshes::Grid< 1,MeshReal, Device, MeshIndex >, Real, Index > +class LaxFridrichsEnergyBase { public: - typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType; - typedef typename MeshType::CoordinatesType CoordinatesType; + typedef Real RealType; - typedef Device DeviceType; typedef Index IndexType; + typedef Mesh MeshType; + typedef typename MeshType::DeviceType DeviceType; + typedef typename MeshType::CoordinatesType CoordinatesType; typedef Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocityX; - MeshFunctionType velocityY; - MeshFunctionType pressure; + static const int Dimensions = MeshType::getMeshDimensions(); + typedef Functions::VectorField< Dimensions, MeshFunctionType > VelocityFieldType; + typedef SharedPointer< MeshFunctionType > MeshFunctionPointer; + typedef SharedPointer< VelocityFieldType > VelocityFieldPointer; + + static String getType() + { + return String( "LaxFridrichsEnergy< " ) + + MeshType::getType() + ", " + + TNL::getType< Real >() + ", " + + TNL::getType< Index >() + " >"; + } void setTau(const Real& tau) { this->tau = tau; }; - - void setVelocityX(const MeshFunctionType& velocityX) + + void setVelocity( const VelocityFieldPointer& velocity ) { - this->velocityX = velocityX; + this->velocity = velocity; }; - - void setVelocityY(const MeshFunctionType& velocityY) - { - this->velocityY = velocityY; - }; - - void setPressure(const MeshFunctionType& pressure) + + void setPressure( const MeshFunctionPointer& pressure ) { this->pressure = pressure; }; + protected: + + RealType tau; + + VelocityFieldPointer velocity; + + MeshFunctionPointer pressure; +}; + +template< typename Mesh, + typename Real = typename Mesh::RealType, + typename Index = typename Mesh::IndexType > +class LaxFridrichsEnergy +{ +}; + +template< typename MeshReal, + typename Device, + typename MeshIndex, + typename Real, + typename Index > +class LaxFridrichsEnergy< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index > + : public LaxFridrichsEnergyBase< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index > +{ + public: + + typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType; + typedef LaxFridrichsEnergyBase< MeshType, Real, Index > BaseType; + + using typename BaseType::RealType; + using typename BaseType::IndexType; + using typename BaseType::DeviceType; + using typename BaseType::CoordinatesType; + using typename BaseType::MeshFunctionType; + using typename BaseType::MeshFunctionPointer; + using typename BaseType::VelocityFieldType; + using typename BaseType::VelocityFieldPointer; + using BaseType::Dimensions; + template< typename MeshFunction, typename MeshEntity > __cuda_callable__ - Real operator()( const MeshFunction& u, + Real operator()( const MeshFunction& e, const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > + const RealType& time = 0.0 ) const + { + static_assert( MeshEntity::entityDimensions == 1, "Wrong mesh entity dimensions." ); + static_assert( MeshFunction::getEntitiesDimensions() == 1, "Wrong preimage function" ); + const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); + + const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); + const IndexType& center = entity.getIndex(); + const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); + const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); + const RealType& pressure_west = this->pressure.template getData< DeviceType >()[ west ]; + const RealType& pressure_east = this->pressure.template getData< DeviceType >()[ east ]; + const RealType& velocity_x_east = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ]; + const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ]; + return 0.5 * ( this->tau * ( e[ west ] - 2.0 * e[ center ] + e[ east ] ) + - ( ( e[ west ] + pressure_west ) * velocity_x_west + - ( e[ east ] + pressure_east ) * velocity_x_east ) * hxInverse ); + + } + + /*template< typename MeshEntity > __cuda_callable__ Index getLinearSystemRowLength( const MeshType& mesh, const IndexType& index, @@ -76,7 +134,7 @@ class LaxFridrichsEnergy< Meshes::Grid< 1,MeshReal, Device, MeshIndex >, Real, I const MeshEntity& entity, const MeshFunctionType& u, Vector& b, - MatrixRow& matrixRow ) const; + MatrixRow& matrixRow ) const;*/ }; template< typename MeshReal, @@ -84,50 +142,58 @@ template< typename MeshReal, typename MeshIndex, typename Real, typename Index > -class LaxFridrichsEnergy< Meshes::Grid< 2,MeshReal, Device, MeshIndex >, Real, Index > +class LaxFridrichsEnergy< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index > + : public LaxFridrichsEnergyBase< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocityX; - MeshFunctionType velocityY; - MeshFunctionType pressure; - - void setTau(const Real& tau) - { - this->tau = tau; - }; - - void setVelocityX(const MeshFunctionType& velocityX) - { - this->velocityX = velocityX; - }; - - void setVelocityY(const MeshFunctionType& velocityY) - { - this->velocityY = velocityY; - }; - - void setPressure(const MeshFunctionType& pressure) - { - this->pressure = pressure; - }; + typedef LaxFridrichsEnergyBase< MeshType, Real, Index > BaseType; + + using typename BaseType::RealType; + using typename BaseType::IndexType; + using typename BaseType::DeviceType; + using typename BaseType::CoordinatesType; + using typename BaseType::MeshFunctionType; + using typename BaseType::MeshFunctionPointer; + using typename BaseType::VelocityFieldType; + using typename BaseType::VelocityFieldPointer; + using BaseType::Dimensions; + template< typename MeshFunction, typename MeshEntity > __cuda_callable__ - Real operator()( const MeshFunction& u, + Real operator()( const MeshFunction& e, const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > + const RealType& time = 0.0 ) const + { + static_assert( MeshEntity::entityDimensions == 2, "Wrong mesh entity dimensions." ); + static_assert( MeshFunction::getEntitiesDimensions() == 2, "Wrong preimage function" ); + const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); + + const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1, 0 >(); + const RealType& hyInverse = entity.getMesh().template getSpaceStepsProducts< 0, -1 >(); + const IndexType& center = entity.getIndex(); + const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); + const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); + const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); + const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); + const RealType& pressure_west = this->pressure.template getData< DeviceType >()[ west ]; + const RealType& pressure_east = this->pressure.template getData< DeviceType >()[ east ]; + const RealType& pressure_north = this->pressure.template getData< DeviceType >()[ north ]; + const RealType& pressure_south = this->pressure.template getData< DeviceType >()[ south ]; + const RealType& velocity_x_east = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ]; + const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ]; + const RealType& velocity_y_north = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ north ]; + const RealType& velocity_y_south = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ]; + + return 0.5 * ( this->tau * ( e[ west ] + e[ east ] + e[ south ] + e[ north ] - 4.0 * e[ center ] ) + - ( ( ( e[ west ] + pressure_west ) * velocity_x_west ) + -( ( e[ east ] + pressure_east ) * velocity_x_east ) ) * hxInverse + - ( ( ( e[ north ] + pressure_north ) * velocity_y_north ) + -( ( e[ south ] + pressure_south ) * velocity_y_south ) ) * hyInverse ); + } + + /*template< typename MeshEntity > __cuda_callable__ Index getLinearSystemRowLength( const MeshType& mesh, const IndexType& index, @@ -142,7 +208,7 @@ class LaxFridrichsEnergy< Meshes::Grid< 2,MeshReal, Device, MeshIndex >, Real, I const MeshEntity& entity, const MeshFunctionType& u, Vector& b, - MatrixRow& matrixRow ) const; + MatrixRow& matrixRow ) const;*/ }; template< typename MeshReal, @@ -150,50 +216,68 @@ template< typename MeshReal, typename MeshIndex, typename Real, typename Index > -class LaxFridrichsEnergy< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, Index > +class LaxFridrichsEnergy< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index > + : public LaxFridrichsEnergyBase< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocityX; - MeshFunctionType velocityY; - MeshFunctionType pressure; - - void setTau(const Real& tau) - { - this->tau = tau; - }; - - void setVelocityX(const MeshFunctionType& velocityX) - { - this->velocityX = velocityX; - }; - - void setVelocityY(const MeshFunctionType& velocityY) - { - this->velocityY = velocityY; - }; - - void setPressure(const MeshFunctionType& pressure) - { - this->pressure = pressure; - }; + typedef LaxFridrichsEnergyBase< MeshType, Real, Index > BaseType; + + using typename BaseType::RealType; + using typename BaseType::IndexType; + using typename BaseType::DeviceType; + using typename BaseType::CoordinatesType; + using typename BaseType::MeshFunctionType; + using typename BaseType::MeshFunctionPointer; + using typename BaseType::VelocityFieldType; + using typename BaseType::VelocityFieldPointer; + using BaseType::Dimensions; template< typename MeshFunction, typename MeshEntity > __cuda_callable__ - Real operator()( const MeshFunction& u, + Real operator()( const MeshFunction& e, const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > + const RealType& time = 0.0 ) const + { + static_assert( MeshEntity::entityDimensions == 3, "Wrong mesh entity dimensions." ); + static_assert( MeshFunction::getEntitiesDimensions() == 3, "Wrong preimage function" ); + const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); + + const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1, 0, 0 >(); + const RealType& hyInverse = entity.getMesh().template getSpaceStepsProducts< 0, -1, 0 >(); + const RealType& hzInverse = entity.getMesh().template getSpaceStepsProducts< 0, 0, -1 >(); + const IndexType& center = entity.getIndex(); + const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); + const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); + const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); + const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); + const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); + const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); + + const RealType& pressure_west = this->pressure.template getData< DeviceType >()[ west ]; + const RealType& pressure_east = this->pressure.template getData< DeviceType >()[ east ]; + const RealType& pressure_north = this->pressure.template getData< DeviceType >()[ north ]; + const RealType& pressure_south = this->pressure.template getData< DeviceType >()[ south ]; + const RealType& pressure_up = this->pressure.template getData< DeviceType >()[ up ]; + const RealType& pressure_down = this->pressure.template getData< DeviceType >()[ down ]; + + const RealType& velocity_x_east = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ]; + const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ]; + const RealType& velocity_y_north = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ north ]; + const RealType& velocity_y_south = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ]; + const RealType& velocity_z_up = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ up ]; + const RealType& velocity_z_down = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ down ]; + + return 0.5 * ( this->tau * ( e[ west ] + e[ east ] + e[ south ] + e[ north ] + e[ up ] + e[ down ] - 6.0 * e[ center ] ) + - ( ( ( e[ west ] + pressure_west ) * velocity_x_west ) + -( ( e[ east ] + pressure_east ) * velocity_x_east ) ) * hxInverse + - ( ( ( e[ north ] + pressure_north ) * velocity_y_north ) + -( ( e[ south ] + pressure_south ) * velocity_y_south ) ) * hyInverse + - ( ( ( e[ up ] + pressure_up ) * velocity_z_up ) + -( ( e[ down ] + pressure_down ) * velocity_z_down ) ) * hzInverse ); + } + + /*template< typename MeshEntity > __cuda_callable__ Index getLinearSystemRowLength( const MeshType& mesh, const IndexType& index, @@ -208,12 +292,7 @@ class LaxFridrichsEnergy< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, I const MeshEntity& entity, const MeshFunctionType& u, Vector& b, - MatrixRow& matrixRow ) const; + MatrixRow& matrixRow ) const;*/ }; } //namespace TNL - - -#include "LaxFridrichsEnergy_impl.h" - -#endif /* LaxFridrichsEnergy_H */ diff --git a/examples/inviscid-flow/LaxFridrichsEnergy_impl.h b/examples/inviscid-flow/LaxFridrichsEnergy_impl.h deleted file mode 100644 index e4f9b2e176..0000000000 --- a/examples/inviscid-flow/LaxFridrichsEnergy_impl.h +++ /dev/null @@ -1,350 +0,0 @@ -#ifndef LaxFridrichsEnergy_IMPL_H -#define LaxFridrichsEnergy_IMPL_H - -namespace TNL { - -/**** - * 1D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsEnergy< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsEnergy< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsEnergy< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 1, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 1, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& hxSquareInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsEnergy< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsEnergy< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - matrixRow.setElement( 0, west, - lambdaX ); - matrixRow.setElement( 1, center, 2.0 * lambdaX ); - matrixRow.setElement( 2, east, - lambdaX ); -} - -/**** - * 2D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsEnergy< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsEnergy< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsEnergy< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 2, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 2, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - //energy - const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1, 0 >(); - const RealType& hyInverse = entity.getMesh().template getSpaceStepsProducts< 0, -1 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); - return (0.5 * this->tau) * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] - 4.0 * u[ center ] ) - - 0.5 * hxInverse * ((( u[ west ] + this->pressure[ west ] ) * this->velocityX[ west ] ) - -(( u[ east ] + this->pressure[ east ] ) * this->velocityX[ east ] )) - - 0.5 * hyInverse * ((( u[ north ] + this->pressure[ north ] ) * this->velocityY[ north ] ) - -(( u[ south ] + this->pressure[ south ] ) * this->velocityY[ south ] )); -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsEnergy< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsEnergy< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2, 0 >(); - const RealType& lambdaY = tau * entity.getMesh().template getSpaceStepsProducts< 0, -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); - matrixRow.setElement( 0, south, -lambdaY ); - matrixRow.setElement( 1, west, -lambdaX ); - matrixRow.setElement( 2, center, 2.0 * ( lambdaX + lambdaY ) ); - matrixRow.setElement( 3, east, -lambdaX ); - matrixRow.setElement( 4, north, -lambdaY ); -} - -/**** - * 3D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsEnergy< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsEnergy< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsEnergy< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 3, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 3, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse + - ( u[ south ] - 2.0 * u[ center ] + u[ north ] ) * hySquareInverse + - ( u[ up ] - 2.0 * u[ center ] + u[ down ] ) * hzSquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsEnergy< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsEnergy< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - matrixRow.setElement( 0, down, -lambdaZ ); - matrixRow.setElement( 1, south, -lambdaY ); - matrixRow.setElement( 2, west, -lambdaX ); - matrixRow.setElement( 3, center, 2.0 * ( lambdaX + lambdaY + lambdaZ ) ); - matrixRow.setElement( 4, east, -lambdaX ); - matrixRow.setElement( 5, north, -lambdaY ); - matrixRow.setElement( 6, up, -lambdaZ ); -} - -} //namespace TNL - -#endif /* LaxFridrichsEnergyIMPL_H */ - diff --git a/examples/inviscid-flow/LaxFridrichsMomentumBase.h b/examples/inviscid-flow/LaxFridrichsMomentumBase.h new file mode 100644 index 0000000000..331796c099 --- /dev/null +++ b/examples/inviscid-flow/LaxFridrichsMomentumBase.h @@ -0,0 +1,59 @@ +/*************************************************************************** + LaxFridrichsMomentumBase.h - description + ------------------- + begin : Feb 17, 2017 + copyright : (C) 2017 by Tomas Oberhuber + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/* See Copyright Notice in tnl/Copyright */ + + +#pragma once + +namespace TNL { + +template< typename Mesh, + typename Real = typename Mesh::RealType, + typename Index = typename Mesh::IndexType > +class LaxFridrichsMomentumBase +{ + public: + + typedef Real RealType; + typedef Index IndexType; + typedef Mesh MeshType; + typedef typename MeshType::DeviceType DeviceType; + typedef typename MeshType::CoordinatesType CoordinatesType; + typedef Functions::MeshFunction< MeshType > MeshFunctionType; + static const int Dimensions = MeshType::getMeshDimensions(); + typedef Functions::VectorField< Dimensions, MeshFunctionType > VelocityFieldType; + typedef SharedPointer< MeshFunctionType > MeshFunctionPointer; + typedef SharedPointer< VelocityFieldType > VelocityFieldPointer; + + + void setTau(const Real& tau) + { + this->tau = tau; + }; + + void setVelocity( const VelocityFieldPointer& velocity ) + { + this->velocity = velocity; + }; + + void setPressure( const MeshFunctionPointer& pressure ) + { + this->pressure = pressure; + }; + + protected: + + RealType tau; + + VelocityFieldPointer velocity; + + MeshFunctionPointer pressure; +}; + +} //namespace TNL diff --git a/examples/inviscid-flow/LaxFridrichsMomentumX.h b/examples/inviscid-flow/LaxFridrichsMomentumX.h index b82758fc7c..df63e366e0 100644 --- a/examples/inviscid-flow/LaxFridrichsMomentumX.h +++ b/examples/inviscid-flow/LaxFridrichsMomentumX.h @@ -1,8 +1,19 @@ -#ifndef LaxFridrichsMomentumX_H -#define LaxFridrichsMomentumX_H +/*************************************************************************** + LaxFridrichsMomentumX.h - description + ------------------- + begin : Feb 18, 2017 + copyright : (C) 2017 by Tomas Oberhuber + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/* See Copyright Notice in tnl/Copyright */ + + +#pragma once #include <TNL/Containers/Vector.h> #include <TNL/Meshes/Grid.h> +#include "LaxFridrichsMomentumBase.h" namespace TNL { @@ -18,51 +29,58 @@ template< typename MeshReal, typename MeshIndex, typename Real, typename Index > -class LaxFridrichsMomentumX< Meshes::Grid< 1,MeshReal, Device, MeshIndex >, Real, Index > +class LaxFridrichsMomentumX< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index > + : public LaxFridrichsMomentumBase< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocityX; - MeshFunctionType velocityY; - MeshFunctionType pressure; - - void setTau(const Real& tau) - { - this->tau = tau; - }; - - void setVelocityX(const MeshFunctionType& velocityX) - { - this->velocityX = velocityX; - }; - void setVelocityY(const MeshFunctionType& velocityY) - { - this->velocityY = velocityY; - }; - - void setPressure(const MeshFunctionType& pressure) + typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType; + typedef LaxFridrichsMomentumBase< MeshType, Real, Index > BaseType; + + using typename BaseType::RealType; + using typename BaseType::IndexType; + using typename BaseType::DeviceType; + using typename BaseType::CoordinatesType; + using typename BaseType::MeshFunctionType; + using typename BaseType::MeshFunctionPointer; + using typename BaseType::VelocityFieldType; + using typename BaseType::VelocityFieldPointer; + using BaseType::Dimensions; + + static String getType() { - this->pressure = pressure; - }; - + return String( "LaxFridrichsMomentumX< " ) + + MeshType::getType() + ", " + + TNL::getType< Real >() + ", " + + TNL::getType< Index >() + " >"; + } + template< typename MeshFunction, typename MeshEntity > __cuda_callable__ Real operator()( const MeshFunction& u, const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > + const RealType& time = 0.0 ) const + { + static_assert( MeshEntity::entityDimensions == 1, "Wrong mesh entity dimensions." ); + static_assert( MeshFunction::getEntitiesDimensions() == 1, "Wrong preimage function" ); + const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); + + const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); + const IndexType& center = entity.getIndex(); + const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); + const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); + const RealType& pressure_west = this->pressure.template getData< DeviceType >()[ west ]; + const RealType& pressure_east = this->pressure.template getData< DeviceType >()[ east ]; + const RealType& velocity_x_east = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ]; + const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ]; + + return 0.5 * ( this->tau * ( u[ west ] + u[ east ] - 2.0 * u[ center ] ) + - ( ( u[ west ] * velocity_x_west + pressure_west ) + - ( u[ east ] * velocity_x_east + pressure_east ) ) * hxInverse ); + } + + /*template< typename MeshEntity > __cuda_callable__ Index getLinearSystemRowLength( const MeshType& mesh, const IndexType& index, @@ -77,7 +95,7 @@ class LaxFridrichsMomentumX< Meshes::Grid< 1,MeshReal, Device, MeshIndex >, Real const MeshEntity& entity, const MeshFunctionType& u, Vector& b, - MatrixRow& matrixRow ) const; + MatrixRow& matrixRow ) const;*/ }; template< typename MeshReal, @@ -85,50 +103,65 @@ template< typename MeshReal, typename MeshIndex, typename Real, typename Index > -class LaxFridrichsMomentumX< Meshes::Grid< 2,MeshReal, Device, MeshIndex >, Real, Index > +class LaxFridrichsMomentumX< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index > + : public LaxFridrichsMomentumBase< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocityX; - MeshFunctionType velocityY; - MeshFunctionType pressure; - - void setTau(const Real& tau) - { - this->tau = tau; - }; - - void setVelocityX(const MeshFunctionType& velocityX) + typedef LaxFridrichsMomentumBase< MeshType, Real, Index > BaseType; + + using typename BaseType::RealType; + using typename BaseType::IndexType; + using typename BaseType::DeviceType; + using typename BaseType::CoordinatesType; + using typename BaseType::MeshFunctionType; + using typename BaseType::MeshFunctionPointer; + using typename BaseType::VelocityFieldType; + using typename BaseType::VelocityFieldPointer; + using BaseType::Dimensions; + + static String getType() { - this->velocityX = velocityX; - }; - - void setVelocityY(const MeshFunctionType& velocityY) - { - this->velocityY = velocityY; - }; - - void setPressure(const MeshFunctionType& pressure) - { - this->pressure = pressure; - }; + return String( "LaxFridrichsMomentumX< " ) + + MeshType::getType() + ", " + + TNL::getType< Real >() + ", " + + TNL::getType< Index >() + " >"; + } template< typename MeshFunction, typename MeshEntity > __cuda_callable__ Real operator()( const MeshFunction& u, const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > + const RealType& time = 0.0 ) const + { + static_assert( MeshEntity::entityDimensions == 2, "Wrong mesh entity dimensions." ); + static_assert( MeshFunction::getEntitiesDimensions() == 2, "Wrong preimage function" ); + const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); + + + const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1, 0 >(); + const RealType& hyInverse = entity.getMesh().template getSpaceStepsProducts< 0, -1 >(); + const IndexType& center = entity.getIndex(); + const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); + const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); + const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); + const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); + + const RealType& pressure_west = this->pressure.template getData< DeviceType >()[ west ]; + const RealType& pressure_east = this->pressure.template getData< DeviceType >()[ east ]; + const RealType& velocity_x_east = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ]; + const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ]; + const RealType& velocity_y_north = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ north ]; + const RealType& velocity_y_south = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ]; + + return 0.5 * ( this->tau * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] - 4.0 * u[ center ] ) + - ( ( u[ west ] * velocity_x_west + pressure_west ) + - ( u[ east ] * velocity_x_east + pressure_east ) )* hxInverse + - ( ( u[ north ] * velocity_y_north ) + - ( u[ south ] * velocity_y_south ) )* hyInverse ); + } + + /*template< typename MeshEntity > __cuda_callable__ Index getLinearSystemRowLength( const MeshType& mesh, const IndexType& index, @@ -143,7 +176,7 @@ class LaxFridrichsMomentumX< Meshes::Grid< 2,MeshReal, Device, MeshIndex >, Real const MeshEntity& entity, const MeshFunctionType& u, Vector& b, - MatrixRow& matrixRow ) const; + MatrixRow& matrixRow ) const;*/ }; template< typename MeshReal, @@ -152,49 +185,74 @@ template< typename MeshReal, typename Real, typename Index > class LaxFridrichsMomentumX< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, Index > + : public LaxFridrichsMomentumBase< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocityX; - MeshFunctionType velocityY; - MeshFunctionType pressure; - - void setTau(const Real& tau) - { - this->tau = tau; - }; - - void setVelocityX(const MeshFunctionType& velocityX) + typedef LaxFridrichsMomentumBase< MeshType, Real, Index > BaseType; + + using typename BaseType::RealType; + using typename BaseType::IndexType; + using typename BaseType::DeviceType; + using typename BaseType::CoordinatesType; + using typename BaseType::MeshFunctionType; + using typename BaseType::MeshFunctionPointer; + using typename BaseType::VelocityFieldType; + using typename BaseType::VelocityFieldPointer; + using BaseType::Dimensions; + + static String getType() { - this->velocityX = velocityX; - }; - - void setVelocityY(const MeshFunctionType& velocityY) - { - this->velocityY = velocityY; - }; - - void setPressure(const MeshFunctionType& pressure) - { - this->pressure = pressure; - }; + return String( "LaxFridrichsMomentumX< " ) + + MeshType::getType() + ", " + + TNL::getType< Real >() + ", " + + TNL::getType< Index >() + " >"; + } template< typename MeshFunction, typename MeshEntity > __cuda_callable__ Real operator()( const MeshFunction& u, const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > + const RealType& time = 0.0 ) const + { + static_assert( MeshEntity::entityDimensions == 3, "Wrong mesh entity dimensions." ); + static_assert( MeshFunction::getEntitiesDimensions() == 3, "Wrong preimage function" ); + const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); + + const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1, 0, 0 >(); + const RealType& hyInverse = entity.getMesh().template getSpaceStepsProducts< 0, -1, 0 >(); + const RealType& hzInverse = entity.getMesh().template getSpaceStepsProducts< 0, 0, -1 >(); + const IndexType& center = entity.getIndex(); + const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); + const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); + const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); + const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); + const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); + const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); + + const RealType& pressure_west = this->pressure.template getData< DeviceType >()[ west ]; + const RealType& pressure_east = this->pressure.template getData< DeviceType >()[ east ]; + //const RealType& pressure_north = this->pressure.template getData< DeviceType >()[ north ]; + //const RealType& pressure_south = this->pressure.template getData< DeviceType >()[ south ]; + //const RealType& pressure_up = this->pressure.template getData< DeviceType >()[ up ]; + //const RealType& pressure_down = this->pressure.template getData< DeviceType >()[ down ]; + + const RealType& velocity_x_east = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ]; + const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ]; + const RealType& velocity_y_north = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ north ]; + const RealType& velocity_y_south = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ]; + const RealType& velocity_z_up = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ up ]; + const RealType& velocity_z_down = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ down ]; + return 0.5 * ( this->tau * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] + u[ up ] + u[ down ] - 6.0 * u[ center ] ) + - ( ( u[ west ] * velocity_x_west + pressure_west ) + - ( u[ east ] * velocity_x_east + pressure_east ) )* hxInverse + - ( ( u[ north ] * velocity_y_north ) + - ( u[ south ] * velocity_y_south ) )* hyInverse + - ( ( u[ up ] * velocity_z_up ) + - ( u[ down ] * velocity_z_down ) )* hzInverse ); + } + + /*template< typename MeshEntity > __cuda_callable__ Index getLinearSystemRowLength( const MeshType& mesh, const IndexType& index, @@ -209,12 +267,9 @@ class LaxFridrichsMomentumX< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real const MeshEntity& entity, const MeshFunctionType& u, Vector& b, - MatrixRow& matrixRow ) const; + MatrixRow& matrixRow ) const;*/ }; } // namespace TNL -#include "LaxFridrichsMomentumX_impl.h" - -#endif /* LaxFridrichsMomentumX_H */ diff --git a/examples/inviscid-flow/LaxFridrichsMomentumX_impl.h b/examples/inviscid-flow/LaxFridrichsMomentumX_impl.h deleted file mode 100644 index bffb1fe02b..0000000000 --- a/examples/inviscid-flow/LaxFridrichsMomentumX_impl.h +++ /dev/null @@ -1,351 +0,0 @@ -#ifndef LaxFridrichsMomentumX_IMPL_H -#define LaxFridrichsMomentumX_IMPL_H - -namespace TNL { - -/**** - * 1D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsMomentumX< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsMomentumX< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsMomentumX< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 1, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 1, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& hxSquareInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsMomentumX< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsMomentumX< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - matrixRow.setElement( 0, west, - lambdaX ); - matrixRow.setElement( 1, center, 2.0 * lambdaX ); - matrixRow.setElement( 2, east, - lambdaX ); -} - -/**** - * 2D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsMomentumX< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsMomentumX< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsMomentumX< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 2, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 2, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - - //rhoVelX - const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1, 0 >(); - const RealType& hyInverse = entity.getMesh().template getSpaceStepsProducts< 0, -1 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); - return (0.5 * this->tau) * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] - 4.0 * u[ center ] ) - - 0.5 * hxInverse * (( u[ west ] * this->velocityX[ west ] + this->pressure[ west ] ) - -( u[ east ] * this->velocityX[ east ] + this->pressure[ east ] )) - - 0.5 * hyInverse * (( u[ north ] * this->velocityY[ north ] ) - -( u[ south ] * this->velocityY[ south ] )); -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsMomentumX< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsMomentumX< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2, 0 >(); - const RealType& lambdaY = tau * entity.getMesh().template getSpaceStepsProducts< 0, -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); - matrixRow.setElement( 0, south, -lambdaY ); - matrixRow.setElement( 1, west, -lambdaX ); - matrixRow.setElement( 2, center, 2.0 * ( lambdaX + lambdaY ) ); - matrixRow.setElement( 3, east, -lambdaX ); - matrixRow.setElement( 4, north, -lambdaY ); -} - -/**** - * 3D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsMomentumX< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsMomentumX< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsMomentumX< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 3, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 3, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse + - ( u[ south ] - 2.0 * u[ center ] + u[ north ] ) * hySquareInverse + - ( u[ up ] - 2.0 * u[ center ] + u[ down ] ) * hzSquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsMomentumX< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsMomentumX< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - matrixRow.setElement( 0, down, -lambdaZ ); - matrixRow.setElement( 1, south, -lambdaY ); - matrixRow.setElement( 2, west, -lambdaX ); - matrixRow.setElement( 3, center, 2.0 * ( lambdaX + lambdaY + lambdaZ ) ); - matrixRow.setElement( 4, east, -lambdaX ); - matrixRow.setElement( 5, north, -lambdaY ); - matrixRow.setElement( 6, up, -lambdaZ ); -} - -} //namespace TNL - -#endif /* LaxFridrichsMomentumXIMPL_H */ - diff --git a/examples/inviscid-flow/LaxFridrichsMomentumY.h b/examples/inviscid-flow/LaxFridrichsMomentumY.h index 309051829f..5b67af0841 100644 --- a/examples/inviscid-flow/LaxFridrichsMomentumY.h +++ b/examples/inviscid-flow/LaxFridrichsMomentumY.h @@ -1,8 +1,19 @@ -#ifndef LaxFridrichsMomentumY_H -#define LaxFridrichsMomentumY_H +/*************************************************************************** + LaxFridrichsMomentumY.h - description + ------------------- + begin : Feb 18, 2017 + copyright : (C) 2017 by Tomas Oberhuber + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/* See Copyright Notice in tnl/Copyright */ + + +#pragma once #include <TNL/Containers/Vector.h> #include <TNL/Meshes/Grid.h> +#include "LaxFridrichsMomentumBase.h" namespace TNL { @@ -18,50 +29,47 @@ template< typename MeshReal, typename MeshIndex, typename Real, typename Index > -class LaxFridrichsMomentumY< Meshes::Grid< 1,MeshReal, Device, MeshIndex >, Real, Index > +class LaxFridrichsMomentumY< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index > + : public LaxFridrichsMomentumBase< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocityX; - MeshFunctionType velocityY; - MeshFunctionType pressure; - - void setTau(const Real& tau) - { - this->tau = tau; - }; - - void setVelocityX(const MeshFunctionType& velocityX) - { - this->velocityX = velocityX; - }; - void setVelocityY(const MeshFunctionType& velocityY) - { - this->velocityY = velocityY; - }; - - void setPressure(const MeshFunctionType& pressure) + typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType; + typedef LaxFridrichsMomentumBase< MeshType, Real, Index > BaseType; + + using typename BaseType::RealType; + using typename BaseType::IndexType; + using typename BaseType::DeviceType; + using typename BaseType::CoordinatesType; + using typename BaseType::MeshFunctionType; + using typename BaseType::MeshFunctionPointer; + using typename BaseType::VelocityFieldType; + using typename BaseType::VelocityFieldPointer; + using BaseType::Dimensions; + + static String getType() { - this->pressure = pressure; - }; + return String( "LaxFridrichsMomentumY< " ) + + MeshType::getType() + ", " + + TNL::getType< Real >() + ", " + + TNL::getType< Index >() + " >"; + } + template< typename MeshFunction, typename MeshEntity > __cuda_callable__ Real operator()( const MeshFunction& u, const MeshEntity& entity, - const RealType& time = 0.0 ) const; + const RealType& time = 0.0 ) const + { + static_assert( MeshEntity::entityDimensions == 1, "Wrong mesh entity dimensions." ); + static_assert( MeshFunction::getEntitiesDimensions() == 1, "Wrong preimage function" ); + const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); + + return 0.0; + } - template< typename MeshEntity > + /*template< typename MeshEntity > __cuda_callable__ Index getLinearSystemRowLength( const MeshType& mesh, const IndexType& index, @@ -76,7 +84,7 @@ class LaxFridrichsMomentumY< Meshes::Grid< 1,MeshReal, Device, MeshIndex >, Real const MeshEntity& entity, const MeshFunctionType& u, Vector& b, - MatrixRow& matrixRow ) const; + MatrixRow& matrixRow ) const;*/ }; template< typename MeshReal, @@ -84,50 +92,65 @@ template< typename MeshReal, typename MeshIndex, typename Real, typename Index > -class LaxFridrichsMomentumY< Meshes::Grid< 2,MeshReal, Device, MeshIndex >, Real, Index > +class LaxFridrichsMomentumY< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index > + : public LaxFridrichsMomentumBase< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocityX; - MeshFunctionType velocityY; - MeshFunctionType pressure; - - void setTau(const Real& tau) - { - this->tau = tau; - }; - - void setVelocityX(const MeshFunctionType& velocityX) - { - this->velocityX = velocityX; - }; - - void setVelocityY(const MeshFunctionType& velocityY) - { - this->velocityY = velocityY; - }; - - void setPressure(const MeshFunctionType& pressure) + typedef LaxFridrichsMomentumBase< MeshType, Real, Index > BaseType; + + using typename BaseType::RealType; + using typename BaseType::IndexType; + using typename BaseType::DeviceType; + using typename BaseType::CoordinatesType; + using typename BaseType::MeshFunctionType; + using typename BaseType::MeshFunctionPointer; + using typename BaseType::VelocityFieldType; + using typename BaseType::VelocityFieldPointer; + using BaseType::Dimensions; + + static String getType() { - this->pressure = pressure; - }; + return String( "LaxFridrichsMomentumY< " ) + + MeshType::getType() + ", " + + TNL::getType< Real >() + ", " + + TNL::getType< Index >() + " >"; + } template< typename MeshFunction, typename MeshEntity > __cuda_callable__ Real operator()( const MeshFunction& u, const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > + const RealType& time = 0.0 ) const + { + static_assert( MeshEntity::entityDimensions == 2, "Wrong mesh entity dimensions." ); + static_assert( MeshFunction::getEntitiesDimensions() == 2, "Wrong preimage function" ); + const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); + + + const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1, 0 >(); + const RealType& hyInverse = entity.getMesh().template getSpaceStepsProducts< 0, -1 >(); + const IndexType& center = entity.getIndex(); + const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); + const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); + const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); + const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); + + const RealType& pressure_north = this->pressure.template getData< DeviceType >()[ north ]; + const RealType& pressure_south = this->pressure.template getData< DeviceType >()[ south ]; + const RealType& velocity_x_east = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ]; + const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ]; + const RealType& velocity_y_north = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ north ]; + const RealType& velocity_y_south = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ]; + + return 0.5 * ( this->tau * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] - 4.0 * u[ center ] ) + - ( ( u[ west ] * velocity_x_west ) + - ( u[ east ] * velocity_x_east ) )* hxInverse + - ( ( u[ north ] * velocity_y_north + pressure_north ) + - ( u[ south ] * velocity_y_south + pressure_south ) )* hyInverse ); + } + + /*template< typename MeshEntity > __cuda_callable__ Index getLinearSystemRowLength( const MeshType& mesh, const IndexType& index, @@ -142,7 +165,7 @@ class LaxFridrichsMomentumY< Meshes::Grid< 2,MeshReal, Device, MeshIndex >, Real const MeshEntity& entity, const MeshFunctionType& u, Vector& b, - MatrixRow& matrixRow ) const; + MatrixRow& matrixRow ) const;*/ }; template< typename MeshReal, @@ -151,49 +174,69 @@ template< typename MeshReal, typename Real, typename Index > class LaxFridrichsMomentumY< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, Index > + : public LaxFridrichsMomentumBase< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, 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 Functions::MeshFunction< MeshType > MeshFunctionType; - enum { Dimensions = MeshType::getMeshDimensions() }; - - static String getType(); - Real tau; - MeshFunctionType velocityX; - MeshFunctionType velocityY; - MeshFunctionType pressure; - - void setTau(const Real& tau) - { - this->tau = tau; - }; - - void setVelocityX(const MeshFunctionType& velocityX) + typedef LaxFridrichsMomentumBase< MeshType, Real, Index > BaseType; + + using typename BaseType::RealType; + using typename BaseType::IndexType; + using typename BaseType::DeviceType; + using typename BaseType::CoordinatesType; + using typename BaseType::MeshFunctionType; + using typename BaseType::MeshFunctionPointer; + using typename BaseType::VelocityFieldType; + using typename BaseType::VelocityFieldPointer; + using BaseType::Dimensions; + + static String getType() { - this->velocityX = velocityX; - }; - - void setVelocityY(const MeshFunctionType& velocityY) - { - this->velocityY = velocityY; - }; - - void setPressure(const MeshFunctionType& pressure) - { - this->pressure = pressure; - }; + return String( "LaxFridrichsMomentumY< " ) + + MeshType::getType() + ", " + + TNL::getType< Real >() + ", " + + TNL::getType< Index >() + " >"; + } template< typename MeshFunction, typename MeshEntity > __cuda_callable__ Real operator()( const MeshFunction& u, const MeshEntity& entity, - const RealType& time = 0.0 ) const; - - template< typename MeshEntity > + const RealType& time = 0.0 ) const + { + static_assert( MeshEntity::entityDimensions == 3, "Wrong mesh entity dimensions." ); + static_assert( MeshFunction::getEntitiesDimensions() == 3, "Wrong preimage function" ); + const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); + + const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1, 0, 0 >(); + const RealType& hyInverse = entity.getMesh().template getSpaceStepsProducts< 0, -1, 0 >(); + const RealType& hzInverse = entity.getMesh().template getSpaceStepsProducts< 0, 0, -1 >(); + const IndexType& center = entity.getIndex(); + const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); + const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); + const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); + const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); + const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); + const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); + + const RealType& pressure_north = this->pressure.template getData< DeviceType >()[ north ]; + const RealType& pressure_south = this->pressure.template getData< DeviceType >()[ south ]; + const RealType& velocity_x_east = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ]; + const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ]; + const RealType& velocity_y_north = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ north ]; + const RealType& velocity_y_south = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ]; + const RealType& velocity_z_up = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ up ]; + const RealType& velocity_z_down = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ down ]; + return 0.5 * ( this->tau * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] + u[ up ] + u[ down ] - 6.0 * u[ center ] ) + - ( ( u[ west ] * velocity_x_west ) + - ( u[ east ] * velocity_x_east ) )* hxInverse + - ( ( u[ north ] * velocity_y_north + pressure_north ) + - ( u[ south ] * velocity_y_south + pressure_south ) )* hyInverse + - ( ( u[ up ] * velocity_z_up ) + - ( u[ down ] * velocity_z_down ) )* hzInverse ); + } + + /*template< typename MeshEntity > __cuda_callable__ Index getLinearSystemRowLength( const MeshType& mesh, const IndexType& index, @@ -208,11 +251,9 @@ class LaxFridrichsMomentumY< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real const MeshEntity& entity, const MeshFunctionType& u, Vector& b, - MatrixRow& matrixRow ) const; + MatrixRow& matrixRow ) const;*/ }; -} // namespace TNL -#include "LaxFridrichsMomentumY_impl.h" +} // namespace TNL -#endif /* LaxFridrichsMomentumY_H */ diff --git a/examples/inviscid-flow/LaxFridrichsMomentumY_impl.h b/examples/inviscid-flow/LaxFridrichsMomentumY_impl.h deleted file mode 100644 index 0624e19aa8..0000000000 --- a/examples/inviscid-flow/LaxFridrichsMomentumY_impl.h +++ /dev/null @@ -1,351 +0,0 @@ -#ifndef LaxFridrichsMomentumY_IMPL_H -#define LaxFridrichsMomentumY_IMPL_H - -namespace TNL { - -/**** - * 1D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsMomentumY< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsMomentumY< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsMomentumY< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 1, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 1, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& hxSquareInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsMomentumY< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsMomentumY< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); - matrixRow.setElement( 0, west, - lambdaX ); - matrixRow.setElement( 1, center, 2.0 * lambdaX ); - matrixRow.setElement( 2, east, - lambdaX ); -} - -/**** - * 2D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsMomentumY< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsMomentumY< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsMomentumY< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 2, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 2, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - - //rhoVelY - const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1, 0 >(); - const RealType& hyInverse = entity.getMesh().template getSpaceStepsProducts< 0, -1 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); - return (0.5 * this->tau) * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] - 4.0 * u[ center ] ) - - 0.5 * hyInverse * (( u[ west ] * this->velocityX[ west ] ) - -( u[ east ] * this->velocityX[ east ] )) - - 0.5 * hxInverse * (( u[ north ] * this->velocityY[ north ] + this->pressure[ north ] ) - -( u[ south ] * this->velocityY[ south ] + this->pressure[ south ])); -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsMomentumY< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsMomentumY< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); - const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2, 0 >(); - const RealType& lambdaY = tau * entity.getMesh().template getSpaceStepsProducts< 0, -2 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1 >(); - matrixRow.setElement( 0, south, -lambdaY ); - matrixRow.setElement( 1, west, -lambdaX ); - matrixRow.setElement( 2, center, 2.0 * ( lambdaX + lambdaY ) ); - matrixRow.setElement( 3, east, -lambdaX ); - matrixRow.setElement( 4, north, -lambdaY ); -} - -/**** - * 3D problem - */ -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -String -LaxFridrichsMomentumY< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getType() -{ - return String( "LaxFridrichsMomentumY< " ) + - MeshType::getType() + ", " + - TNL::getType< Real >() + ", " + - TNL::getType< Index >() + " >"; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshFunction, typename MeshEntity > -__cuda_callable__ -Real -LaxFridrichsMomentumY< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -operator()( const MeshFunction& u, - const MeshEntity& entity, - const Real& time ) const -{ - /**** - * Implement your explicit form of the differential operator here. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - static_assert( MeshEntity::entityDimensions == 3, "Wrong mesh entity dimensions." ); - static_assert( MeshFunction::getEntitiesDimensions() == 3, "Wrong preimage function" ); - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse + - ( u[ south ] - 2.0 * u[ center ] + u[ north ] ) * hySquareInverse + - ( u[ up ] - 2.0 * u[ center ] + u[ down ] ) * hzSquareInverse; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > -template< typename MeshEntity > -__cuda_callable__ -Index -LaxFridrichsMomentumY< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -getLinearSystemRowLength( const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity ) const -{ - /**** - * Return a number of non-zero elements in a line (associated with given grid element) of - * the linear system. - * The following example is the Laplace operator approximated - * by the Finite difference method. - */ - - return 2*Dimensions + 1; -} - -template< typename MeshReal, - typename Device, - typename MeshIndex, - typename Real, - typename Index > - template< typename MeshEntity, typename Vector, typename MatrixRow > -__cuda_callable__ -void -LaxFridrichsMomentumY< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: -updateLinearSystem( const RealType& time, - const RealType& tau, - const MeshType& mesh, - const IndexType& index, - const MeshEntity& entity, - const MeshFunctionType& u, - Vector& b, - MatrixRow& matrixRow ) const -{ - /**** - * Setup the non-zero elements of the linear system here. - * The following example is the Laplace operator appriximated - * by the Finite difference method. - */ - - const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); - 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 >(); - const IndexType& center = entity.getIndex(); - const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); - const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); - const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); - const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); - const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); - const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); - matrixRow.setElement( 0, down, -lambdaZ ); - matrixRow.setElement( 1, south, -lambdaY ); - matrixRow.setElement( 2, west, -lambdaX ); - matrixRow.setElement( 3, center, 2.0 * ( lambdaX + lambdaY + lambdaZ ) ); - matrixRow.setElement( 4, east, -lambdaX ); - matrixRow.setElement( 5, north, -lambdaY ); - matrixRow.setElement( 6, up, -lambdaZ ); -} - -} //namespace TNL - -#endif /* LaxFridrichsMomentumYIMPL_H */ - diff --git a/examples/inviscid-flow/LaxFridrichsMomentumZ.h b/examples/inviscid-flow/LaxFridrichsMomentumZ.h new file mode 100644 index 0000000000..ba1ea656a7 --- /dev/null +++ b/examples/inviscid-flow/LaxFridrichsMomentumZ.h @@ -0,0 +1,239 @@ +/*************************************************************************** + LaxFridrichsMomentumZ.h - description + ------------------- + begin : Feb 18, 2017 + copyright : (C) 2017 by Tomas Oberhuber + email : tomas.oberhuber@fjfi.cvut.cz + ***************************************************************************/ + +/* See Copyright Notice in tnl/Copyright */ + + +#pragma once + +#include <TNL/Containers/Vector.h> +#include <TNL/Meshes/Grid.h> +#include "LaxFridrichsMomentumBase.h" + +namespace TNL { + +template< typename Mesh, + typename Real = typename Mesh::RealType, + typename Index = typename Mesh::IndexType > +class LaxFridrichsMomentumZ +{ +}; + +template< typename MeshReal, + typename Device, + typename MeshIndex, + typename Real, + typename Index > +class LaxFridrichsMomentumZ< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index > + : public LaxFridrichsMomentumBase< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index > +{ + public: + + typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType; + typedef LaxFridrichsMomentumBase< MeshType, Real, Index > BaseType; + + using typename BaseType::RealType; + using typename BaseType::IndexType; + using typename BaseType::DeviceType; + using typename BaseType::CoordinatesType; + using typename BaseType::MeshFunctionType; + using typename BaseType::MeshFunctionPointer; + using typename BaseType::VelocityFieldType; + using typename BaseType::VelocityFieldPointer; + using BaseType::Dimensions; + + static String getType() + { + return String( "LaxFridrichsMomentumZ< " ) + + MeshType::getType() + ", " + + TNL::getType< Real >() + ", " + + TNL::getType< Index >() + " >"; + } + + + template< typename MeshFunction, typename MeshEntity > + __cuda_callable__ + Real operator()( const MeshFunction& u, + const MeshEntity& entity, + const RealType& time = 0.0 ) const + { + static_assert( MeshEntity::entityDimensions == 1, "Wrong mesh entity dimensions." ); + static_assert( MeshFunction::getEntitiesDimensions() == 1, "Wrong preimage function" ); + const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); + + return 0.0; + } + + /*template< typename MeshEntity > + __cuda_callable__ + Index getLinearSystemRowLength( const MeshType& mesh, + const IndexType& index, + const MeshEntity& entity ) const; + + template< typename MeshEntity, typename Vector, typename MatrixRow > + __cuda_callable__ + void updateLinearSystem( const RealType& time, + const RealType& tau, + const MeshType& mesh, + const IndexType& index, + const MeshEntity& entity, + const MeshFunctionType& u, + Vector& b, + MatrixRow& matrixRow ) const;*/ +}; + +template< typename MeshReal, + typename Device, + typename MeshIndex, + typename Real, + typename Index > +class LaxFridrichsMomentumZ< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index > + : public LaxFridrichsMomentumBase< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index > +{ + public: + typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType; + typedef LaxFridrichsMomentumBase< MeshType, Real, Index > BaseType; + + using typename BaseType::RealType; + using typename BaseType::IndexType; + using typename BaseType::DeviceType; + using typename BaseType::CoordinatesType; + using typename BaseType::MeshFunctionType; + using typename BaseType::MeshFunctionPointer; + using typename BaseType::VelocityFieldType; + using typename BaseType::VelocityFieldPointer; + using BaseType::Dimensions; + + static String getType() + { + return String( "LaxFridrichsMomentumZ< " ) + + MeshType::getType() + ", " + + TNL::getType< Real >() + ", " + + TNL::getType< Index >() + " >"; + } + + template< typename MeshFunction, typename MeshEntity > + __cuda_callable__ + Real operator()( const MeshFunction& u, + const MeshEntity& entity, + const RealType& time = 0.0 ) const + { + static_assert( MeshEntity::entityDimensions == 2, "Wrong mesh entity dimensions." ); + static_assert( MeshFunction::getEntitiesDimensions() == 2, "Wrong preimage function" ); + const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); + + return 0.0; + } + + /*template< typename MeshEntity > + __cuda_callable__ + Index getLinearSystemRowLength( const MeshType& mesh, + const IndexType& index, + const MeshEntity& entity ) const; + + template< typename MeshEntity, typename Vector, typename MatrixRow > + __cuda_callable__ + void updateLinearSystem( const RealType& time, + const RealType& tau, + const MeshType& mesh, + const IndexType& index, + const MeshEntity& entity, + const MeshFunctionType& u, + Vector& b, + MatrixRow& matrixRow ) const;*/ +}; + +template< typename MeshReal, + typename Device, + typename MeshIndex, + typename Real, + typename Index > +class LaxFridrichsMomentumZ< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, Index > + : public LaxFridrichsMomentumBase< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index > +{ + public: + typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType; + typedef LaxFridrichsMomentumBase< MeshType, Real, Index > BaseType; + + using typename BaseType::RealType; + using typename BaseType::IndexType; + using typename BaseType::DeviceType; + using typename BaseType::CoordinatesType; + using typename BaseType::MeshFunctionType; + using typename BaseType::MeshFunctionPointer; + using typename BaseType::VelocityFieldType; + using typename BaseType::VelocityFieldPointer; + using BaseType::Dimensions; + + static String getType() + { + return String( "LaxFridrichsMomentumZ< " ) + + MeshType::getType() + ", " + + TNL::getType< Real >() + ", " + + TNL::getType< Index >() + " >"; + } + + template< typename MeshFunction, typename MeshEntity > + __cuda_callable__ + Real operator()( const MeshFunction& u, + const MeshEntity& entity, + const RealType& time = 0.0 ) const + { + static_assert( MeshEntity::entityDimensions == 3, "Wrong mesh entity dimensions." ); + static_assert( MeshFunction::getEntitiesDimensions() == 3, "Wrong preimage function" ); + const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); + + const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1, 0, 0 >(); + const RealType& hyInverse = entity.getMesh().template getSpaceStepsProducts< 0, -1, 0 >(); + const RealType& hzInverse = entity.getMesh().template getSpaceStepsProducts< 0, 0, -1 >(); + const IndexType& center = entity.getIndex(); + const IndexType& east = neighbourEntities.template getEntityIndex< 1, 0, 0 >(); + const IndexType& west = neighbourEntities.template getEntityIndex< -1, 0, 0 >(); + const IndexType& north = neighbourEntities.template getEntityIndex< 0, 1, 0 >(); + const IndexType& south = neighbourEntities.template getEntityIndex< 0, -1, 0 >(); + const IndexType& up = neighbourEntities.template getEntityIndex< 0, 0, 1 >(); + const IndexType& down = neighbourEntities.template getEntityIndex< 0, 0, -1 >(); + + const RealType& pressure_up = this->pressure.template getData< DeviceType >()[ up ]; + const RealType& pressure_down = this->pressure.template getData< DeviceType >()[ down ]; + const RealType& velocity_x_east = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ]; + const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ]; + const RealType& velocity_y_north = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ north ]; + const RealType& velocity_y_south = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ]; + const RealType& velocity_z_up = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ up ]; + const RealType& velocity_z_down = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ down ]; + return 0.5 * ( this->tau * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] + u[ up ] + u[ down ] - 6.0 * u[ center ] ) + - ( ( u[ west ] * velocity_x_west ) + - ( u[ east ] * velocity_x_east ) )* hxInverse + - ( ( u[ north ] * velocity_y_north ) + - ( u[ south ] * velocity_y_south ) )* hyInverse + - ( ( u[ up ] * velocity_z_up + pressure_up ) + - ( u[ down ] * velocity_z_down + pressure_down ) )* hzInverse ); + } + + /*template< typename MeshEntity > + __cuda_callable__ + Index getLinearSystemRowLength( const MeshType& mesh, + const IndexType& index, + const MeshEntity& entity ) const; + + template< typename MeshEntity, typename Vector, typename MatrixRow > + __cuda_callable__ + void updateLinearSystem( const RealType& time, + const RealType& tau, + const MeshType& mesh, + const IndexType& index, + const MeshEntity& entity, + const MeshFunctionType& u, + Vector& b, + MatrixRow& matrixRow ) const;*/ +}; + + +} // namespace TNL + diff --git a/examples/inviscid-flow/eulerBuildConfigTag.h b/examples/inviscid-flow/eulerBuildConfigTag.h index fef4dfffce..2bb81e8cb3 100644 --- a/examples/inviscid-flow/eulerBuildConfigTag.h +++ b/examples/inviscid-flow/eulerBuildConfigTag.h @@ -21,7 +21,7 @@ template<> struct ConfigTagReal< eulerBuildConfigTag, long double > { enum { ena template<> struct ConfigTagIndex< eulerBuildConfigTag, short int >{ enum { enabled = false }; }; template<> struct ConfigTagIndex< eulerBuildConfigTag, long int >{ enum { enabled = false }; }; -template< int Dimensions > struct ConfigTagDimensions< eulerBuildConfigTag, Dimensions >{ enum { enabled = ( Dimensions == 2 ) }; }; +template< int Dimensions > struct ConfigTagDimensions< eulerBuildConfigTag, Dimensions >{ enum { enabled = true }; }; /**** * Use of Grid is enabled for allowed dimensions and Real, Device and Index types. @@ -43,7 +43,7 @@ template<> struct ConfigTagTimeDiscretisation< eulerBuildConfigTag, ImplicitTime /**** * Only the Runge-Kutta-Merson solver is enabled by default. */ -template<> struct ConfigTagExplicitSolver< eulerBuildConfigTag, ExplicitEulerSolverTag >{ enum { enabled = false }; }; +template<> struct ConfigTagExplicitSolver< eulerBuildConfigTag, ExplicitEulerSolverTag >{ enum { enabled = true }; }; } // namespace Solvers } // namespace TNL diff --git a/examples/inviscid-flow/eulerProblem.h b/examples/inviscid-flow/eulerProblem.h index d12273c894..5e1374d646 100644 --- a/examples/inviscid-flow/eulerProblem.h +++ b/examples/inviscid-flow/eulerProblem.h @@ -60,11 +60,8 @@ class eulerProblem: typedef typename DifferentialOperator::Continuity Continuity; typedef typename DifferentialOperator::MomentumX MomentumX; typedef typename DifferentialOperator::MomentumY MomentumY; + typedef typename DifferentialOperator::MomentumZ MomentumZ; typedef typename DifferentialOperator::Energy Energy; - typedef typename DifferentialOperator::Velocity Velocity; - typedef typename DifferentialOperator::VelocityX VelocityX; - typedef typename DifferentialOperator::VelocityY VelocityY; - typedef typename DifferentialOperator::Pressure Pressure; static String getTypeStatic(); diff --git a/examples/inviscid-flow/eulerProblem_impl.h b/examples/inviscid-flow/eulerProblem_impl.h index 2258a05028..d3bbc1c312 100644 --- a/examples/inviscid-flow/eulerProblem_impl.h +++ b/examples/inviscid-flow/eulerProblem_impl.h @@ -26,10 +26,7 @@ #include "LaxFridrichsEnergy.h" #include "LaxFridrichsMomentumX.h" #include "LaxFridrichsMomentumY.h" -#include "EulerPressureGetter.h" -#include "EulerVelXGetter.h" -#include "EulerVelYGetter.h" -#include "EulerVelGetter.h" +#include "LaxFridrichsMomentumZ.h" namespace TNL { @@ -248,14 +245,13 @@ getExplicitRHS( const RealType& time, SharedPointer< Continuity > continuity; SharedPointer< MomentumX > momentumX; SharedPointer< MomentumY > momentumY; - //SharedPointer< MomentumZ > momentumZ; + SharedPointer< MomentumZ > momentumZ; SharedPointer< Energy > energy; //this->bindDofs( mesh, _u ); //rho continuity->setTau(tau); - continuity->setVelocityX( *( *velocity )[ 0 ] ); - continuity->setVelocityY( *( *velocity )[ 1 ] ); + continuity->setVelocity( this->velocity ); Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, Continuity, BoundaryCondition, RightHandSide > explicitUpdaterContinuity; explicitUpdaterContinuity.template update< typename Mesh::Cell >( time, mesh, @@ -267,9 +263,8 @@ getExplicitRHS( const RealType& time, //rhoVelocityX momentumX->setTau(tau); - momentumX->setVelocityX( *( *velocity )[ 0 ] ); - momentumX->setVelocityY( *( *velocity )[ 1 ] ); - momentumX->setPressure( *pressure ); + momentumX->setVelocity( this->velocity ); + momentumX->setPressure( this->pressure ); Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, MomentumX, BoundaryCondition, RightHandSide > explicitUpdaterMomentumX; explicitUpdaterMomentumX.template update< typename Mesh::Cell >( time, mesh, @@ -280,24 +275,41 @@ getExplicitRHS( const RealType& time, ( *this->conservativeVariablesRHS->getMomentum() )[ 0 ] ); //, fuRhoVelocityX ); //rhoVelocityY - momentumY->setTau(tau); - momentumY->setVelocityX( *( *velocity )[ 0 ] ); - momentumY->setVelocityY( *( *velocity )[ 1 ] ); - momentumY->setPressure( *pressure ); - Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, MomentumY, BoundaryCondition, RightHandSide > explicitUpdaterMomentumY; - explicitUpdaterMomentumY.template update< typename Mesh::Cell >( time, - mesh, - momentumY, - this->boundaryConditionPointer, - this->rightHandSidePointer, - ( *this->conservativeVariables->getMomentum() )[ 1 ], // uRhoVelocityX, - ( *this->conservativeVariablesRHS->getMomentum() )[ 1 ] ); //, fuRhoVelocityX ); + if( Dimensions > 1 ) + { + momentumY->setTau(tau); + momentumY->setVelocity( this->velocity ); + momentumY->setPressure( this->pressure ); + Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, MomentumY, BoundaryCondition, RightHandSide > explicitUpdaterMomentumY; + explicitUpdaterMomentumY.template update< typename Mesh::Cell >( time, + mesh, + momentumY, + this->boundaryConditionPointer, + this->rightHandSidePointer, + ( *this->conservativeVariables->getMomentum() )[ 1 ], // uRhoVelocityX, + ( *this->conservativeVariablesRHS->getMomentum() )[ 1 ] ); //, fuRhoVelocityX ); + } + + if( Dimensions > 2 ) + { + momentumY->setTau(tau); + momentumY->setVelocity( this->velocity ); + momentumY->setPressure( this->pressure ); + Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, MomentumZ, BoundaryCondition, RightHandSide > explicitUpdaterMomentumZ; + explicitUpdaterMomentumZ.template update< typename Mesh::Cell >( time, + mesh, + momentumZ, + this->boundaryConditionPointer, + this->rightHandSidePointer, + ( *this->conservativeVariables->getMomentum() )[ 2 ], // uRhoVelocityX, + ( *this->conservativeVariablesRHS->getMomentum() )[ 2 ] ); //, fuRhoVelocityX ); + } + //energy energy->setTau(tau); - energy->setVelocityX( *( *velocity )[ 0 ] ); - energy->setVelocityY( *( *velocity )[ 1 ] ); - energy->setPressure( *pressure ); + energy->setVelocity( this->velocity ); + energy->setPressure( this->pressure ); Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, Energy, BoundaryCondition, RightHandSide > explicitUpdaterEnergy; explicitUpdaterEnergy.template update< typename Mesh::Cell >( time, mesh, @@ -306,6 +318,11 @@ getExplicitRHS( const RealType& time, this->rightHandSidePointer, this->conservativeVariables->getEnergy(), // uRhoVelocityX, this->conservativeVariablesRHS->getEnergy() ); //, fuRhoVelocityX ); + + /*this->conservativeVariablesRHS->getDensity()->write( "density", "gnuplot" ); + this->conservativeVariablesRHS->getEnergy()->write( "energy", "gnuplot" ); + this->conservativeVariablesRHS->getMomentum()->write( "momentum", "gnuplot", 0.05 ); + getchar();*/ /* BoundaryConditionsSetter< MeshFunctionType, BoundaryCondition > boundaryConditionsSetter; diff --git a/examples/inviscid-flow/run-euler b/examples/inviscid-flow/run-euler index 211dbbf387..7af3493805 100644 --- a/examples/inviscid-flow/run-euler +++ b/examples/inviscid-flow/run-euler @@ -10,22 +10,25 @@ tnl-grid-setup --dimensions 2 \ tnl-init --test-function sin-wave \ --output-file init.tnl -tnl-euler-2d --discontinuity-placement-0 0.3 \ +tnl-euler-2d --initial-condition riemann-problem \ + --discontinuity-placement-0 0.3 \ --discontinuity-placement-1 0.3 \ --discontinuity-placement-2 0.3 \ --left-density 1.0 \ --right-density 1.0 \ - --left-velocity-0 -2.0 \ - --left-velocity-1 -2.0 \ - --left-velocity-2 -2.0 \ - --right-velocity-0 2.0 \ - --right-velocity-1 2.0 \ - --right-velocity-2 2.0 \ + --left-velocity-0 -0.2 \ + --left-velocity-1 -0.2 \ + --left-velocity-2 -0.2 \ + --right-velocity-0 0.2 \ + --right-velocity-1 0.2 \ + --right-velocity-2 0.2 \ --left-pressure 0.4 \ --right-pressure 0.4 \ --time-discretisation explicit \ + --boundary-conditions-type neumann \ --boundary-conditions-constant 0 \ - --discrete-solver merson \ + --discrete-solver euler \ + --time-step 0.0001 \ --snapshot-period 0.01 \ --final-time 1.0 diff --git a/src/TNL/Functions/VectorField.h b/src/TNL/Functions/VectorField.h index c7dadcf7df..5f98b65fcc 100644 --- a/src/TNL/Functions/VectorField.h +++ b/src/TNL/Functions/VectorField.h @@ -183,6 +183,7 @@ class VectorField< Size, MeshFunction< Mesh, MeshEntityDimensions, Real > > VectorType v; for( int i = 0; i < Size; i++ ) v[ i ] = ( *this->vectorField[ i ] )[ index ]; + return v; } template< typename EntityType > @@ -192,6 +193,7 @@ class VectorField< Size, MeshFunction< Mesh, MeshEntityDimensions, Real > > VectorType v; for( int i = 0; i < Size; i++ ) v[ i ] = ( *this->vectorField[ i ] )( meshEntity ); + return v; } bool save( File& file ) const -- GitLab