diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index b0d5d7f7d3f32843b3981c157aee52e5bff13fd5..0f00654e037cfb1ab253fa4276a57fc1f0f25126 100755
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1,3 +1,5 @@
add_subdirectory( heat-equation )
add_subdirectory( navier-stokes )
+add_subdirectory( advection )
+add_subdirectory( inviscid-flow )
#add_subdirectory( mean-curvature-flow )
diff --git a/examples/advection/CMakeLists.txt b/examples/advection/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a7dfa1398d30ff61a3f5540e9b591c589d25b0ea
--- /dev/null
+++ b/examples/advection/CMakeLists.txt
@@ -0,0 +1,20 @@
+set( tnl_heat_equation_SOURCES
+ advection.cpp
+ advection.cu )
+
+IF( BUILD_CUDA )
+ CUDA_ADD_EXECUTABLE( tnl-advection{debugExt} advection.cu)
+ target_link_libraries( tnl-advection${debugExt} tnl${debugExt}-${tnlVersion} ${CUSPARSE_LIBRARY} )
+ELSE( BUILD_CUDA )
+ ADD_EXECUTABLE( tnl-advection${debugExt} advection.cpp)
+ target_link_libraries( tnl-advection${debugExt} tnl${debugExt}-${tnlVersion} )
+ENDIF( BUILD_CUDA )
+
+
+INSTALL( TARGETS tnl-advection${debugExt}
+ RUNTIME DESTINATION bin
+ PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE )
+
+INSTALL( FILES tnl-run-advection
+ ${tnl_heat_equation_SOURCES}
+ DESTINATION share/tnl-${tnlVersion}/examples/advection )
diff --git a/examples/advection/LaxFridrichs.h b/examples/advection/LaxFridrichs.h
new file mode 100644
index 0000000000000000000000000000000000000000..85e445704e63a21eb79a35e21f2103a7e6873bfd
--- /dev/null
+++ b/examples/advection/LaxFridrichs.h
@@ -0,0 +1,218 @@
+#ifndef LaxFridrichs_H
+#define LaxFridrichs_H
+
+#include <core/vectors/tnlVector.h>
+#include <mesh/tnlGrid.h>
+
+template< typename Mesh,
+ typename Real = typename Mesh::RealType,
+ typename Index = typename Mesh::IndexType >
+class LaxFridrichs
+{
+};
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class LaxFridrichs< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index >
+{
+ public:
+ typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef tnlMeshFunction< MeshType > MeshFunctionType;
+ enum { Dimensions = MeshType::getMeshDimensions() };
+ Real tau;
+ double artificalViscosity;
+ double advectionSpeedX;
+ double advectionSpeedY;
+
+ void setAdvectionSpeedY(const double advectionSpeed)
+ {
+ this->advectionSpeedY = advectionSpeed;
+ }
+
+
+ void setAdvectionSpeedX(const double advectionSpeed)
+ {
+ this->advectionSpeedX = advectionSpeed;
+ }
+
+ void setViscosity(const double artificalViscosity)
+ {
+ this->artificalViscosity = artificalViscosity;
+ }
+
+ void setTau(const Real& tau)
+ {
+ this->tau = tau;
+ };
+
+ static tnlString getType();
+
+ template< typename MeshFunction, typename MeshEntity >
+ __cuda_callable__
+ Real operator()( const MeshFunction& u,
+ const MeshEntity& entity,
+ const RealType& time = 0.0 ) const;
+
+ __cuda_callable__
+ template< typename MeshEntity >
+ 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 LaxFridrichs< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index >
+{
+ public:
+ typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef tnlMeshFunction< MeshType > MeshFunctionType;
+ enum { Dimensions = MeshType::getMeshDimensions() };
+ Real tau;
+ double artificalViscosity;
+ double advectionSpeedX;
+ double advectionSpeedY;
+
+ void setAdvectionSpeedY(const double advectionSpeed)
+ {
+ this->advectionSpeedY = advectionSpeed;
+ }
+
+
+ void setAdvectionSpeedX(const double advectionSpeed)
+ {
+ this->advectionSpeedX = advectionSpeed;
+ }
+
+ void setViscosity(const double artificalViscosity)
+ {
+ this->artificalViscosity = artificalViscosity;
+ }
+
+ void setTau(const Real& tau)
+ {
+ this->tau = tau;
+ };
+
+ static tnlString getType();
+
+ template< typename MeshFunction, typename MeshEntity >
+ __cuda_callable__
+ Real operator()( const MeshFunction& u,
+ const MeshEntity& entity,
+ const RealType& time = 0.0 ) const;
+
+ __cuda_callable__
+ template< typename MeshEntity >
+ 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 LaxFridrichs< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, Index >
+{
+ public:
+ typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef tnlMeshFunction< MeshType > MeshFunctionType;
+ enum { Dimensions = MeshType::getMeshDimensions() };
+ Real tau;
+ double artificalViscosity;
+ double advectionSpeedX;
+ double advectionSpeedY;
+
+ void setAdvectionSpeedY(const double advectionSpeed)
+ {
+ this->advectionSpeedY = advectionSpeed;
+ }
+
+
+ void setAdvectionSpeedX(const double advectionSpeed)
+ {
+ this->advectionSpeedX = advectionSpeed;
+ }
+
+ void setViscosity(const double artificalViscosity)
+ {
+ this->artificalViscosity = artificalViscosity;
+ }
+
+ void setTau(const Real& tau)
+ {
+ this->tau = tau;
+ };
+
+ static tnlString getType();
+
+ template< typename MeshFunction, typename MeshEntity >
+ __cuda_callable__
+ Real operator()( const MeshFunction& u,
+ const MeshEntity& entity,
+ const RealType& time = 0.0 ) const;
+
+ __cuda_callable__
+ template< typename MeshEntity >
+ 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;
+};
+
+
+#include "LaxFridrichs_impl.h"
+
+#endif /* LaxFridrichs_H */
diff --git a/examples/advection/LaxFridrichs_impl.h b/examples/advection/LaxFridrichs_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..a82fca1d3613ec5d95b5ba8fca64cfaf9a3cb0d8
--- /dev/null
+++ b/examples/advection/LaxFridrichs_impl.h
@@ -0,0 +1,369 @@
+#ifndef LaxFridrichs_IMPL_H
+#define LaxFridrichs_IMPL_H
+
+/****
+ * 1D problem
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString
+LaxFridrichs< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+ return tnlString( "LaxFridrichs< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename MeshFunction, typename MeshEntity >
+__cuda_callable__
+Real
+LaxFridrichs< tnlGrid< 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& 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 >();
+ double a;
+ a = this->advectionSpeedX;
+ return (0.5 / this->tau ) * this->artificalViscosity *
+ ( u[ west ]
+ - 2.0 * u[ center ]
+ + u[ east ] )
+ - (a * ( u[ east ]
+ - u[west] ) )
+ * hxInverse * 0.5;
+
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename MeshEntity >
+__cuda_callable__
+Index
+LaxFridrichs< tnlGrid< 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
+LaxFridrichs< tnlGrid< 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 >
+tnlString
+LaxFridrichs< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+ return tnlString( "LaxFridrichs< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename MeshFunction, typename MeshEntity >
+__cuda_callable__
+Real
+LaxFridrichs< tnlGrid< 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& 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 >();
+ double a;
+ double b;
+ a = this->advectionSpeedX;
+ b = this->advectionSpeedY;
+ return ( 0.25 / this->tau ) * this->artificalViscosity * (
+ u[ west ]
+ + u[ east ]
+ + u[ south ]
+ + u[ north ]
+ - 4 * u[ center ] )
+ - (a * ( u[ east ]
+ - u[west] ) )
+ * hxInverse * 0.5
+ - (b * ( u[ south ]
+ - u[ south ] ) )
+ * hyInverse * 0.5;
+
+
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename MeshEntity >
+__cuda_callable__
+Index
+LaxFridrichs< tnlGrid< 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
+LaxFridrichs< tnlGrid< 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 >
+tnlString
+LaxFridrichs< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+ return tnlString( "LaxFridrichs< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename MeshFunction, typename MeshEntity >
+__cuda_callable__
+Real
+LaxFridrichs< tnlGrid< 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
+LaxFridrichs< tnlGrid< 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
+LaxFridrichs< tnlGrid< 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 );
+}
+
+#endif /* LaxFridrichsIMPL_H */
+
diff --git a/examples/advection/advection.cpp b/examples/advection/advection.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..374d4714086168f7ebf1ed2a62664c4aaae0c4b7
--- /dev/null
+++ b/examples/advection/advection.cpp
@@ -0,0 +1 @@
+#include "advection.h"
diff --git a/examples/advection/advection.cu b/examples/advection/advection.cu
new file mode 100644
index 0000000000000000000000000000000000000000..374d4714086168f7ebf1ed2a62664c4aaae0c4b7
--- /dev/null
+++ b/examples/advection/advection.cu
@@ -0,0 +1 @@
+#include "advection.h"
diff --git a/examples/advection/advection.h b/examples/advection/advection.h
new file mode 100644
index 0000000000000000000000000000000000000000..7c722df204506bd0977dd5eec1b231e84fd1f8e2
--- /dev/null
+++ b/examples/advection/advection.h
@@ -0,0 +1,116 @@
+#include <tnlConfig.h>
+#include <solvers/tnlSolver.h>
+#include <solvers/tnlBuildConfigTags.h>
+#include <operators/tnlDirichletBoundaryConditions.h>
+#include <operators/tnlNeumannBoundaryConditions.h>
+#include <functions/tnlConstantFunction.h>
+#include "advectionProblem.h"
+#include "LaxFridrichs.h"
+#include "advectionRhs.h"
+#include "advectionBuildConfigTag.h"
+
+typedef advectionBuildConfigTag 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 advectionConfig
+{
+ public:
+ static void configSetup( tnlConfigDescription & config )
+ {
+ config.addDelimiter( "advection settings:" );
+ config.addEntry< tnlString >( "boundary-conditions-type", "Choose the boundary conditions type.", "dirichlet");
+ config.addEntryEnum< tnlString >( "dirichlet" );
+ config.addEntryEnum< tnlString >( "neumann" );
+ config.addEntry< double >( "boundary-conditions-constant", "This sets a value in case of the constant boundary conditions." );
+ config.addEntry< double >( "artifical-viscosity", "This sets value of artifical viscosity (default 1)", 1.0);
+ config.addEntry< tnlString >( "begin", "choose begin type", "sin");
+ config.addEntryEnum< tnlString >( "sin");
+ config.addEntryEnum< tnlString >( "sin_square");
+ config.addEntry< double >( "advection-speedX", "This sets value of advection speed in X direction (default 1)" , 1.0);
+ config.addEntry< double >( "advection-speedY", "This sets value of advection speed in Y direction (default 1)" , 1.0);
+ config.addEntry< tnlString >( "move", "choose movement type", "advection");
+ config.addEntryEnum< tnlString >( "advection");
+ config.addEntryEnum< tnlString >( "rotation");
+
+ /****
+ * Add definition of your solver command line arguments.
+ */
+
+ }
+};
+
+template< typename Real,
+ typename Device,
+ typename Index,
+ typename MeshType,
+ typename ConfigTag,
+ typename SolverStarter >
+class advectionSetter
+{
+ public:
+
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ static bool run( const tnlParameterContainer & parameters )
+ {
+ enum { Dimensions = MeshType::getMeshDimensions() };
+ typedef LaxFridrichs< MeshType, Real, Index > ApproximateOperator;
+ typedef advectionRhs< MeshType, Real > RightHandSide;
+ typedef tnlStaticVector < 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 tnlVector.
+ */
+ tnlString boundaryConditionsType = parameters.getParameter< tnlString >( "boundary-conditions-type" );
+ if( parameters.checkParameter( "boundary-conditions-constant" ) )
+ {
+ typedef tnlConstantFunction< Dimensions, Real > ConstantFunction;
+ if( boundaryConditionsType == "dirichlet" )
+ {
+ typedef tnlDirichletBoundaryConditions< MeshType, ConstantFunction, MeshType::getMeshDimensions(), Real, Index > BoundaryConditions;
+ typedef advectionProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
+ SolverStarter solverStarter;
+ return solverStarter.template run< Problem >( parameters );
+ }
+ typedef tnlNeumannBoundaryConditions< MeshType, ConstantFunction, Real, Index > BoundaryConditions;
+ typedef advectionProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
+ SolverStarter solverStarter;
+ return solverStarter.template run< Problem >( parameters );
+ }
+ typedef tnlMeshFunction< MeshType > MeshFunction;
+ if( boundaryConditionsType == "dirichlet" )
+ {
+ typedef tnlDirichletBoundaryConditions< MeshType, MeshFunction, MeshType::getMeshDimensions(), Real, Index > BoundaryConditions;
+ typedef advectionProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
+ SolverStarter solverStarter;
+ return solverStarter.template run< Problem >( parameters );
+ }
+ typedef tnlNeumannBoundaryConditions< MeshType, MeshFunction, Real, Index > BoundaryConditions;
+ typedef advectionProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
+ SolverStarter solverStarter;
+ return solverStarter.template run< Problem >( parameters );
+ }
+
+};
+
+int main( int argc, char* argv[] )
+{
+ tnlSolver< advectionSetter, advectionConfig, BuildConfig > solver;
+ if( ! solver. run( argc, argv ) )
+ return EXIT_FAILURE;
+ return EXIT_SUCCESS;
+}
+
diff --git a/examples/advection/advectionBuildConfigTag.h b/examples/advection/advectionBuildConfigTag.h
new file mode 100644
index 0000000000000000000000000000000000000000..a8ea6ddef4eba8fa8802a3ae41f7ceefe1628273
--- /dev/null
+++ b/examples/advection/advectionBuildConfigTag.h
@@ -0,0 +1,43 @@
+#ifndef advectionBUILDCONFIGTAG_H_
+#define advectionBUILDCONFIGTAG_H_
+
+#include <solvers/tnlBuildConfigTags.h>
+
+class advectionBuildConfigTag{};
+
+/****
+ * Turn off support for float and long double.
+ */
+template<> struct tnlConfigTagReal< advectionBuildConfigTag, float > { enum { enabled = false }; };
+template<> struct tnlConfigTagReal< advectionBuildConfigTag, long double > { enum { enabled = false }; };
+
+/****
+ * Turn off support for short int and long int indexing.
+ */
+template<> struct tnlConfigTagIndex< advectionBuildConfigTag, short int >{ enum { enabled = false }; };
+template<> struct tnlConfigTagIndex< advectionBuildConfigTag, long int >{ enum { enabled = false }; };
+
+/****
+ * Use of tnlGrid is enabled for allowed dimensions and Real, Device and Index types.
+ */
+
+template< int Dimensions, typename Real, typename Device, typename Index >
+ struct tnlConfigTagMesh< advectionBuildConfigTag, tnlGrid< Dimensions, Real, Device, Index > >
+ { enum { enabled = tnlConfigTagDimensions< advectionBuildConfigTag, Dimensions >::enabled &&
+ tnlConfigTagReal< advectionBuildConfigTag, Real >::enabled &&
+ tnlConfigTagDevice< advectionBuildConfigTag, Device >::enabled &&
+ tnlConfigTagIndex< advectionBuildConfigTag, Index >::enabled }; };
+
+/****
+ * Please, chose your preferred time discretisation here.
+ */
+template<> struct tnlConfigTagTimeDiscretisation< advectionBuildConfigTag, tnlExplicitTimeDiscretisationTag >{ enum { enabled = true }; };
+template<> struct tnlConfigTagTimeDiscretisation< advectionBuildConfigTag, tnlSemiImplicitTimeDiscretisationTag >{ enum { enabled = true }; };
+template<> struct tnlConfigTagTimeDiscretisation< advectionBuildConfigTag, tnlImplicitTimeDiscretisationTag >{ enum { enabled = true }; };
+
+/****
+ * Only the Runge-Kutta-Merson solver is enabled by default.
+ */
+template<> struct tnlConfigTagExplicitSolver< advectionBuildConfigTag, tnlExplicitEulerSolverTag >{ enum { enabled = true }; };
+
+#endif /* advectionBUILDCONFIGTAG_H_ */
diff --git a/examples/advection/advectionProblem.h b/examples/advection/advectionProblem.h
new file mode 100644
index 0000000000000000000000000000000000000000..64ae07480c584cacf3647d09c5bf075ad0f67ce8
--- /dev/null
+++ b/examples/advection/advectionProblem.h
@@ -0,0 +1,83 @@
+#ifndef advectionPROBLEM_H_
+#define advectionPROBLEM_H_
+
+#include <problems/tnlPDEProblem.h>
+#include <functions/tnlMeshFunction.h>
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+class advectionProblem:
+ public tnlPDEProblem< 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 tnlMeshFunction< Mesh > MeshFunctionType;
+ typedef tnlPDEProblem< Mesh, RealType, DeviceType, IndexType > BaseType;
+
+ using typename BaseType::MeshType;
+ using typename BaseType::DofVectorType;
+ using typename BaseType::MeshDependentDataType;
+ tnlString velocityType;
+ static tnlString getTypeStatic();
+
+ tnlString getPrologHeader() const;
+
+ void writeProlog( tnlLogger& logger,
+ const tnlParameterContainer& parameters ) const;
+
+ bool setup( const tnlParameterContainer& parameters );
+
+ bool setInitialCondition( const tnlParameterContainer& parameters,
+ const MeshType& mesh,
+ DofVectorType& dofs,
+ MeshDependentDataType& meshDependentData );
+
+ template< typename Matrix >
+ bool setupLinearSystem( const MeshType& mesh,
+ Matrix& matrix );
+
+ bool makeSnapshot( const RealType& time,
+ const IndexType& step,
+ const MeshType& mesh,
+ DofVectorType& dofs,
+ MeshDependentDataType& meshDependentData );
+
+ IndexType getDofs( const MeshType& mesh ) const;
+
+ void bindDofs( const MeshType& mesh,
+ DofVectorType& dofs );
+
+ void getExplicitRHS( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ DofVectorType& _u,
+ DofVectorType& _fu,
+ MeshDependentDataType& meshDependentData );
+
+ template< typename Matrix >
+ void assemblyLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ DofVectorType& dofs,
+ Matrix& matrix,
+ DofVectorType& rightHandSide,
+ MeshDependentDataType& meshDependentData );
+
+ protected:
+
+ DifferentialOperator differentialOperator;
+ BoundaryCondition boundaryCondition;
+ RightHandSide rightHandSide;
+};
+
+#include "advectionProblem_impl.h"
+
+#endif /* advectionPROBLEM_H_ */
diff --git a/examples/advection/advectionProblem_impl.h b/examples/advection/advectionProblem_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..7d7c324d5ae1d6d2b90ee73827101743461bba2a
--- /dev/null
+++ b/examples/advection/advectionProblem_impl.h
@@ -0,0 +1,317 @@
+#ifndef advectionPROBLEM_IMPL_H_
+#define advectionPROBLEM_IMPL_H_
+
+#include <core/mfilename.h>
+#include <matrices/tnlMatrixSetter.h>
+#include <solvers/pde/tnlExplicitUpdater.h>
+#include <solvers/pde/tnlLinearSystemAssembler.h>
+#include <solvers/pde/tnlBackwardTimeDiscretisation.h>
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+tnlString
+advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+getTypeStatic()
+{
+ return tnlString( "advectionProblem< " ) + Mesh :: getTypeStatic() + " >";
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+tnlString
+advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+getPrologHeader() const
+{
+ return tnlString( "advection" );
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+void
+advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+writeProlog( tnlLogger& logger, const tnlParameterContainer& 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
+advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+setup( const tnlParameterContainer& parameters )
+{
+ if( ! this->boundaryCondition.setup( parameters, "boundary-conditions-" ) ||
+ ! this->rightHandSide.setup( parameters, "right-hand-side-" ) )
+ return false;
+ return true;
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+typename advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::IndexType
+advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+getDofs( const MeshType& mesh ) const
+{
+ /****
+ * Return number of DOFs (degrees of freedom) i.e. number
+ * of unknowns to be resolved by the main solver.
+ */
+ return mesh.template getEntitiesCount< typename MeshType::Cell >();
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+void
+advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+bindDofs( const MeshType& mesh,
+ DofVectorType& dofVector )
+{
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+bool
+advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+setInitialCondition( const tnlParameterContainer& parameters,
+ const MeshType& mesh,
+ DofVectorType& dofs,
+ MeshDependentDataType& meshDependentData )
+{
+ typedef typename MeshType::Cell Cell;
+ int count = mesh.template getEntitiesCount< Cell >();
+ const RealType& size = 1;//mesh.template getSpaceStepsProducts< -1, 0 >();
+ const tnlString& beginChoice = parameters.getParameter< tnlString >( "begin" );
+ int dimensions = 2; //vyresit!!!!
+ if (beginChoice == "sin_square")
+ {
+ double constantFunction;
+ if (dimensions == 1)
+ {
+ dofs[0] = 0;
+ double expValue;
+ for (long int i = 1; i < count; i++)
+ {
+ expValue = exp(-pow(size*i-2,2));
+ if ((i>0.4*count) && (i<0.5*count)) constantFunction=1; else constantFunction=0;
+ if (expValue>constantFunction) dofs[i] = expValue; else dofs[i] = constantFunction;
+ };
+ dofs[count] = 0;
+ };
+ if (dimensions == 2)
+ {
+ double expValue;
+ for (long int i = 1; i < count; i++)
+ for (long int j = 1; j < count; j++)
+ {
+ expValue = exp(-pow(size*i-2,2)-pow(size*j-2,2));
+ if ((i>0.4*count) && (i<0.5*count)&&((j>0.4*count) && (j<0.5*count))) constantFunction=1; else constantFunction=0;
+ if (expValue>constantFunction) dofs[(i-1) * count -1+ j] = expValue; else dofs[(i-1) * count -1+ j] = constantFunction;
+ };
+ };
+ }
+ else if (beginChoice == "sin")
+ {
+ if (dimensions == 1)
+ {
+ dofs[0] = 0;
+ for (long int i = 1; i < count; i++)
+ {
+ dofs[i] = exp(-pow(size*i-2,2));
+ };
+ dofs[count] = 0;
+ };
+ if (dimensions == 2)
+ {
+ for (long int i = 1; i < count; i++)
+ for (long int j = 1; j < count; j++)
+ {
+ dofs[(i-1) * count -1+ j] = exp(-pow(size*i-2,2)-pow(size*j-2,2));
+ };
+ };
+ };
+//setting velocity field
+
+
+ /*const tnlString& initialConditionFile = parameters.getParameter< tnlString >( "initial-condition" );
+ if( ! dofs.load( initialConditionFile ) )
+ {
+ cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." << endl;
+ return false;
+ }
+ return true;*/
+ dofs.save( "dofs.tnl" );
+ this->velocityType = parameters.getParameter< tnlString >( "move" );
+ const double artificalViscosity = parameters.getParameter< double >( "artifical-viscosity" );
+ differentialOperator.setViscosity(artificalViscosity);
+ const double advectionSpeedX = parameters.getParameter< double >( "advection-speedX" );
+ differentialOperator.setAdvectionSpeedX(advectionSpeedX);
+ const double advectionSpeedY = parameters.getParameter< double >( "advection-speedY" );
+ differentialOperator.setAdvectionSpeedY(advectionSpeedY);
+ return true;
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+ template< typename Matrix >
+bool
+advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+setupLinearSystem( const MeshType& mesh,
+ Matrix& matrix )
+{
+ const IndexType dofs = this->getDofs( mesh );
+ typedef typename Matrix::CompressedRowsLengthsVector CompressedRowsLengthsVectorType;
+ CompressedRowsLengthsVectorType rowLengths;
+ if( ! rowLengths.setSize( dofs ) )
+ return false;
+ tnlMatrixSetter< 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
+advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+makeSnapshot( const RealType& time,
+ const IndexType& step,
+ const MeshType& mesh,
+ DofVectorType& dofs,
+ MeshDependentDataType& meshDependentData )
+{
+ cout << endl << "Writing output at time " << time << " step " << step << "." << endl;
+ this->bindDofs( mesh, dofs );
+ tnlString fileName;
+ FileNameBaseNumberEnding( "u-", step, 5, ".tnl", fileName );
+ if( ! dofs.save( fileName ) )
+ return false;
+ return true;
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+void
+advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+getExplicitRHS( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ DofVectorType& _u,
+ DofVectorType& _fu,
+ MeshDependentDataType& meshDependentData )
+{
+ /****
+ * If you use an explicit solver like tnlEulerSolver or tnlMersonSolver, you
+ * need to implement this method. Compute the right-hand side of
+ *
+ * d/dt u(x) = fu( x, u )
+ *
+ * You may use supporting mesh dependent data if you need.
+ */
+ typedef typename MeshType::Cell Cell;
+ int count = mesh.template getEntitiesCount< Cell >();
+ const RealType& size = 1;//mesh.template getSpaceStepsProducts< -1, 0 >();
+ if (this->velocityType == "rotation")
+ {
+ double radius;
+ for (int i =1; i < count; i++)
+ for (int j =1; j < count; j++)
+ {
+ radius = sqrt(pow(i-1-(count/2.0),2) + pow(j-1-(count/2.0),2));
+ if (radius != 0.0)
+ _fu[(i-1)*count+j-1] =(0.25*tau)*differentialOperator.artificalViscosity* //smoothening part
+ (_u[(i-1)*count-2+j]+_u[(i-1)*count+j]+
+ _u[i*count+j-1]+_u[(i-2)*count+j-1]-
+ 4.0*_u[(i-1)*count+j-1])
+ -((1.0/(2.0*count))*differentialOperator.advectionSpeedX //X addition
+ *radius*(-1)*((j-1-(count/2.0))/radius)
+ *(_u[(i-1)*count+j]-_u[(i-1)*count+j-2]))
+ -((1.0/(2.0*count))*differentialOperator.advectionSpeedY //Y addition
+ *radius*((i-1-(count/2.0))/radius)
+ *(_u[i*count+j-1]-_u[(i-2)*count+j-1]))
+ ;}
+ }
+ else if (this->velocityType == "advection")
+ {
+ this->bindDofs( mesh, _u );
+ tnlExplicitUpdater< Mesh, MeshFunctionType, DifferentialOperator, BoundaryCondition, RightHandSide > explicitUpdater;
+ MeshFunctionType u( mesh, _u );
+ MeshFunctionType fu( mesh, _fu );
+ explicitUpdater.template update< typename Mesh::Cell >( time,
+ mesh,
+ this->differentialOperator,
+ this->boundaryCondition,
+ this->rightHandSide,
+ u,
+ fu );
+ tnlBoundaryConditionsSetter< MeshFunctionType, BoundaryCondition > boundaryConditionsSetter;
+ boundaryConditionsSetter.template apply< typename Mesh::Cell >(
+ this->boundaryCondition,
+ time + tau,
+ u );
+ }
+}
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+ template< typename Matrix >
+void
+advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+assemblyLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ DofVectorType& _u,
+ Matrix& matrix,
+ DofVectorType& b,
+ MeshDependentDataType& meshDependentData )
+{
+ tnlLinearSystemAssembler< Mesh,
+ MeshFunctionType,
+ DifferentialOperator,
+ BoundaryCondition,
+ RightHandSide,
+ tnlBackwardTimeDiscretisation,
+ Matrix,
+ DofVectorType > systemAssembler;
+
+ tnlMeshFunction< Mesh > u( mesh, _u );
+ systemAssembler.template assembly< typename Mesh::Cell >( time,
+ tau,
+ mesh,
+ this->differentialOperator,
+ this->boundaryCondition,
+ this->rightHandSide,
+ u,
+ matrix,
+ b );
+}
+
+#endif /* advectionPROBLEM_IMPL_H_ */
diff --git a/examples/advection/advectionRhs.h b/examples/advection/advectionRhs.h
new file mode 100644
index 0000000000000000000000000000000000000000..3bd2c294e142cac0779f66bdc99da51dc572222b
--- /dev/null
+++ b/examples/advection/advectionRhs.h
@@ -0,0 +1,29 @@
+#ifndef advectionRHS_H_
+#define advectionRHS_H_
+#include<functions/tnlDomain.h>
+template< typename Mesh, typename Real >class advectionRhs
+ : public tnlDomain< Mesh::meshDimensions, MeshDomain >
+ {
+ public:
+
+ typedef Mesh MeshType;
+ typedef Real RealType;
+
+ bool setup( const tnlParameterContainer& parameters,
+ const tnlString& 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;
+ };
+};
+
+#endif /* advectionRHS_H_ */
diff --git a/examples/advection/tnl-run-advection b/examples/advection/tnl-run-advection
new file mode 100644
index 0000000000000000000000000000000000000000..97122dd9ca95ecc0da694bd54a9b97ab8c383310
--- /dev/null
+++ b/examples/advection/tnl-run-advection
@@ -0,0 +1,20 @@
+#!/usr/bin/env bash
+
+tnl-grid-setup --dimensions 1 \
+ --origin-x 0.0 \
+ --proportions-x 10.0 \
+ --size-x 500 \
+
+#tnl-init --test-function sin-wave \
+# --output-file init.tnl
+./advection --time-discretisation explicit \
+ --time-step 1.0e-3 \
+ --boundary-conditions-constant 0 \
+ --discrete-solver euler \
+ --snapshot-period 0.1 \
+ --final-time 1.0 \
+ --artifical-viscosity 1.0 \
+ --begin sin \
+ --advection-speedX 2.0 \
+
+tnl-view --mesh mesh.tnl --input-files *tnl
diff --git a/examples/advection/tnl-run-advection1 b/examples/advection/tnl-run-advection1
new file mode 100644
index 0000000000000000000000000000000000000000..fb5e81faf2e19ab8e6b19c02c0c441dd96a9c807
--- /dev/null
+++ b/examples/advection/tnl-run-advection1
@@ -0,0 +1,24 @@
+#!/usr/bin/env bash
+
+tnl-grid-setup --dimensions 2 \
+ --origin-x 0.0 \
+ --origin-y 0.0 \
+ --proportions-x 10.0 \
+ --proportions-y 10.0 \
+ --size-x 500 \
+ --size-y 500 \
+
+#tnl-init --test-function sin-wave \
+# --output-file init.tnl
+./advection --time-discretisation explicit \
+ --time-step 1.0e-3 \
+ --boundary-conditions-constant 0 \
+ --discrete-solver euler \
+ --snapshot-period 0.1 \
+ --final-time 1.0 \
+ --artifical-viscosity 0.2 \
+ --begin sin_square \
+ --advection-speedX 2.0 \
+ --advection-speedY 2.0 \
+
+tnl-view --mesh mesh.tnl --input-files *tnl
diff --git a/examples/advection/tnl-run-advectionrot1 b/examples/advection/tnl-run-advectionrot1
new file mode 100644
index 0000000000000000000000000000000000000000..5d3e3990a98ef64233a6b629b4002cc50b26a923
--- /dev/null
+++ b/examples/advection/tnl-run-advectionrot1
@@ -0,0 +1,25 @@
+#!/usr/bin/env bash
+
+tnl-grid-setup --dimensions 2 \
+ --origin-x 0.0 \
+ --origin-y 0.0 \
+ --proportions-x 10.0 \
+ --proportions-y 10.0 \
+ --size-x 100 \
+ --size-y 100 \
+
+#tnl-init --test-function sin-wave \
+# --output-file init.tnl
+./advection --time-discretisation explicit \
+ --time-step 1.0e-2 \
+ --boundary-conditions-constant 0 \
+ --discrete-solver euler \
+ --snapshot-period 0.1 \
+ --final-time 1.0 \
+ --artifical-viscosity 5.0 \
+ --begin sin_square \
+ --advection-speedX 2.0 \
+ --advection-speedY 2.0 \
+ --move rotation \
+
+tnl-view --mesh mesh.tnl --input-files *tnl
diff --git a/examples/inviscid-flow/1d/CMakeLists.txt b/examples/inviscid-flow/1d/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..6fe0766eec90a76d7888c56f615122804abf4bcc
--- /dev/null
+++ b/examples/inviscid-flow/1d/CMakeLists.txt
@@ -0,0 +1,20 @@
+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/LaxFridrichs.h b/examples/inviscid-flow/1d/LaxFridrichs.h
new file mode 100644
index 0000000000000000000000000000000000000000..db63defd3f6585702840e7bb575e225c6667c917
--- /dev/null
+++ b/examples/inviscid-flow/1d/LaxFridrichs.h
@@ -0,0 +1,143 @@
+#ifndef LaxFridrichs_H
+#define LaxFridrichs_H
+
+#include <core/vectors/tnlVector.h>
+#include <mesh/tnlGrid.h>
+
+template< typename Mesh,
+ typename Real = typename Mesh::RealType,
+ typename Index = typename Mesh::IndexType >
+class LaxFridrichs
+{
+};
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class LaxFridrichs< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index >
+{
+ public:
+ typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef tnlMeshFunction< MeshType > MeshFunctionType;
+ enum { Dimensions = MeshType::getMeshDimensions() };
+
+ static tnlString getType();
+
+ template< typename MeshFunction, typename MeshEntity >
+ __cuda_callable__
+ Real operator()( const MeshFunction& u,
+ const MeshEntity& entity,
+ const RealType& time = 0.0 ) const;
+
+ __cuda_callable__
+ template< typename MeshEntity >
+ 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 LaxFridrichs< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index >
+{
+ public:
+ typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef tnlMeshFunction< MeshType > MeshFunctionType;
+ enum { Dimensions = MeshType::getMeshDimensions() };
+
+ static tnlString getType();
+
+ template< typename MeshFunction, typename MeshEntity >
+ __cuda_callable__
+ Real operator()( const MeshFunction& u,
+ const MeshEntity& entity,
+ const RealType& time = 0.0 ) const;
+
+ __cuda_callable__
+ template< typename MeshEntity >
+ 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 LaxFridrichs< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, Index >
+{
+ public:
+ typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef tnlMeshFunction< MeshType > MeshFunctionType;
+ enum { Dimensions = MeshType::getMeshDimensions() };
+
+ static tnlString getType();
+
+ template< typename MeshFunction, typename MeshEntity >
+ __cuda_callable__
+ Real operator()( const MeshFunction& u,
+ const MeshEntity& entity,
+ const RealType& time = 0.0 ) const;
+
+ __cuda_callable__
+ template< typename MeshEntity >
+ 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;
+};
+
+
+#include "LaxFridrichs_impl.h"
+
+#endif /* LaxFridrichs_H */
diff --git a/examples/inviscid-flow/1d/LaxFridrichs_impl.h b/examples/inviscid-flow/1d/LaxFridrichs_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..c1fa87b34c637da6edfb77e0fafda1a4c47c8399
--- /dev/null
+++ b/examples/inviscid-flow/1d/LaxFridrichs_impl.h
@@ -0,0 +1,344 @@
+#ifndef LaxFridrichs_IMPL_H
+#define LaxFridrichs_IMPL_H
+
+/****
+ * 1D problem
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString
+LaxFridrichs< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+ return tnlString( "LaxFridrichs< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename MeshFunction, typename MeshEntity >
+__cuda_callable__
+Real
+LaxFridrichs< tnlGrid< 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< -2 >();
+ 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
+LaxFridrichs< tnlGrid< 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
+LaxFridrichs< tnlGrid< 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 >
+tnlString
+LaxFridrichs< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+ return tnlString( "LaxFridrichs< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename MeshFunction, typename MeshEntity >
+__cuda_callable__
+Real
+LaxFridrichs< tnlGrid< 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
+LaxFridrichs< tnlGrid< 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
+LaxFridrichs< tnlGrid< 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 >
+tnlString
+LaxFridrichs< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+ return tnlString( "LaxFridrichs< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename MeshFunction, typename MeshEntity >
+__cuda_callable__
+Real
+LaxFridrichs< tnlGrid< 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
+LaxFridrichs< tnlGrid< 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
+LaxFridrichs< tnlGrid< 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 );
+}
+
+#endif /* LaxFridrichsIMPL_H */
+
diff --git a/examples/inviscid-flow/1d/euler.cpp b/examples/inviscid-flow/1d/euler.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4d76005cb1f70724be978ff0fa6fec63c4a8a76f
--- /dev/null
+++ b/examples/inviscid-flow/1d/euler.cpp
@@ -0,0 +1 @@
+#include "euler.h"
diff --git a/examples/inviscid-flow/1d/euler.cu b/examples/inviscid-flow/1d/euler.cu
new file mode 100644
index 0000000000000000000000000000000000000000..4d76005cb1f70724be978ff0fa6fec63c4a8a76f
--- /dev/null
+++ b/examples/inviscid-flow/1d/euler.cu
@@ -0,0 +1 @@
+#include "euler.h"
diff --git a/examples/inviscid-flow/1d/euler.h b/examples/inviscid-flow/1d/euler.h
new file mode 100644
index 0000000000000000000000000000000000000000..cba3c12b65fa683360aa5076e861b7e69d6ec16a
--- /dev/null
+++ b/examples/inviscid-flow/1d/euler.h
@@ -0,0 +1,115 @@
+#include <tnlConfig.h>
+#include <solvers/tnlSolver.h>
+#include <solvers/tnlBuildConfigTags.h>
+#include <operators/tnlDirichletBoundaryConditions.h>
+#include <operators/tnlNeumannBoundaryConditions.h>
+#include <functions/tnlConstantFunction.h>
+#include "eulerProblem.h"
+#include "LaxFridrichs.h"
+#include "eulerRhs.h"
+#include "eulerBuildConfigTag.h"
+
+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( tnlConfigDescription & config )
+ {
+ config.addDelimiter( "euler settings:" );
+ config.addEntry< tnlString >( "boundary-conditions-type", "Choose the boundary conditions type.", "dirichlet");
+ config.addEntryEnum< tnlString >( "dirichlet" );
+ config.addEntryEnum< tnlString >( "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 tnlParameterContainer & parameters )
+ {
+ enum { Dimensions = MeshType::getMeshDimensions() };
+ typedef LaxFridrichs< MeshType, Real, Index > ApproximateOperator;
+ typedef eulerRhs< MeshType, Real > RightHandSide;
+ typedef tnlStaticVector < 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 tnlVector.
+ */
+ tnlString boundaryConditionsType = parameters.getParameter< tnlString >( "boundary-conditions-type" );
+ if( parameters.checkParameter( "boundary-conditions-constant" ) )
+ {
+ typedef tnlConstantFunction< Dimensions, Real > ConstantFunction;
+ if( boundaryConditionsType == "dirichlet" )
+ {
+ typedef tnlDirichletBoundaryConditions< MeshType, ConstantFunction, MeshType::getMeshDimensions(), Real, Index > BoundaryConditions;
+ typedef eulerProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
+ SolverStarter solverStarter;
+ return solverStarter.template run< Problem >( parameters );
+ }
+ typedef tnlNeumannBoundaryConditions< MeshType, ConstantFunction, Real, Index > BoundaryConditions;
+ typedef eulerProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
+ SolverStarter solverStarter;
+ return solverStarter.template run< Problem >( parameters );
+ }
+ typedef tnlMeshFunction< MeshType > MeshFunction;
+ if( boundaryConditionsType == "dirichlet" )
+ {
+ typedef tnlDirichletBoundaryConditions< MeshType, MeshFunction, MeshType::getMeshDimensions(), Real, Index > BoundaryConditions;
+ typedef eulerProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
+ SolverStarter solverStarter;
+ return solverStarter.template run< Problem >( parameters );
+ }
+ typedef tnlNeumannBoundaryConditions< 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[] )
+{
+ tnlSolver< 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
new file mode 100644
index 0000000000000000000000000000000000000000..dea2c8445e7bf7dfbd537e1b98ca40b669b674b6
--- /dev/null
+++ b/examples/inviscid-flow/1d/eulerBuildConfigTag.h
@@ -0,0 +1,43 @@
+#ifndef eulerBUILDCONFIGTAG_H_
+#define eulerBUILDCONFIGTAG_H_
+
+#include <solvers/tnlBuildConfigTags.h>
+
+class eulerBuildConfigTag{};
+
+/****
+ * Turn off support for float and long double.
+ */
+template<> struct tnlConfigTagReal< eulerBuildConfigTag, float > { enum { enabled = false }; };
+template<> struct tnlConfigTagReal< eulerBuildConfigTag, long double > { enum { enabled = false }; };
+
+/****
+ * Turn off support for short int and long int indexing.
+ */
+template<> struct tnlConfigTagIndex< eulerBuildConfigTag, short int >{ enum { enabled = false }; };
+template<> struct tnlConfigTagIndex< eulerBuildConfigTag, long int >{ enum { enabled = false }; };
+
+/****
+ * Use of tnlGrid is enabled for allowed dimensions and Real, Device and Index types.
+ */
+
+template< int Dimensions, typename Real, typename Device, typename Index >
+ struct tnlConfigTagMesh< eulerBuildConfigTag, tnlGrid< Dimensions, Real, Device, Index > >
+ { enum { enabled = tnlConfigTagDimensions< eulerBuildConfigTag, Dimensions >::enabled &&
+ tnlConfigTagReal< eulerBuildConfigTag, Real >::enabled &&
+ tnlConfigTagDevice< eulerBuildConfigTag, Device >::enabled &&
+ tnlConfigTagIndex< eulerBuildConfigTag, Index >::enabled }; };
+
+/****
+ * Please, chose your preferred time discretisation here.
+ */
+template<> struct tnlConfigTagTimeDiscretisation< eulerBuildConfigTag, tnlExplicitTimeDiscretisationTag >{ enum { enabled = true }; };
+template<> struct tnlConfigTagTimeDiscretisation< eulerBuildConfigTag, tnlSemiImplicitTimeDiscretisationTag >{ enum { enabled = true }; };
+template<> struct tnlConfigTagTimeDiscretisation< eulerBuildConfigTag, tnlImplicitTimeDiscretisationTag >{ enum { enabled = false }; };
+
+/****
+ * Only the Runge-Kutta-Merson solver is enabled by default.
+ */
+template<> struct tnlConfigTagExplicitSolver< eulerBuildConfigTag, tnlExplicitEulerSolverTag >{ enum { enabled = false }; };
+
+#endif /* eulerBUILDCONFIGTAG_H_ */
diff --git a/examples/inviscid-flow/1d/eulerProblem.h b/examples/inviscid-flow/1d/eulerProblem.h
new file mode 100644
index 0000000000000000000000000000000000000000..8c2e5020fc2417de5038ffbf73ba1f8f5b62c09a
--- /dev/null
+++ b/examples/inviscid-flow/1d/eulerProblem.h
@@ -0,0 +1,90 @@
+#ifndef eulerPROBLEM_H_
+#define eulerPROBLEM_H_
+
+#include <problems/tnlPDEProblem.h>
+#include <functions/tnlMeshFunction.h>
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+class eulerProblem:
+ public tnlPDEProblem< 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 tnlMeshFunction< Mesh > MeshFunctionType;
+ typedef tnlPDEProblem< Mesh, RealType, DeviceType, IndexType > BaseType;
+
+ using typename BaseType::MeshType;
+ using typename BaseType::DofVectorType;
+ using typename BaseType::MeshDependentDataType;
+
+ static tnlString getTypeStatic();
+ tnlSharedVector< RealType, DeviceType, IndexType > rho;
+ tnlSharedVector< RealType, DeviceType, IndexType > rhoVel;
+ tnlSharedVector< RealType, DeviceType, IndexType > energy;
+ tnlVector< RealType, DeviceType, IndexType > data;
+ tnlSharedVector< RealType, DeviceType, IndexType > pressure;
+ tnlSharedVector< RealType, DeviceType, IndexType > velocity;
+ double gamma;
+
+ tnlString getPrologHeader() const;
+
+ void writeProlog( tnlLogger& logger,
+ const tnlParameterContainer& parameters ) const;
+
+ bool setup( const tnlParameterContainer& parameters );
+
+ bool setInitialCondition( const tnlParameterContainer& parameters,
+ const MeshType& mesh,
+ DofVectorType& dofs,
+ MeshDependentDataType& meshDependentData );
+
+ template< typename Matrix >
+ bool setupLinearSystem( const MeshType& mesh,
+ Matrix& matrix );
+
+ bool makeSnapshot( const RealType& time,
+ const IndexType& step,
+ const MeshType& mesh,
+ DofVectorType& dofs,
+ MeshDependentDataType& meshDependentData );
+
+ IndexType getDofs( const MeshType& mesh ) const;
+
+ void bindDofs( const MeshType& mesh,
+ DofVectorType& dofs );
+
+ void getExplicitRHS( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ DofVectorType& _u,
+ DofVectorType& _fu,
+ MeshDependentDataType& meshDependentData );
+
+ template< typename Matrix >
+ void assemblyLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ DofVectorType& dofs,
+ Matrix& matrix,
+ DofVectorType& rightHandSide,
+ MeshDependentDataType& meshDependentData );
+
+ protected:
+
+ DifferentialOperator differentialOperator;
+ BoundaryCondition boundaryCondition;
+ RightHandSide rightHandSide;
+};
+
+#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
new file mode 100644
index 0000000000000000000000000000000000000000..06c647f9594586c9181c37b955ffe8d360456504
--- /dev/null
+++ b/examples/inviscid-flow/1d/eulerProblem_impl.h
@@ -0,0 +1,338 @@
+#ifndef eulerPROBLEM_IMPL_H_
+#define eulerPROBLEM_IMPL_H_
+
+#include <core/mfilename.h>
+#include <matrices/tnlMatrixSetter.h>
+#include <solvers/pde/tnlExplicitUpdater.h>
+#include <solvers/pde/tnlLinearSystemAssembler.h>
+#include <solvers/pde/tnlBackwardTimeDiscretisation.h>
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+tnlString
+eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+getTypeStatic()
+{
+ return tnlString( "eulerProblem< " ) + Mesh :: getTypeStatic() + " >";
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+tnlString
+eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+getPrologHeader() const
+{
+ return tnlString( "euler" );
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+void
+eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+writeProlog( tnlLogger& logger, const tnlParameterContainer& 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 tnlParameterContainer& parameters )
+{
+ if( ! this->boundaryCondition.setup( parameters, "boundary-conditions-" ) ||
+ ! this->rightHandSide.setup( parameters, "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 MeshType& 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 MeshType& mesh,
+ DofVectorType& dofVector )
+{
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+bool
+eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+setInitialCondition( const tnlParameterContainer& parameters,
+ const MeshType& mesh,
+ DofVectorType& dofs,
+ MeshDependentDataType& meshDependentData )
+{
+ typedef typename MeshType::Cell Cell;
+ double gamma = parameters.getParameter< double >( "gamma" );
+ double rhoL = parameters.getParameter< double >( "left-density" );
+ double velL = parameters.getParameter< double >( "left-velocity" );
+ double preL = parameters.getParameter< double >( "left-pressure" );
+ double eL = ( preL / (gamma - 1) ) + 0.5 * rhoL * velL * velL;
+ double rhoR = parameters.getParameter< double >( "right-density" );
+ double velR = parameters.getParameter< double >( "right-velocity" );
+ double preR = parameters.getParameter< double >( "right-pressure" );
+ double eR = ( preR / (gamma - 1) ) + 0.5 * rhoR * velR * velR;
+ double x0 = parameters.getParameter< double >( "riemann-border" );
+ cout << gamma << " " << rhoL << " " << velL << " " << preL << " " << eL << " " << rhoR << " " << velR << " " << preR << " " << eR << " " << x0 << " " << gamma << endl;
+ int count = mesh.template getEntitiesCount< Cell >();
+ this->rho.bind(dofs,0,count);
+ this->rhoVel.bind(dofs,count,count);
+ this->energy.bind(dofs,2 * count,count);
+ this->data.setSize(2*count);
+ this->pressure.bind(this->data,0,count);
+ this->velocity.bind(this->data,count,count);
+ for(long int i = 0; i < count; i++)
+ if (i < x0 * count )
+ {
+ this->rho[i] = rhoL;
+ this->rhoVel[i] = rhoL * velL;
+ this->energy[i] = eL;
+ this->velocity[i] = velL;
+ this->pressure[i] = preL;
+ }
+ else
+ {
+ this->rho[i] = rhoR;
+ this->rhoVel[i] = rhoR * velR;
+ this->energy[i] = eR;
+ this->velocity[i] = velR;
+ this->pressure[i] = preR;
+ };
+ this->gamma = gamma;
+ cout << "dofs = " << dofs << endl;
+ getchar();
+
+
+ /*
+ const tnlString& initialConditionFile = parameters.getParameter< tnlString >( "initial-condition" );
+ if( ! dofs.load( initialConditionFile ) )
+ {
+ cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." << endl;
+ return false;
+ }
+ */
+ return true;
+}
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+ template< typename Matrix >
+bool
+eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+setupLinearSystem( const MeshType& mesh,
+ Matrix& matrix )
+{
+ const IndexType dofs = this->getDofs( mesh );
+ typedef typename Matrix::CompressedRowsLengthsVector CompressedRowsLengthsVectorType;
+ CompressedRowsLengthsVectorType rowLengths;
+ if( ! rowLengths.setSize( dofs ) )
+ return false;
+ tnlMatrixSetter< 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 MeshType& mesh,
+ DofVectorType& dofs,
+ MeshDependentDataType& meshDependentData )
+{
+ cout << endl << "Writing output at time " << time << " step " << step << "." << endl;
+ this->bindDofs( mesh, dofs );
+ tnlString fileName;
+ ofstream vysledek;
+ cout << "pressure:" << endl;
+ for (int i = 0; i<100; i++) cout << this->pressure[i] << " " ;
+ vysledek.open("pressure" + to_string(step) + ".txt");
+ for (int i = 0; i<101; i++)
+ vysledek << 0.01*i << " " << pressure[i] << endl;
+ vysledek.close();
+ cout << " " << endl;
+ cout << "velocity:" << endl;
+ for (int i = 0; i<100; i++) cout << this->velocity[i] << " " ;
+ vysledek.open("velocity" + to_string(step) + ".txt");
+ for (int i = 0; i<101; i++)
+ vysledek << 0.01*i << " " << pressure[i] << endl;
+ vysledek.close();
+ cout << "energy:" << endl;
+ for (int i = 0; i<100; i++) cout << this->energy[i] << " " ;
+ vysledek.open("energy" + to_string(step) + ".txt");
+ for (int i = 0; i<101; i++)
+ vysledek << 0.01*i << " " << energy[i] << endl;
+ vysledek.close();
+ cout << " " << endl;
+ cout << "density:" << endl;
+ for (int i = 0; i<100; i++) cout << this->rho[i] << " " ;
+ vysledek.open("density" + to_string(step) + ".txt");
+ for (int i = 0; i<101; i++)
+ vysledek << 0.01*i << " " << rho[i] << endl;
+ vysledek.close();
+ getchar();
+
+ FileNameBaseNumberEnding( "rho-", step, 5, ".tnl", fileName );
+ if( ! rho.save( fileName ) )
+ return false;
+ FileNameBaseNumberEnding( "rhoVel-", step, 5, ".tnl", fileName );
+ if( ! rhoVel.save( fileName ) )
+ return false;
+ FileNameBaseNumberEnding( "energy-", step, 5, ".tnl", fileName );
+ if( ! energy.save( fileName ) )
+ 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 MeshType& mesh,
+ DofVectorType& _u,
+ DofVectorType& _fu,
+ MeshDependentDataType& meshDependentData )
+{
+ /*
+ W[1] [0 ... count-1 ]
+ W[2] [count ... 2*count-1 ]
+ W[3] [2*count ... 3*count-1 ]
+ V this->velocity[]
+ p this->pressure[]
+ */
+ typedef typename MeshType::Cell Cell;
+ int count = mesh.template getEntitiesCount< Cell >();
+ const RealType& size = 1;//mesh.template getSpaceStepsProducts< -1, 0 >();
+ for (long int i = 1; i < count - 1; i++)
+ _fu[i] = 1.0 / (2.0*tau) * (this->rho[i-1]+this->rho[i+1]-2.0*this->rho[i])-(1.0/(2.0 * size)) * (this->rho[i+1]
+ * this->velocity[i+1] - this->rho[i-1] * this->velocity[i-1]);
+ _fu[count-1] = _fu[count-2];
+ for (long int i = 1; i < count - 1; i++)
+ _fu[count + i] =
+ 1.0 / (2.0 * tau) * (this->rhoVel[i+1] + this->rhoVel[i+1] -
+ 2.0 * this->rhoVel[i])-(1.0 / (2.0 * size)) * ((this->rhoVel[i+1]
+ * this->velocity[i + 1] + this->pressure[i + 1]) - (this->rhoVel[i-1] * this->velocity[i - 1]
+ + this->pressure[i - 1]));
+ _fu[2*count-1] = _fu[2*count-2];
+ for (long int i = 1; i < count - 1; i++)
+ _fu[i + 2 * count] = 1.0 / (2.0*tau) * (this->energy[i-1]
+ + this->energy[i+1]-2.0*this->energy[i])
+ - (1.0/(2.0 * size)) * ((this->energy[i+1]
+ + this->pressure[i + 1]) * this->velocity[i + 1] - (this->energy[i-1] + this->pressure[i -1]) * this->velocity[i - 1]);
+ _fu[3 * count-1] = _fu[3 * count-2];
+ for (long int i = 0; i <= count; i++) //pressure
+ this->pressure[i] = (this->gamma - 1 ) * ( this->energy[i] - 0.5 * this->rhoVel[i] * this->velocity[i]);
+ for (long int i = 0; i <= count ; i++) //velocity
+ this->velocity[i] = this->rhoVel[i]/this->rho[i];
+
+ /****
+ * If you use an explicit solver like tnlEulerSolver or tnlMersonSolver, you
+ * need to implement this method. Compute the right-hand side of
+ *
+ * d/dt u(x) = fu( x, u )
+ *
+ * You may use supporting mesh dependent data if you need.
+
+
+ this->bindDofs( mesh, _u );
+ tnlExplicitUpdater< Mesh, MeshFunctionType, DifferentialOperator, BoundaryCondition, RightHandSide > explicitUpdater;
+ MeshFunctionType u( mesh, _u );
+ MeshFunctionType fu( mesh, _fu );
+ explicitUpdater.template update< typename Mesh::Cell >( time,
+ mesh,
+ this->differentialOperator,
+ this->boundaryCondition,
+ this->rightHandSide,
+ u,
+ fu );
+ tnlBoundaryConditionsSetter< MeshFunctionType, BoundaryCondition > boundaryConditionsSetter;
+ boundaryConditionsSetter.template apply< typename Mesh::Cell >(
+ this->boundaryCondition,
+ time + tau,
+ u ); */
+ }
+
+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 MeshType& mesh,
+ DofVectorType& _u,
+ Matrix& matrix,
+ DofVectorType& b,
+ MeshDependentDataType& meshDependentData )
+{
+ tnlLinearSystemAssembler< Mesh,
+ MeshFunctionType,
+ DifferentialOperator,
+ BoundaryCondition,
+ RightHandSide,
+ tnlBackwardTimeDiscretisation,
+ Matrix,
+ DofVectorType > systemAssembler;
+
+ tnlMeshFunction< Mesh > u( mesh, _u );
+ systemAssembler.template assembly< typename Mesh::Cell >( time,
+ tau,
+ mesh,
+ this->differentialOperator,
+ this->boundaryCondition,
+ this->rightHandSide,
+ u,
+ matrix,
+ b );
+}
+
+#endif /* eulerPROBLEM_IMPL_H_ */
diff --git a/examples/inviscid-flow/1d/eulerRhs.h b/examples/inviscid-flow/1d/eulerRhs.h
new file mode 100644
index 0000000000000000000000000000000000000000..582dc36413ae36d78aa06a03be2ad745325263d6
--- /dev/null
+++ b/examples/inviscid-flow/1d/eulerRhs.h
@@ -0,0 +1,29 @@
+#ifndef eulerRHS_H_
+#define eulerRHS_H_
+#include<functions/tnlDomain.h>
+template< typename Mesh, typename Real >class eulerRhs
+ : public tnlDomain< Mesh::meshDimensions, MeshDomain >
+ {
+ public:
+
+ typedef Mesh MeshType;
+ typedef Real RealType;
+
+ bool setup( const tnlParameterContainer& parameters,
+ const tnlString& 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;
+ };
+};
+
+#endif /* eulerRHS_H_ */
diff --git a/examples/inviscid-flow/1d/tnl-run-euler-1d b/examples/inviscid-flow/1d/tnl-run-euler-1d
new file mode 100644
index 0000000000000000000000000000000000000000..f68b98a8406dea2f2142526283ec60248551c56d
--- /dev/null
+++ b/examples/inviscid-flow/1d/tnl-run-euler-1d
@@ -0,0 +1,19 @@
+#!/usr/bin/env bash
+
+tnl-grid-setup --dimensions 2 \
+ --origin-x 0.0 \
+ --origin-y 0.0 \
+ --proportions-x 1.0 \
+ --proportions-y 1.0 \
+ --size-x 100 \
+ --size-y 100
+
+tnl-init --test-function sin-wave \
+ --output-file init.tnl
+./euler --time-discretisation explicit \
+ --boundary-conditions-constant 0 \
+ --discrete-solver merson \
+ --snapshot-period 0.01 \
+ --final-time 1.0
+
+tnl-view --mesh mesh.tnl --input-files *tnl
diff --git a/examples/inviscid-flow/2d/CMakeLists.txt b/examples/inviscid-flow/2d/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..315b26d110246c4512d87f13bb178c9d5a07f539
--- /dev/null
+++ b/examples/inviscid-flow/2d/CMakeLists.txt
@@ -0,0 +1,21 @@
+set( tnl_heat_equation_SOURCES
+ euler.cpp
+ euler.cu )
+
+IF( BUILD_CUDA )
+ CUDA_ADD_EXECUTABLE(tnl-euler-2d${debugExt} euler.cu)
+ target_link_libraries (tnl-euler-2d${debugExt} tnl${debugExt}-${tnlVersion} ${CUSPARSE_LIBRARY} )
+ELSE( BUILD_CUDA )
+ ADD_EXECUTABLE(tnl-euler-2d${debugExt} euler.cpp)
+ target_link_libraries (tnl-euler-2d${debugExt} tnl${debugExt}-${tnlVersion} )
+ENDIF( BUILD_CUDA )
+
+
+INSTALL( TARGETS tnl-euler-2d${debugExt}
+ RUNTIME DESTINATION bin
+ PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE )
+
+INSTALL( FILES tnl-run-euler-2d
+ ${tnl_heat_equation_SOURCES}
+ DESTINATION share/tnl-${tnlVersion}/examples/inviscid-flow-2d )
+
diff --git a/examples/inviscid-flow/2d/LaxFridrichs.h b/examples/inviscid-flow/2d/LaxFridrichs.h
new file mode 100644
index 0000000000000000000000000000000000000000..db63defd3f6585702840e7bb575e225c6667c917
--- /dev/null
+++ b/examples/inviscid-flow/2d/LaxFridrichs.h
@@ -0,0 +1,143 @@
+#ifndef LaxFridrichs_H
+#define LaxFridrichs_H
+
+#include <core/vectors/tnlVector.h>
+#include <mesh/tnlGrid.h>
+
+template< typename Mesh,
+ typename Real = typename Mesh::RealType,
+ typename Index = typename Mesh::IndexType >
+class LaxFridrichs
+{
+};
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class LaxFridrichs< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index >
+{
+ public:
+ typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef tnlMeshFunction< MeshType > MeshFunctionType;
+ enum { Dimensions = MeshType::getMeshDimensions() };
+
+ static tnlString getType();
+
+ template< typename MeshFunction, typename MeshEntity >
+ __cuda_callable__
+ Real operator()( const MeshFunction& u,
+ const MeshEntity& entity,
+ const RealType& time = 0.0 ) const;
+
+ __cuda_callable__
+ template< typename MeshEntity >
+ 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 LaxFridrichs< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index >
+{
+ public:
+ typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef tnlMeshFunction< MeshType > MeshFunctionType;
+ enum { Dimensions = MeshType::getMeshDimensions() };
+
+ static tnlString getType();
+
+ template< typename MeshFunction, typename MeshEntity >
+ __cuda_callable__
+ Real operator()( const MeshFunction& u,
+ const MeshEntity& entity,
+ const RealType& time = 0.0 ) const;
+
+ __cuda_callable__
+ template< typename MeshEntity >
+ 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 LaxFridrichs< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, Index >
+{
+ public:
+ typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef tnlMeshFunction< MeshType > MeshFunctionType;
+ enum { Dimensions = MeshType::getMeshDimensions() };
+
+ static tnlString getType();
+
+ template< typename MeshFunction, typename MeshEntity >
+ __cuda_callable__
+ Real operator()( const MeshFunction& u,
+ const MeshEntity& entity,
+ const RealType& time = 0.0 ) const;
+
+ __cuda_callable__
+ template< typename MeshEntity >
+ 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;
+};
+
+
+#include "LaxFridrichs_impl.h"
+
+#endif /* LaxFridrichs_H */
diff --git a/examples/inviscid-flow/2d/LaxFridrichs_impl.h b/examples/inviscid-flow/2d/LaxFridrichs_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..7f9c73731154d797334ccb86a26e7640743a7ae7
--- /dev/null
+++ b/examples/inviscid-flow/2d/LaxFridrichs_impl.h
@@ -0,0 +1,495 @@
+#ifndef LaxFridrichs_IMPL_H
+#define LaxFridrichs_IMPL_H
+
+/****
+ * 1D problem
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString
+LaxFridrichs< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+ return tnlString( "LaxFridrichs< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename MeshFunction, typename MeshEntity >
+__cuda_callable__
+Real
+LaxFridrichs< tnlGrid< 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();
+ /*
+ //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 >();
+ typedef typename MeshType::Cell Cell;
+ const int size = mesh.template getEntitiesCount< Cell >()/5;
+ return (0.5 * this->tau)( u[ west ] - 2.0 * u[ center ] + u[ east ] )
+ - 0.5 * hxInverse * ( u[ west ] * u[ west + 3*size ] - u[ east ] * u[ east + 3*size ] );
+ */
+ /*
+ //rhoVel
+ 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 >();
+ typedef typename MeshType::Cell Cell;
+ const int size = mesh.template getEntitiesCount< Cell >()/5;
+ return (0.5 * this->tau)( u[ west ] - 2.0 * u[ center ] + u[ east ] )
+ - 0.5 * hxInverse *
+ (( u[ west ] * u[ west + 2*size ] + u[ west + 3*size ] )
+ - (u[ east ] * u[ east + 2*size ] + u[ east + 3*size ] ));
+ */
+ /*
+ //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 >();
+ typedef typename MeshType::Cell Cell;
+ const int size = mesh.template getEntitiesCount< Cell >()/5;
+ return (0.5 * this->tau)( u[ west ] - 2.0 * u[ center ] + u[ east ] )
+ - 0.5 * hxInverse *
+ (( u[ west ] + u[ west + 2*size ] ) * u[ west + size ]
+ - (u[ east ] + u[ east + 2*size ] ) * u[ east + size ] );
+ */
+ /*
+ //vel
+ const IndexType& center = entity.getIndex();
+ typedef typename MeshType::Cell Cell;
+ const int size = mesh.template getEntitiesCount< Cell >()/5;
+ return ( fu[ center - 2*size] - fu[ center - 3*size]);
+ */
+ /*
+ //pressure
+ const IndexType& center = entity.getIndex();
+ typedef typename MeshType::Cell Cell;
+ const int size = mesh.template getEntitiesCount< Cell >()/5;
+ return ( this->gamma - 1 ) * ( fu[ center - 2*size ] - 0.5 * fu[ center - 3* size ] * fu[ center - 4*size]);
+ */
+
+
+
+
+ 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 ( 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
+LaxFridrichs< tnlGrid< 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
+LaxFridrichs< tnlGrid< 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 >
+tnlString
+LaxFridrichs< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+ return tnlString( "LaxFridrichs< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename MeshFunction, typename MeshEntity >
+__cuda_callable__
+Real
+LaxFridrichs< tnlGrid< 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 >();
+ typedef typename MeshType::Cell Cell;
+ const int size = mesh.template getEntitiesCount< Cell >()/8;
+ return (0.5 * this->tau)( u[ west ] + u[ east ] + u[ south ] + u[ north ] - 4.0 * u[ center ] )
+ - 0.5 * hxInverse * ( u[ west ] * u[ west + 4*size ] - u[ east ] * u[ east + 4*size ] )
+ - 0.5 * hyInverse * ( u[ north ] * u[ north + 5*size ] - u[ south ] * u[ south + 5*size ] );
+ */
+ /*
+ //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 >();
+ typedef typename MeshType::Cell Cell;
+ const int size = mesh.template getEntitiesCount< Cell >()/8;
+ return (0.5 * this->tau)( u[ west ] + u[ east ] + u[ south ] + u[ north ] - 4.0 * u[ center ] )
+ - 0.5 * hxInverse * (( u[ west ] * u[ west + 3*size ] + u[ west + 6*size ] )
+ -( u[ east ] * u[ east + 3*size ] + u[ east + 6*size ] ))
+ - 0.5 * hyInverse * (( u[ north ] * u[ north + 4*size ] + u[ north + 6*size ] )
+ -( u[ south ] * u[ south + 4*size ] + u[ south + 6*size ] ));
+ */
+ /*
+ //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 >();
+ typedef typename MeshType::Cell Cell;
+ const int size = mesh.template getEntitiesCount< Cell >()/8;
+ return (0.5 * this->tau)( u[ west ] + u[ east ] + u[ south ] + u[ north ] - 4.0 * u[ center ] )
+ - 0.5 * hxInverse * (( u[ west ] * u[ west + 3*size ] + u[ west + 5*size ] )
+ -( u[ east ] * u[ east + 3*size ] + u[ east + 5*size ] ))
+ - 0.5 * hyInverse * (( u[ north ] * u[ north + 2*size ] + u[ north + 5*size ] )
+ -( u[ south ] * u[ south + 2*size ] + u[ south + 5*size ] ));
+ */
+ /*
+ //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 >();
+ typedef typename MeshType::Cell Cell;
+ const int size = mesh.template getEntitiesCount< Cell >()/8;
+ return (0.5 * this->tau)( u[ west ] + u[ east ] + u[ south ] + u[ north ] - 4.0 * u[ center ] )
+ - 0.5 * hxInverse * (( u[ west ] + u[ west + 4*size ] ) * u[ west + size ] )
+ -( u[ east ] + u[ east + 4*size ] ) * u[ east + size ] ))
+ - 0.5 * hyInverse * (( u[ north ] + u[ north + 4*size ] ) * u[ north + 2*size ] )
+ -( u[ south ] + u[ south + 4*size ] ) * u[ south + 2*size ] ));
+ */
+ /*
+ //velX
+ const IndexType& center = entity.getIndex();
+ typedef typename MeshType::Cell Cell;
+ const int size = mesh.template getEntitiesCount< Cell >()/8;
+ return ( fu[ center - 3*size] - fu[ center - 5*size]);
+ */
+ /*
+ //velY
+ const IndexType& center = entity.getIndex();
+ typedef typename MeshType::Cell Cell;
+ const int size = mesh.template getEntitiesCount< Cell >()/8;
+ return ( fu[ center - 3*size] - fu[ center - 5*size]);
+ */
+ /*
+ //vel
+ const IndexType& center = entity.getIndex();
+ typedef typename MeshType::Cell Cell;
+ const int size = mesh.template getEntitiesCount< Cell >()/8;
+ return sqrt(pow( fu[ center - 2*size],2)+pow( fu[ center - size],2));
+ */
+ /*
+ //pressure
+ const IndexType& center = entity.getIndex();
+ typedef typename MeshType::Cell Cell;
+ const int size = mesh.template getEntitiesCount< Cell >()/5;
+ return ( this->gamma - 1 ) * ( fu[ center - 4*size ] - 0.5 * fu[ center - 7 * size ] * pow( fu[ center - size],2 ));
+ */
+
+
+
+ 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
+LaxFridrichs< tnlGrid< 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
+LaxFridrichs< tnlGrid< 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 >
+tnlString
+LaxFridrichs< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+ return tnlString( "LaxFridrichs< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename MeshFunction, typename MeshEntity >
+__cuda_callable__
+Real
+LaxFridrichs< tnlGrid< 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
+LaxFridrichs< tnlGrid< 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
+LaxFridrichs< tnlGrid< 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 );
+}
+
+#endif /* LaxFridrichsIMPL_H */
+
diff --git a/examples/inviscid-flow/2d/euler-cuda.cu b/examples/inviscid-flow/2d/euler-cuda.cu
new file mode 100644
index 0000000000000000000000000000000000000000..4d76005cb1f70724be978ff0fa6fec63c4a8a76f
--- /dev/null
+++ b/examples/inviscid-flow/2d/euler-cuda.cu
@@ -0,0 +1 @@
+#include "euler.h"
diff --git a/examples/inviscid-flow/2d/euler.cpp b/examples/inviscid-flow/2d/euler.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4d76005cb1f70724be978ff0fa6fec63c4a8a76f
--- /dev/null
+++ b/examples/inviscid-flow/2d/euler.cpp
@@ -0,0 +1 @@
+#include "euler.h"
diff --git a/examples/inviscid-flow/2d/euler.h b/examples/inviscid-flow/2d/euler.h
new file mode 100644
index 0000000000000000000000000000000000000000..ce35381216dbaaa170a3934247cb7e4d49a216d6
--- /dev/null
+++ b/examples/inviscid-flow/2d/euler.h
@@ -0,0 +1,117 @@
+#include <tnlConfig.h>
+#include <solvers/tnlSolver.h>
+#include <solvers/tnlBuildConfigTags.h>
+#include <operators/tnlDirichletBoundaryConditions.h>
+#include <operators/tnlNeumannBoundaryConditions.h>
+#include <functions/tnlConstantFunction.h>
+#include "eulerProblem.h"
+#include "LaxFridrichs.h"
+#include "eulerRhs.h"
+#include "eulerBuildConfigTag.h"
+
+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( tnlConfigDescription & config )
+ {
+ config.addDelimiter( "euler2D settings:" );
+ config.addEntry< tnlString >( "boundary-conditions-type", "Choose the boundary conditions type.", "dirichlet");
+ config.addEntryEnum< tnlString >( "dirichlet" );
+ config.addEntryEnum< tnlString >( "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-velocityX", "This sets a value of left_x velocity." );
+ config.addEntry< double >( "left-velocityY", "This sets a value of left_y 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-velocityX", "This sets a value of right_x velocity." );
+ config.addEntry< double >( "right-velocityY", "This sets a value of right_y 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 tnlParameterContainer & parameters )
+ {
+ enum { Dimensions = MeshType::getMeshDimensions() };
+ typedef LaxFridrichs< MeshType, Real, Index > ApproximateOperator;
+ typedef eulerRhs< MeshType, Real > RightHandSide;
+ typedef tnlStaticVector < 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 tnlVector.
+ */
+ tnlString boundaryConditionsType = parameters.getParameter< tnlString >( "boundary-conditions-type" );
+ if( parameters.checkParameter( "boundary-conditions-constant" ) )
+ {
+ typedef tnlConstantFunction< Dimensions, Real > ConstantFunction;
+ if( boundaryConditionsType == "dirichlet" )
+ {
+ typedef tnlDirichletBoundaryConditions< MeshType, ConstantFunction, MeshType::getMeshDimensions(), Real, Index > BoundaryConditions;
+ typedef eulerProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
+ SolverStarter solverStarter;
+ return solverStarter.template run< Problem >( parameters );
+ }
+ typedef tnlNeumannBoundaryConditions< MeshType, ConstantFunction, Real, Index > BoundaryConditions;
+ typedef eulerProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
+ SolverStarter solverStarter;
+ return solverStarter.template run< Problem >( parameters );
+ }
+ typedef tnlMeshFunction< MeshType > MeshFunction;
+ if( boundaryConditionsType == "dirichlet" )
+ {
+ typedef tnlDirichletBoundaryConditions< MeshType, MeshFunction, MeshType::getMeshDimensions(), Real, Index > BoundaryConditions;
+ typedef eulerProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
+ SolverStarter solverStarter;
+ return solverStarter.template run< Problem >( parameters );
+ }
+ typedef tnlNeumannBoundaryConditions< 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[] )
+{
+ tnlSolver< eulerSetter, eulerConfig, BuildConfig > solver;
+ if( ! solver. run( argc, argv ) )
+ return EXIT_FAILURE;
+ return EXIT_SUCCESS;
+}
+
diff --git a/examples/inviscid-flow/2d/eulerBuildConfigTag.h b/examples/inviscid-flow/2d/eulerBuildConfigTag.h
new file mode 100644
index 0000000000000000000000000000000000000000..d242a88650f98f2cac9e0ea68473a10e7c29469f
--- /dev/null
+++ b/examples/inviscid-flow/2d/eulerBuildConfigTag.h
@@ -0,0 +1,43 @@
+#ifndef eulerBUILDCONFIGTAG_H_
+#define eulerBUILDCONFIGTAG_H_
+
+#include <solvers/tnlBuildConfigTags.h>
+
+class eulerBuildConfigTag{};
+
+/****
+ * Turn off support for float and long double.
+ */
+template<> struct tnlConfigTagReal< eulerBuildConfigTag, float > { enum { enabled = false }; };
+template<> struct tnlConfigTagReal< eulerBuildConfigTag, long double > { enum { enabled = false }; };
+
+/****
+ * Turn off support for short int and long int indexing.
+ */
+template<> struct tnlConfigTagIndex< eulerBuildConfigTag, short int >{ enum { enabled = false }; };
+template<> struct tnlConfigTagIndex< eulerBuildConfigTag, long int >{ enum { enabled = false }; };
+
+/****
+ * Use of tnlGrid is enabled for allowed dimensions and Real, Device and Index types.
+ */
+
+template< int Dimensions, typename Real, typename Device, typename Index >
+ struct tnlConfigTagMesh< eulerBuildConfigTag, tnlGrid< Dimensions, Real, Device, Index > >
+ { enum { enabled = tnlConfigTagDimensions< eulerBuildConfigTag, Dimensions >::enabled &&
+ tnlConfigTagReal< eulerBuildConfigTag, Real >::enabled &&
+ tnlConfigTagDevice< eulerBuildConfigTag, Device >::enabled &&
+ tnlConfigTagIndex< eulerBuildConfigTag, Index >::enabled }; };
+
+/****
+ * Please, chose your preferred time discretisation here.
+ */
+template<> struct tnlConfigTagTimeDiscretisation< eulerBuildConfigTag, tnlExplicitTimeDiscretisationTag >{ enum { enabled = true }; };
+template<> struct tnlConfigTagTimeDiscretisation< eulerBuildConfigTag, tnlSemiImplicitTimeDiscretisationTag >{ enum { enabled = true }; };
+template<> struct tnlConfigTagTimeDiscretisation< eulerBuildConfigTag, tnlImplicitTimeDiscretisationTag >{ enum { enabled = true }; };
+
+/****
+ * Only the Runge-Kutta-Merson solver is enabled by default.
+ */
+template<> struct tnlConfigTagExplicitSolver< eulerBuildConfigTag, tnlExplicitEulerSolverTag >{ enum { enabled = false }; };
+
+#endif /* eulerBUILDCONFIGTAG_H_ */
diff --git a/examples/inviscid-flow/2d/eulerProblem.h b/examples/inviscid-flow/2d/eulerProblem.h
new file mode 100644
index 0000000000000000000000000000000000000000..6a514c1da74cf0149feea852fdc227b08b1355cc
--- /dev/null
+++ b/examples/inviscid-flow/2d/eulerProblem.h
@@ -0,0 +1,93 @@
+#ifndef eulerPROBLEM_H_
+#define eulerPROBLEM_H_
+
+#include <problems/tnlPDEProblem.h>
+#include <functions/tnlMeshFunction.h>
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+class eulerProblem:
+ public tnlPDEProblem< 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 tnlMeshFunction< Mesh > MeshFunctionType;
+ typedef tnlPDEProblem< Mesh, RealType, DeviceType, IndexType > BaseType;
+
+ using typename BaseType::MeshType;
+ using typename BaseType::DofVectorType;
+ using typename BaseType::MeshDependentDataType;
+ tnlSharedVector< RealType, DeviceType, IndexType > rho;
+ tnlSharedVector< RealType, DeviceType, IndexType > rhoVelX;
+ tnlSharedVector< RealType, DeviceType, IndexType > rhoVelY;
+ tnlSharedVector< RealType, DeviceType, IndexType > energy;
+ tnlVector< RealType, DeviceType, IndexType > data;
+ tnlSharedVector< RealType, DeviceType, IndexType > pressure;
+ tnlSharedVector< RealType, DeviceType, IndexType > velocity;
+ tnlSharedVector< RealType, DeviceType, IndexType > velocityX;
+ tnlSharedVector< RealType, DeviceType, IndexType > velocityY;
+ double gamma;
+
+ static tnlString getTypeStatic();
+
+ tnlString getPrologHeader() const;
+
+ void writeProlog( tnlLogger& logger,
+ const tnlParameterContainer& parameters ) const;
+
+ bool setup( const tnlParameterContainer& parameters );
+
+ bool setInitialCondition( const tnlParameterContainer& parameters,
+ const MeshType& mesh,
+ DofVectorType& dofs,
+ MeshDependentDataType& meshDependentData );
+
+ template< typename Matrix >
+ bool setupLinearSystem( const MeshType& mesh,
+ Matrix& matrix );
+
+ bool makeSnapshot( const RealType& time,
+ const IndexType& step,
+ const MeshType& mesh,
+ DofVectorType& dofs,
+ MeshDependentDataType& meshDependentData );
+
+ IndexType getDofs( const MeshType& mesh ) const;
+
+ void bindDofs( const MeshType& mesh,
+ DofVectorType& dofs );
+
+ void getExplicitRHS( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ DofVectorType& _u,
+ DofVectorType& _fu,
+ MeshDependentDataType& meshDependentData );
+
+ template< typename Matrix >
+ void assemblyLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ DofVectorType& dofs,
+ Matrix& matrix,
+ DofVectorType& rightHandSide,
+ MeshDependentDataType& meshDependentData );
+
+ protected:
+
+ DifferentialOperator differentialOperator;
+ BoundaryCondition boundaryCondition;
+ RightHandSide rightHandSide;
+};
+
+#include "eulerProblem_impl.h"
+
+#endif /* eulerPROBLEM_H_ */
diff --git a/examples/inviscid-flow/2d/eulerProblem_impl.h b/examples/inviscid-flow/2d/eulerProblem_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..aedadddb10efa559e2ec138974cd36f2b9368833
--- /dev/null
+++ b/examples/inviscid-flow/2d/eulerProblem_impl.h
@@ -0,0 +1,355 @@
+#ifndef eulerPROBLEM_IMPL_H_
+#define eulerPROBLEM_IMPL_H_
+
+#include <core/mfilename.h>
+#include <matrices/tnlMatrixSetter.h>
+#include <solvers/pde/tnlExplicitUpdater.h>
+#include <solvers/pde/tnlLinearSystemAssembler.h>
+#include <solvers/pde/tnlBackwardTimeDiscretisation.h>
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+tnlString
+eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+getTypeStatic()
+{
+ return tnlString( "eulerProblem< " ) + Mesh :: getTypeStatic() + " >";
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+tnlString
+eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+getPrologHeader() const
+{
+ return tnlString( "euler2D" );
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+void
+eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+writeProlog( tnlLogger& logger, const tnlParameterContainer& 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 tnlParameterContainer& parameters )
+{
+ if( ! this->boundaryCondition.setup( parameters, "boundary-conditions-" ) ||
+ ! this->rightHandSide.setup( parameters, "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 MeshType& mesh ) const
+{
+ /****
+ * Return number of DOFs (degrees of freedom) i.e. number
+ * of unknowns to be resolved by the main solver.
+ */
+ return 4*mesh.template getEntitiesCount< typename MeshType::Cell >();
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+void
+eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+bindDofs( const MeshType& mesh,
+ DofVectorType& dofVector )
+{
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+bool
+eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+setInitialCondition( const tnlParameterContainer& parameters,
+ const MeshType& mesh,
+ DofVectorType& dofs,
+ MeshDependentDataType& meshDependentData )
+{
+ typedef typename MeshType::Cell Cell;
+ double gamma = parameters.getParameter< double >( "gamma" );
+ double rhoL = parameters.getParameter< double >( "left-density" );
+ double velLX = parameters.getParameter< double >( "left-velocityX" );
+ double velLY = parameters.getParameter< double >( "left-velocityY" );
+ double preL = parameters.getParameter< double >( "left-pressure" );
+ double eL = ( preL / (gamma - 1) ) + 0.5 * rhoL * pow(velLX,2)+pow(velLY,2);
+ double rhoR = parameters.getParameter< double >( "right-density" );
+ double velRX = parameters.getParameter< double >( "right-velocityX" );
+ double velRY = parameters.getParameter< double >( "right-velocityY" );
+ double preR = parameters.getParameter< double >( "right-pressure" );
+ double eR = ( preR / (gamma - 1) ) + 0.5 * rhoR * pow(velRX,2)+pow(velRY,2);
+ double x0 = parameters.getParameter< double >( "riemann-border" );
+ int size = mesh.template getEntitiesCount< Cell >();
+ int size2 = pow(size,2);
+ this->rho.bind(dofs,0,size2);
+ this->rhoVelX.bind(dofs,size2,size2);
+ this->rhoVelY.bind(dofs,2*size2,size2);
+ this->energy.bind(dofs,3*size2,size2);
+ this->data.setSize(4*size2);
+ this->pressure.bind(this->data,0,size2);
+ this->velocity.bind(this->data,size2,size2);
+ this->velocityX.bind(this->data,2*size2,size2);
+ this->velocityY.bind(this->data,3*size2,size2);
+ for(long int j = 0; j < size; j++)
+ for(long int i = 0; i < size; i++)
+ if ((i < x0 * size)&&(j < x0 * size) )
+ {
+ this->rho[j*size+i] = rhoL;
+ this->rhoVelX[j*size+i] = rhoL * velLX;
+ this->rhoVelY[j*size+i] = rhoL * velLY;
+ this->energy[j*size+i] = eL;
+ this->velocity[j*size+i] = sqrt(pow(velLX,2)+pow(velLY,2));
+ this->velocityX[j*size+i] = velLX;
+ this->velocityY[j*size+i] = velLY;
+ this->pressure[j*size+i] = preL;
+ }
+ else
+ {
+ this->rho[j*size+i] = rhoR;
+ this->rhoVelX[j*size+i] = rhoR * velRX;
+ this->rhoVelY[j*size+i] = rhoR * velRY;
+ this->energy[j*size+i] = eR;
+ this->velocity[j*size+i] = sqrt(pow(velRX,2)+pow(velRY,2));
+ this->velocityX[j*size+i] = velRX;
+ this->velocityY[j*size+i] = velRY;
+ this->pressure[j*size+i] = preR;
+ };
+ this->gamma = gamma;
+ return true;
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+ template< typename Matrix >
+bool
+eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+setupLinearSystem( const MeshType& mesh,
+ Matrix& matrix )
+{
+ const IndexType dofs = this->getDofs( mesh );
+ typedef typename Matrix::CompressedRowsLengthsVector CompressedRowsLengthsVectorType;
+ CompressedRowsLengthsVectorType rowLengths;
+ if( ! rowLengths.setSize( dofs ) )
+ return false;
+ tnlMatrixSetter< 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 MeshType& mesh,
+ DofVectorType& dofs,
+ MeshDependentDataType& meshDependentData )
+{
+ cout << endl << "Writing output at time " << time << " step " << step << "." << endl;
+ this->bindDofs( mesh, dofs );
+ tnlString fileName;
+ FileNameBaseNumberEnding( "rho-", step, 5, ".tnl", fileName );
+ if( ! this->rho.save( fileName ) )
+ return false;
+ FileNameBaseNumberEnding( "rhoVelX-", step, 5, ".tnl", fileName );
+ if( ! this->rhoVelX.save( fileName ) )
+ return false;
+ FileNameBaseNumberEnding( "rhoVelY-", step, 5, ".tnl", fileName );
+ if( ! this->rhoVelY.save( fileName ) )
+ return false;
+ FileNameBaseNumberEnding( "energy-", step, 5, ".tnl", fileName );
+ if( ! this->energy.save( fileName ) )
+ return false;
+ FileNameBaseNumberEnding( "velocityX-", step, 5, ".tnl", fileName );
+ if( ! this->velocityX.save( fileName ) )
+ return false;
+ FileNameBaseNumberEnding( "velocityY-", step, 5, ".tnl", fileName );
+ if( ! this->velocityY.save( fileName ) )
+ return false;
+ FileNameBaseNumberEnding( "velocity-", step, 5, ".tnl", fileName );
+ if( ! this->velocity.save( fileName ) )
+ return false;
+ FileNameBaseNumberEnding( "pressure-", step, 5, ".tnl", fileName );
+ if( ! this->pressure.save( fileName ) )
+ 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 MeshType& mesh,
+ DofVectorType& _u,
+ DofVectorType& _fu,
+ MeshDependentDataType& meshDependentData )
+{
+ typedef typename MeshType::Cell Cell;
+ const RealType& cellSize = 1;// mesh.template getSpaceStepsProducts< -1, 0 >();
+ int size = mesh.template getEntitiesCount< Cell >()/4;
+// prepsat na SWEN
+
+ for (long int j = 1; j < size - 1; j++)
+ {
+ for (long int i = 1; i < size - 1; i++)
+ {
+ _fu[j*size+i] = 1.0 / (4.0*tau) *
+ (this->rho[j*size+i-1]+this->rho[j*size+i+1]+this->rho[(j-1)*size+i]+this->rho[(j+1)*size+i]-4.0*this->rho[j*size+i])
+ -(1.0/(2.0*cellSize))*(this->rho[j*size+i+1]*this->velocityX[j*size+i+1]-this->rho[j*size+i-1]*this->velocityX[j*size+i-1])
+ -(1.0/(2.0*cellSize))*(this->rho[(j+1)*size+i]*this->velocityY[(j+1)*size+i]-this->rho[(j-1)*size+i]*this->velocityY[(j-1)*size+i]);
+ _fu[pow(size,2)-size+i]=_fu[pow(size,2)-2*size+i];
+ };
+ _fu[j*size-1] = _fu[j*size-1];
+ };
+ for (long int j = 1; j < size - 1; j++)
+ {
+ for (long int i = 1; i < size - 1; i++)
+ {
+ _fu[pow(size,2) + j*size+i] = 1.0 / (4.0 * tau) *
+ (this->rhoVelX[j*size+i+1]+this->rhoVelX[j*size+i-1]+this->rhoVelX[(j-1)*size+i]+this->rhoVelX[(j+1)*size+i]-4.0*this->rhoVelX[j*size+i])
+ -(1.0/(2.0*cellSize))*((this->rhoVelX[j*size+i+1]*this->velocityX[j*size+i+1]+this->pressure[j*size+i+1])-(this->rhoVelX[j*size+i-1]*this->velocityX[j*size+i-1]+this->pressure[j*size+i-1]))
+ -(1.0/(2.0*cellSize))*((this->rhoVelX[(j+1)*size+i]*this->velocityY[(j+1)*size+i])-(this->rhoVelX[(j-1)*size+i]*this->velocityY[(j-1)*size+i]));
+ _fu[2*pow(size,2)-size+i]=_fu[2*pow(size,2)-2*size+i];
+ };
+ _fu[pow(size,2)+j*size] = _fu[pow(size,2)+j*size-1];
+ };
+ for (long int j = 1; j < size - 1; j++)
+ {
+ for (long int i = 1; i < size - 1; i++)
+ {
+ _fu[2*pow(size,2) + j*size+i] = 1.0 / (4.0 * tau) *
+ (this->rhoVelY[j*size+i+1]+this->rhoVelY[j*size+i-1]+this->rhoVelY[(j-1)*size+i]+this->rhoVelY[(j+1)*size+i]-4.0*this->rhoVelY[j*size+i])
+ -(1.0/(2.0*cellSize))*((this->rhoVelY[(j+1)*size+i]*this->velocityY[(j+1)*size+i]+this->pressure[(j+1)*size+i])-(this->rhoVelY[(j-1)*size+i]*this->velocityY[(j-1)*size+i]+this->pressure[(j-1)*size+i]))
+ -(1.0/(2.0*cellSize))*((this->rhoVelY[j*size+i+1]*this->velocityX[j*size+i+1])-(this->rhoVelY[j*size+i-1]*this->velocityX[j*size+i-1]));
+ _fu[3*pow(size,2)-size+i]=_fu[3*pow(size,2)-2*size+i];
+ };
+ _fu[2*pow(size,2)+j*size] = _fu[2*pow(size,2)+j*size-1];
+ };
+ for (long int j = 1; j < size - 1; j++)
+ {
+ for (long int i = 1; i < size - 1; i++)
+ {
+ _fu[3*pow(size,2) + j*size+i] = 1.0 / (4.0 * tau) *
+ (this->energy[j*size+i+1]+this->energy[j*size+i-1]+this->energy[(j-1)*size+i]+this->energy[(j+1)*size+i]-4.0*this->energy[j*size+i])
+ -(1.0/(2.0*cellSize))*((this->energy[j*size+i+1]+this->pressure[j*size+i+1])*this->velocityX[j*size+i+1]-(this->energy[j*size+i-1]+this->pressure[j*size+i-1])*this->velocityX[j*size+i-1])
+ -(1.0/(2.0*cellSize))*((this->energy[(j+1)*size+i]+this->pressure[(j+1)*size+i])*this->velocityY[(j+1)*size+i]-(this->energy[(j-1)*size+i]+this->pressure[(j-1)*size+i])*this->velocityY[(j-1)*size+i]);
+ _fu[4*pow(size,2)-size+i] = _fu[4*pow(size,2)-2*size+i];
+ };
+ _fu[3*pow(size,2)+j*size] = _fu[3*pow(size,2)+j*size-1];
+ };
+
+ for (long int j = 1; j < size; j++) //pressure
+ for (long int i = 1; i < size; i++)
+ this->pressure[j*size+i] = (this->gamma - 1 ) * ( this->energy[j*size+i] - 0.5 * this->rho[j*size+i] * (pow(this->velocityX[j*size+i],2) + pow(this->velocityY[j*size+i],2)));
+ for (long int j = 1; j < size; j++) //velocityX
+ for (long int i = 1; i < size; i++)
+ this->velocityX[j*size+i] = this->rhoVelX[j*size+i]/this->rho[j*size+i];
+ for (long int j = 1; j < size; j++) //velocityY
+ for (long int i = 1; i < size; i++)
+ this->velocityY[j*size+i] = this->rhoVelY[j*size+i]/this->rho[j*size+i];
+ for (long int j = 1; j < size; j++) //velocity
+ for (long int i = 1; i < size; i++)
+ this->velocity[j*size+i] =sqrt(pow(velocityX[j*size+i],2)+pow(velocityY[j*size+i],2));
+
+/*
+
+ this->bindDofs( mesh, _u );
+ tnlExplicitUpdater< Mesh, MeshFunctionType, DifferentialOperator, BoundaryCondition, RightHandSide > explicitUpdater;
+ MeshFunctionType u( mesh, _u );
+ MeshFunctionType fu( mesh, _fu );
+ explicitUpdater.template update< typename Mesh::Cell >( time,
+ mesh,
+ this->differentialOperator,
+ this->boundaryCondition,
+ this->rightHandSide,
+ u,
+ fu );
+ tnlBoundaryConditionsSetter< MeshFunctionType, BoundaryCondition > boundaryConditionsSetter;
+ boundaryConditionsSetter.template apply< typename Mesh::Cell >(
+ this->boundaryCondition,
+ time + tau,
+ u );
+*/
+ }
+
+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 MeshType& mesh,
+ DofVectorType& _u,
+ Matrix& matrix,
+ DofVectorType& b,
+ MeshDependentDataType& meshDependentData )
+{
+ tnlLinearSystemAssembler< Mesh,
+ MeshFunctionType,
+ DifferentialOperator,
+ BoundaryCondition,
+ RightHandSide,
+ tnlBackwardTimeDiscretisation,
+ Matrix,
+ DofVectorType > systemAssembler;
+
+ tnlMeshFunction< Mesh > u( mesh, _u );
+ systemAssembler.template assembly< typename Mesh::Cell >( time,
+ tau,
+ mesh,
+ this->differentialOperator,
+ this->boundaryCondition,
+ this->rightHandSide,
+ u,
+ matrix,
+ b );
+}
+
+#endif /* eulerPROBLEM_IMPL_H_ */
diff --git a/examples/inviscid-flow/2d/eulerRhs.h b/examples/inviscid-flow/2d/eulerRhs.h
new file mode 100644
index 0000000000000000000000000000000000000000..582dc36413ae36d78aa06a03be2ad745325263d6
--- /dev/null
+++ b/examples/inviscid-flow/2d/eulerRhs.h
@@ -0,0 +1,29 @@
+#ifndef eulerRHS_H_
+#define eulerRHS_H_
+#include<functions/tnlDomain.h>
+template< typename Mesh, typename Real >class eulerRhs
+ : public tnlDomain< Mesh::meshDimensions, MeshDomain >
+ {
+ public:
+
+ typedef Mesh MeshType;
+ typedef Real RealType;
+
+ bool setup( const tnlParameterContainer& parameters,
+ const tnlString& 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;
+ };
+};
+
+#endif /* eulerRHS_H_ */
diff --git a/examples/inviscid-flow/2d/run-euler b/examples/inviscid-flow/2d/run-euler
new file mode 100644
index 0000000000000000000000000000000000000000..f68b98a8406dea2f2142526283ec60248551c56d
--- /dev/null
+++ b/examples/inviscid-flow/2d/run-euler
@@ -0,0 +1,19 @@
+#!/usr/bin/env bash
+
+tnl-grid-setup --dimensions 2 \
+ --origin-x 0.0 \
+ --origin-y 0.0 \
+ --proportions-x 1.0 \
+ --proportions-y 1.0 \
+ --size-x 100 \
+ --size-y 100
+
+tnl-init --test-function sin-wave \
+ --output-file init.tnl
+./euler --time-discretisation explicit \
+ --boundary-conditions-constant 0 \
+ --discrete-solver merson \
+ --snapshot-period 0.01 \
+ --final-time 1.0
+
+tnl-view --mesh mesh.tnl --input-files *tnl
diff --git a/examples/inviscid-flow/CMakeLists.txt b/examples/inviscid-flow/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..10af3c9c773bc36bca44191359b4a4e244d1faa4
--- /dev/null
+++ b/examples/inviscid-flow/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_subdirectory( 1d )
+add_subdirectory( 2d )