diff --git a/examples/flow/1DBoundaryConditions.h b/examples/flow/1DBoundaryConditions.h
new file mode 100644
index 0000000000000000000000000000000000000000..0a07664dcbf5b3bd8b034641fa8436200a07085d
--- /dev/null
+++ b/examples/flow/1DBoundaryConditions.h
@@ -0,0 +1,122 @@
+#include <TNL/Functions/FunctionAdapter.h>
+
+#include "1DDensityBoundaryCondition.h"
+#include "1DMomentumXBoundaryCondition.h"
+#include "1DMomentumYBoundaryCondition.h"
+#include "1DEnergyBoundaryCondition.h"
+
+namespace TNL {
+
+template< typename Mesh,
+ typename Function,
+ typename Real = typename Mesh::RealType,
+ typename Index = typename Mesh::IndexType >
+class BoundaryConditions
+{
+ public:
+ typedef Mesh MeshType;
+ typedef Real RealType;
+ typedef Index IndexType;
+ typedef Function FunctionType;
+ typedef Functions::MeshFunction< Mesh > MeshFunctionType;
+ typedef typename Mesh::DeviceType DeviceType;
+
+ typedef TNL::Operators::DensityBoundaryConditions< MeshType, FunctionType, RealType, IndexType > DensityBoundaryConditionsType;
+ typedef TNL::Operators::MomentumXBoundaryConditions< MeshType, FunctionType, RealType, IndexType > MomentumXBoundaryConditionsType;
+ typedef TNL::Operators::MomentumYBoundaryConditions< MeshType, FunctionType, RealType, IndexType > MomentumYBoundaryConditionsType;
+ typedef TNL::Operators::EnergyBoundaryConditions< MeshType, FunctionType, RealType, IndexType > EnergyBoundaryConditionsType;
+ typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
+
+ typedef SharedPointer< DensityBoundaryConditionsType > DensityBoundaryConditionsTypePointer;
+ typedef SharedPointer< MomentumXBoundaryConditionsType > MomentumXBoundaryConditionsTypePointer;
+ typedef SharedPointer< MomentumYBoundaryConditionsType > MomentumYBoundaryConditionsTypePointer;
+ typedef SharedPointer< EnergyBoundaryConditionsType > EnergyBoundaryConditionsTypePointer;
+ typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
+ typedef SharedPointer< MeshType > MeshPointer;
+ typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
+
+ static void configSetup( Config::ConfigDescription& config,
+ const String& prefix = "" )
+ {
+ }
+
+ bool setup( const MeshPointer& meshPointer,
+ const Config::ParameterContainer& parameters,
+ const String& prefix = "" )
+ {
+ this->densityBoundaryConditionsPointer->setup( meshPointer, parameters, prefix);
+ this->momentumXBoundaryConditionsPointer->setup( meshPointer, parameters, prefix);
+ this->momentumYBoundaryConditionsPointer->setup( meshPointer, parameters, prefix);
+ this->energyBoundaryConditionsPointer->setup( meshPointer, parameters, prefix);
+ return true;
+ }
+
+ void setCompressibleConservativeVariables(const CompressibleConservativeVariablesPointer& compressibleConservativeVariables)
+ {
+ this->densityBoundaryConditionsPointer->setCompressibleConservativeVariables(compressibleConservativeVariables);
+ this->momentumXBoundaryConditionsPointer->setCompressibleConservativeVariables(compressibleConservativeVariables);
+ this->momentumYBoundaryConditionsPointer->setCompressibleConservativeVariables(compressibleConservativeVariables);
+ this->energyBoundaryConditionsPointer->setCompressibleConservativeVariables(compressibleConservativeVariables);
+ }
+
+ void setTimestep(const RealType timestep)
+ {
+ this->densityBoundaryConditionsPointer->setTimestep(timestep);
+ this->momentumXBoundaryConditionsPointer->setTimestep(timestep);
+ this->momentumYBoundaryConditionsPointer->setTimestep(timestep);
+ this->energyBoundaryConditionsPointer->setTimestep(timestep);
+ }
+
+ void setGamma(const RealType gamma)
+ {
+ this->densityBoundaryConditionsPointer->setGamma(gamma);
+ this->momentumXBoundaryConditionsPointer->setGamma(gamma);
+ this->momentumYBoundaryConditionsPointer->setGamma(gamma);
+ this->energyBoundaryConditionsPointer->setGamma(gamma);
+ }
+
+ void setPressure(const MeshFunctionPointer& pressure)
+ {
+ this->densityBoundaryConditionsPointer->setPressure(pressure);
+ this->momentumXBoundaryConditionsPointer->setPressure(pressure);
+ this->momentumYBoundaryConditionsPointer->setPressure(pressure);
+ this->energyBoundaryConditionsPointer->setPressure(pressure);
+ }
+
+ void setCavitySpeed(const RealType cavitySpeed)
+ {
+ this->momentumXBoundaryConditionsPointer->setCavitySpeed(cavitySpeed);
+ this->momentumYBoundaryConditionsPointer->setCavitySpeed(cavitySpeed);
+ this->energyBoundaryConditionsPointer->setCavitySpeed(cavitySpeed);
+ }
+
+ DensityBoundaryConditionsTypePointer& getDensityBoundaryCondition()
+ {
+ return this->densityBoundaryConditionsPointer;
+ }
+
+ MomentumXBoundaryConditionsTypePointer& getMomentumXBoundaryCondition()
+ {
+ return this->momentumXBoundaryConditionsPointer;
+ }
+
+ MomentumYBoundaryConditionsTypePointer& getMomentumYBoundaryCondition()
+ {
+ return this->momentumYBoundaryConditionsPointer;
+ }
+
+ EnergyBoundaryConditionsTypePointer& getEnergyBoundaryCondition()
+ {
+ return this->energyBoundaryConditionsPointer;
+ }
+
+
+ protected:
+ DensityBoundaryConditionsTypePointer densityBoundaryConditionsPointer;
+ MomentumXBoundaryConditionsTypePointer momentumXBoundaryConditionsPointer;
+ MomentumYBoundaryConditionsTypePointer momentumYBoundaryConditionsPointer;
+ EnergyBoundaryConditionsTypePointer energyBoundaryConditionsPointer;
+
+};
+
+} //namespace TNL
diff --git a/examples/flow/1DDensityBoundaryCondition.h b/examples/flow/1DDensityBoundaryCondition.h
new file mode 100644
index 0000000000000000000000000000000000000000..b45afdd6d34b1b0d096cf13ac7923251cb413151
--- /dev/null
+++ b/examples/flow/1DDensityBoundaryCondition.h
@@ -0,0 +1,542 @@
+/***************************************************************************
+ IdentityOperator.h - description
+ -------------------
+ begin : Nov 17, 2014
+ copyright : (C) 2014 by Tomas Oberhuber
+ email : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+
+#pragma once
+
+#include <TNL/Functions/FunctionAdapter.h>
+
+namespace TNL {
+namespace Operators {
+
+template< typename Mesh,
+ typename Function,
+ typename Real = typename Mesh::RealType,
+ typename Index = typename Mesh::GlobalIndexType >
+class DensityBoundaryConditions
+{
+
+};
+
+/****
+ * Base
+ */
+template< typename Function >
+class DensityBoundaryConditionsBase
+{
+ public:
+
+ typedef Function FunctionType;
+
+ static void configSetup( const Config::ConfigDescription& config,
+ const String& prefix = "" )
+ {
+ Function::configSetup( config, prefix );
+ }
+
+ template< typename MeshPointer >
+ bool setup( const MeshPointer& meshPointer,
+ const Config::ParameterContainer& parameters,
+ const String& prefix = "" )
+ {
+ return Functions::FunctionAdapter< typename MeshPointer::ObjectType, FunctionType >::setup( this->function, meshPointer, parameters, prefix );
+ }
+
+ static void configSetup( Config::ConfigDescription& config,
+ const String& prefix = "" )
+ {
+ Function::configSetup( config, prefix );
+ };
+
+ bool setup( const Config::ParameterContainer& parameters,
+ const String& prefix = "" )
+ {
+ return this->function.setup( parameters, prefix );
+ };
+
+ void setFunction( const FunctionType& function )
+ {
+ this->function = function;
+ };
+
+ FunctionType& getFunction()
+ {
+ return this->function;
+ }
+
+ const FunctionType& getFunction() const
+ {
+ return this->function;
+ };
+
+ protected:
+
+ FunctionType function;
+
+};
+
+/****
+ * 1D grid
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Function,
+ typename Real,
+ typename Index >
+class DensityBoundaryConditions< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public DensityBoundaryConditionsBase< Function >,
+ public Operator< Meshes::Grid< 1, MeshReal, Device, MeshIndex >,
+ Functions::MeshBoundaryDomain,
+ 1, 1,
+ Real,
+ Index >
+{
+ public:
+
+ typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ typedef Function FunctionType;
+ typedef Functions::MeshFunction< MeshType > MeshFunctionType;
+ typedef Containers::Vector< RealType, DeviceType, IndexType> DofVectorType;
+ typedef Containers::StaticVector< 1, RealType > PointType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef DensityBoundaryConditions< MeshType, Function, Real, Index > ThisType;
+ typedef DensityBoundaryConditionsBase< Function > BaseType;
+ typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
+ typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
+ typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
+
+
+ template< typename EntityType,
+ typename MeshFunction >
+ __cuda_callable__
+ const RealType operator()( const MeshFunction& u,
+ const EntityType& entity,
+ const RealType& time = 0 ) const
+ {
+ const MeshType& mesh = entity.getMesh();
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ if( entity.getCoordinates().x() == 0 )
+ return u[ neighborEntities.template getEntityIndex< 0 >() ];
+ else
+ return u[ neighborEntities.template getEntityIndex< -1 >() ];
+
+ }
+
+
+ template< typename EntityType >
+ __cuda_callable__
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const EntityType& entity ) const
+ {
+ return 2;
+ }
+
+ template< typename PreimageFunction,
+ typename MeshEntity,
+ typename Matrix,
+ typename Vector >
+ __cuda_callable__
+ void setMatrixElements( const PreimageFunction& u,
+ const MeshEntity& entity,
+ const RealType& time,
+ const RealType& tau,
+ Matrix& matrix,
+ Vector& b ) const
+ {
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
+ if( entity.getCoordinates().x() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 1 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ else
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< -1 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ }
+
+ void setTimestep(const RealType timestep )
+ {
+ this->timestep = timestep;
+ }
+
+ void setGamma(const RealType gamma )
+ {
+ this->gamma = gamma;
+ }
+
+ void setCompressibleConservativeVariables(const CompressibleConservativeVariablesPointer& compressibleConservativeVariables)
+ {
+ this->compressibleConservativeVariables = compressibleConservativeVariables;
+ }
+
+ void setPressure(const MeshFunctionPointer& pressure)
+ {
+ this->pressure = pressure;
+ }
+
+ private:
+ CompressibleConservativeVariablesPointer compressibleConservativeVariables;
+ RealType timestep;
+ RealType gamma;
+ MeshFunctionPointer pressure;
+
+};
+
+/****
+ * 2D grid
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Function,
+ typename Real,
+ typename Index >
+class DensityBoundaryConditions< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public DensityBoundaryConditionsBase< Function >,
+ public Operator< Meshes::Grid< 2, MeshReal, Device, MeshIndex >,
+ Functions::MeshBoundaryDomain,
+ 2, 2,
+ Real,
+ Index >
+
+{
+ public:
+
+ typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ typedef Function FunctionType;
+ typedef Functions::MeshFunction< MeshType > MeshFunctionType;
+ typedef Containers::Vector< RealType, DeviceType, IndexType> DofVectorType;
+ typedef Containers::StaticVector< 2, RealType > PointType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef DensityBoundaryConditions< MeshType, Function, Real, Index > ThisType;
+ typedef DensityBoundaryConditionsBase< Function > BaseType;
+ typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
+ typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
+ typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
+
+
+ template< typename EntityType,
+ typename MeshFunction >
+ __cuda_callable__
+ const RealType operator()( const MeshFunction& u,
+ const EntityType& entity,
+ const RealType& time = 0 ) const
+ {
+ const MeshType& mesh = entity.getMesh();
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ if( entity.getCoordinates().x() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ }
+ // The following line is commented to avoid compiler warning
+ //if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ }
+ }
+
+ template< typename EntityType >
+ __cuda_callable__
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const EntityType& entity ) const
+ {
+ return 2;
+ }
+
+ template< typename PreimageFunction,
+ typename MeshEntity,
+ typename Matrix,
+ typename Vector >
+ __cuda_callable__
+ void setMatrixElements( const PreimageFunction& u,
+ const MeshEntity& entity,
+ const RealType& time,
+ const RealType& tau,
+ Matrix& matrix,
+ Vector& b ) const
+ {
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
+ if( entity.getCoordinates().x() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 1, 0 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< -1, 0 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 0, 1 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< 0, -1 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ }
+
+ void setTimestep(const RealType timestep )
+ {
+ this->timestep = timestep;
+ }
+
+ void setGamma(const RealType gamma )
+ {
+ this->gamma = gamma;
+ }
+
+ void setCompressibleConservativeVariables(const CompressibleConservativeVariablesPointer& compressibleConservativeVariables)
+ {
+ this->compressibleConservativeVariables = compressibleConservativeVariables;
+ }
+
+ void setPressure(const MeshFunctionPointer& pressure)
+ {
+ this->pressure = pressure;
+ }
+
+ private:
+ CompressibleConservativeVariablesPointer compressibleConservativeVariables;
+ RealType timestep;
+ RealType gamma;
+ MeshFunctionPointer pressure;
+
+};
+
+/****
+ * 3D grid
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Function,
+ typename Real,
+ typename Index >
+class DensityBoundaryConditions< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public DensityBoundaryConditionsBase< Function >,
+ public Operator< Meshes::Grid< 3, MeshReal, Device, MeshIndex >,
+ Functions::MeshBoundaryDomain,
+ 3, 3,
+ Real,
+ Index >
+{
+ public:
+
+ typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ typedef Function FunctionType;
+ typedef Functions::MeshFunction< MeshType > MeshFunctionType;
+ typedef Containers::Vector< RealType, DeviceType, IndexType> DofVectorType;
+ typedef Containers::StaticVector< 3, RealType > PointType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef DensityBoundaryConditions< MeshType, Function, Real, Index > ThisType;
+ typedef DensityBoundaryConditionsBase< Function > BaseType;
+ typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
+ typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
+ typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
+
+ template< typename EntityType,
+ typename MeshFunction >
+ __cuda_callable__
+ const RealType operator()( const MeshFunction& u,
+ const EntityType& entity,
+ const RealType& time = 0 ) const
+ {
+ const MeshType& mesh = entity.getMesh();
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ if( entity.getCoordinates().x() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 1, 0, 0 >() ] + entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< -1, 0, 0 >() ] + entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 1, 0 >() ] + entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, -1, 0 >() ] + entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().z() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 1 >() ] + entity.getMesh().getSpaceSteps().z() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ // The following line is commented to avoid compiler warning
+ //if( entity.getCoordinates().z() == entity.getMesh().getDimensions().z() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, -1 >() ] + entity.getMesh().getSpaceSteps().z() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ }
+
+
+ template< typename EntityType >
+ __cuda_callable__
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const EntityType& entity ) const
+ {
+ return 2;
+ }
+
+ template< typename PreimageFunction,
+ typename MeshEntity,
+ typename Matrix,
+ typename Vector >
+ __cuda_callable__
+ void setMatrixElements( const PreimageFunction& u,
+ const MeshEntity& entity,
+ const RealType& time,
+ const RealType& tau,
+ Matrix& matrix,
+ Vector& b ) const
+ {
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
+ if( entity.getCoordinates().x() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 1, 0, 0 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< -1, 0, 0 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 0, 1, 0 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< 0, -1, 0 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().z() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 0, 0, 1 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().z() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().z() == entity.getMesh().getDimensions().z() - 1 )
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< 0, 0, -1 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().z() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ }
+
+ void setCompressibleConservativeVariables(const CompressibleConservativeVariablesPointer& compressibleConservativeVariables)
+ {
+ this->compressibleConservativeVariables = compressibleConservativeVariables;
+ }
+
+ void setTimestep(const RealType timestep )
+ {
+ this->timestep = timestep;
+ }
+
+ void setGamma(const RealType gamma )
+ {
+ this->gamma = gamma;
+ }
+
+ void setPressure(const MeshFunctionPointer& pressure)
+ {
+ this->pressure = pressure;
+ }
+
+ private:
+ CompressibleConservativeVariablesPointer compressibleConservativeVariables;
+ RealType timestep;
+ RealType gamma;
+ MeshFunctionPointer pressure;
+
+};
+
+template< typename Mesh,
+ typename Function,
+ typename Real,
+ typename Index >
+std::ostream& operator << ( std::ostream& str, const DensityBoundaryConditions< Mesh, Function, Real, Index >& bc )
+{
+ str << "Neumann boundary conditions: function = " << bc.getFunction();
+ return str;
+}
+
+} // namespace Operators
+} // namespace TNL
+
diff --git a/examples/flow/1DEnergyBoundaryCondition.h b/examples/flow/1DEnergyBoundaryCondition.h
new file mode 100644
index 0000000000000000000000000000000000000000..bce86634a926e35506d027a66384d59686808580
--- /dev/null
+++ b/examples/flow/1DEnergyBoundaryCondition.h
@@ -0,0 +1,588 @@
+/***************************************************************************
+ IdentityOperator.h - description
+ -------------------
+ begin : Nov 17, 2014
+ copyright : (C) 2014 by Tomas Oberhuber
+ email : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+
+#pragma once
+
+#include <TNL/Functions/FunctionAdapter.h>
+#include "CompressibleConservativeVariables.h"
+
+namespace TNL {
+namespace Operators {
+
+template< typename Mesh,
+ typename Function,
+ typename Real = typename Mesh::RealType,
+ typename Index = typename Mesh::GlobalIndexType >
+class EnergyBoundaryConditions
+{
+
+};
+
+/****
+ * Base
+ */
+template< typename Function>
+class EnergyBoundaryConditionsBase
+{
+ public:
+
+ typedef Function FunctionType;
+
+ static void configSetup( const Config::ConfigDescription& config,
+ const String& prefix = "" )
+ {
+ Function::configSetup( config, prefix );
+ }
+
+ template< typename MeshPointer >
+ bool setup( const MeshPointer& meshPointer,
+ const Config::ParameterContainer& parameters,
+ const String& prefix = "" )
+ {
+ return Functions::FunctionAdapter< typename MeshPointer::ObjectType, FunctionType >::setup( this->function, meshPointer, parameters, prefix );
+ }
+
+ static void configSetup( Config::ConfigDescription& config,
+ const String& prefix = "" )
+ {
+ Function::configSetup( config, prefix );
+ };
+
+ bool setup( const Config::ParameterContainer& parameters,
+ const String& prefix = "" )
+ {
+ return this->function.setup( parameters, prefix );
+ };
+
+ void setFunction( const FunctionType& function )
+ {
+ this->function = function;
+ };
+
+ FunctionType& getFunction()
+ {
+ return this->function;
+ }
+
+ const FunctionType& getFunction() const
+ {
+ return this->function;
+ };
+
+ protected:
+
+ FunctionType function;
+
+
+};
+
+/****
+ * 1D grid
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Function,
+ typename Real,
+ typename Index >
+class EnergyBoundaryConditions< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public EnergyBoundaryConditionsBase< Function >,
+ public Operator< Meshes::Grid< 1, MeshReal, Device, MeshIndex >,
+ Functions::MeshBoundaryDomain,
+ 1, 1,
+ Real,
+ Index >
+{
+ public:
+
+ typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ typedef Function FunctionType;
+ typedef Functions::MeshFunction< MeshType > MeshFunctionType;
+ typedef Containers::Vector< RealType, DeviceType, IndexType> DofVectorType;
+ typedef Containers::StaticVector< 1, RealType > PointType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef EnergyBoundaryConditions< MeshType, Function, Real, Index > ThisType;
+ typedef EnergyBoundaryConditionsBase< Function > BaseType;
+ typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
+ typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
+ typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
+
+ template< typename EntityType,
+ typename MeshFunction >
+ __cuda_callable__
+ const RealType operator()( const MeshFunction& u,
+ const EntityType& entity,
+ const RealType& time = 0 ) const
+ {
+ const MeshType& mesh = entity.getMesh();
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ if( entity.getCoordinates().x() == 0 )
+ return ( (* this->pressure)[ neighborEntities.template getEntityIndex< 0 >() ]
+ / ( this->gamma - 1 )
+ )
+ + 0.5
+ * (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0 >()]
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0 >()]
+ + this->timestep
+ )
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0 >()]
+ + this->timestep
+ );
+ else
+ return u[ neighborEntities.template getEntityIndex< -1 >() ];
+
+ }
+
+
+ template< typename EntityType >
+ __cuda_callable__
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const EntityType& entity ) const
+ {
+ return 2;
+ }
+
+ template< typename PreimageFunction,
+ typename MeshEntity,
+ typename Matrix,
+ typename Vector >
+ __cuda_callable__
+ void setMatrixElements( const PreimageFunction& u,
+ const MeshEntity& entity,
+ const RealType& time,
+ const RealType& tau,
+ Matrix& matrix,
+ Vector& b ) const
+ {
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
+ if( entity.getCoordinates().x() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 1 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ else
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< -1 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ }
+
+ void setTimestep(const RealType timestep )
+ {
+ this->timestep = timestep;
+ }
+
+ void setGamma(const RealType gamma )
+ {
+ this->gamma = gamma;
+ }
+
+ void setCompressibleConservativeVariables(const CompressibleConservativeVariablesPointer& compressibleConservativeVariables)
+ {
+ this->compressibleConservativeVariables = compressibleConservativeVariables;
+ }
+
+ void setPressure(const MeshFunctionPointer& pressure)
+ {
+ this->pressure = pressure;
+ }
+
+ void setCavitySpeed(const RealType cavitySpeed)
+ {
+ this->cavitySpeed = cavitySpeed;
+ }
+
+ private:
+ CompressibleConservativeVariablesPointer compressibleConservativeVariables;
+ RealType timestep;
+ RealType cavitySpeed;
+ RealType gamma;
+ MeshFunctionPointer pressure;
+};
+
+/****
+ * 2D grid
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Function,
+ typename Real,
+ typename Index >
+class EnergyBoundaryConditions< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public EnergyBoundaryConditionsBase< Function >,
+ public Operator< Meshes::Grid< 2, MeshReal, Device, MeshIndex >,
+ Functions::MeshBoundaryDomain,
+ 2, 2,
+ Real,
+ Index >
+
+{
+ public:
+
+ typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ typedef Function FunctionType;
+ typedef Functions::MeshFunction< MeshType > MeshFunctionType;
+ typedef Containers::Vector< RealType, DeviceType, IndexType> DofVectorType;
+ typedef Containers::StaticVector< 2, RealType > PointType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef EnergyBoundaryConditions< MeshType, Function, Real, Index > ThisType;
+ typedef EnergyBoundaryConditionsBase< Function > BaseType;
+ typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
+ typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
+ typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
+
+
+ template< typename EntityType,
+ typename MeshFunction >
+ __cuda_callable__
+ const RealType operator()( const MeshFunction& u,
+ const EntityType& entity,
+ const RealType& time = 0 ) const
+ {
+ const MeshType& mesh = entity.getMesh();
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ if( entity.getCoordinates().x() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ }
+ // The following line is commented to avoid compiler warning
+ //if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ return ( (* this->pressure)[ neighborEntities.template getEntityIndex< 0, 0 >() ]
+ / ( this->gamma - 1 )
+ )
+ + 0.5
+ * (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ * (
+ this->cavitySpeed
+ *
+ this->cavitySpeed
+ +
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ + 0
+ )
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ + 0
+ )
+ );
+ }
+ }
+
+ template< typename EntityType >
+ __cuda_callable__
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const EntityType& entity ) const
+ {
+ return 2;
+ }
+
+ template< typename PreimageFunction,
+ typename MeshEntity,
+ typename Matrix,
+ typename Vector >
+ __cuda_callable__
+ void setMatrixElements( const PreimageFunction& u,
+ const MeshEntity& entity,
+ const RealType& time,
+ const RealType& tau,
+ Matrix& matrix,
+ Vector& b ) const
+ {
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
+ if( entity.getCoordinates().x() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 1, 0 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< -1, 0 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 0, 1 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< 0, -1 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ }
+
+ void setTimestep(const RealType timestep )
+ {
+ this->timestep = timestep;
+ }
+
+ void setGamma(const RealType gamma )
+ {
+ this->gamma = gamma;
+ }
+
+ void setCompressibleConservativeVariables(const CompressibleConservativeVariablesPointer& compressibleConservativeVariables)
+ {
+ this->compressibleConservativeVariables = compressibleConservativeVariables;
+ }
+
+ void setPressure(const MeshFunctionPointer& pressure)
+ {
+ this->pressure = pressure;
+ }
+
+ void setCavitySpeed(const RealType cavitySpeed)
+ {
+ this->cavitySpeed = cavitySpeed;
+ }
+
+ private:
+ CompressibleConservativeVariablesPointer compressibleConservativeVariables;
+ RealType timestep;
+ RealType cavitySpeed;
+ RealType gamma;
+ MeshFunctionPointer pressure;
+};
+
+/****
+ * 3D grid
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Function,
+ typename Real,
+ typename Index >
+class EnergyBoundaryConditions< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public EnergyBoundaryConditionsBase< Function >,
+ public Operator< Meshes::Grid< 3, MeshReal, Device, MeshIndex >,
+ Functions::MeshBoundaryDomain,
+ 3, 3,
+ Real,
+ Index >
+{
+ public:
+
+ typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ typedef Function FunctionType;
+ typedef Functions::MeshFunction< MeshType > MeshFunctionType;
+ typedef Containers::Vector< RealType, DeviceType, IndexType> DofVectorType;
+ typedef Containers::StaticVector< 3, RealType > PointType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef EnergyBoundaryConditions< MeshType, Function, Real, Index > ThisType;
+ typedef EnergyBoundaryConditionsBase< Function > BaseType;
+ typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
+ typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
+ typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
+
+ template< typename EntityType,
+ typename MeshFunction >
+ __cuda_callable__
+ const RealType operator()( const MeshFunction& u,
+ const EntityType& entity,
+ const RealType& time = 0 ) const
+ {
+ const MeshType& mesh = entity.getMesh();
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ if( entity.getCoordinates().x() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 1, 0, 0 >() ] + entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< -1, 0, 0 >() ] + entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 1, 0 >() ] + entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, -1, 0 >() ] + entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().z() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 1 >() ] + entity.getMesh().getSpaceSteps().z() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ // The following line is commented to avoid compiler warning
+ //if( entity.getCoordinates().z() == entity.getMesh().getDimensions().z() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, -1 >() ] + entity.getMesh().getSpaceSteps().z() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ }
+
+
+ template< typename EntityType >
+ __cuda_callable__
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const EntityType& entity ) const
+ {
+ return 2;
+ }
+
+ template< typename PreimageFunction,
+ typename MeshEntity,
+ typename Matrix,
+ typename Vector >
+ __cuda_callable__
+ void setMatrixElements( const PreimageFunction& u,
+ const MeshEntity& entity,
+ const RealType& time,
+ const RealType& tau,
+ Matrix& matrix,
+ Vector& b ) const
+ {
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
+ if( entity.getCoordinates().x() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 1, 0, 0 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< -1, 0, 0 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 0, 1, 0 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< 0, -1, 0 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().z() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 0, 0, 1 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().z() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().z() == entity.getMesh().getDimensions().z() - 1 )
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< 0, 0, -1 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().z() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ }
+
+ void setTimestep(const RealType timestep )
+ {
+ this->timestep = timestep;
+ }
+
+ void setGamma(const RealType gamma )
+ {
+ this->gamma = gamma;
+ }
+
+ void setCompressibleConservativeVariables(const CompressibleConservativeVariablesPointer& compressibleConservativeVariables)
+ {
+ this->compressibleConservativeVariables = compressibleConservativeVariables;
+ }
+
+ void setPressure(const MeshFunctionPointer& pressure)
+ {
+ this->pressure = pressure;
+ }
+
+ void setCavitySpeed(const RealType cavitySpeed)
+ {
+ this->cavitySpeed = cavitySpeed;
+ }
+
+ private:
+ CompressibleConservativeVariablesPointer compressibleConservativeVariables;
+ RealType timestep;
+ RealType cavitySpeed;
+ RealType gamma;
+ MeshFunctionPointer pressure;
+};
+
+template< typename Mesh,
+ typename Function,
+ typename Real,
+ typename Index >
+std::ostream& operator << ( std::ostream& str, const EnergyBoundaryConditions< Mesh, Function, Real, Index >& bc )
+{
+ str << "Neumann boundary conditions: function = " << bc.getFunction();
+ return str;
+}
+
+} // namespace Operators
+} // namespace TNL
+
diff --git a/examples/flow/1DMomentumXBoundaryCondition.h b/examples/flow/1DMomentumXBoundaryCondition.h
new file mode 100644
index 0000000000000000000000000000000000000000..101050bc098b1a71dac9b407e783e156cbf1f9e0
--- /dev/null
+++ b/examples/flow/1DMomentumXBoundaryCondition.h
@@ -0,0 +1,563 @@
+/***************************************************************************
+ IdentityOperator.h - description
+ -------------------
+ begin : Nov 17, 2014
+ copyright : (C) 2014 by Tomas Oberhuber
+ email : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+
+#pragma once
+
+#include <TNL/Functions/FunctionAdapter.h>
+
+namespace TNL {
+namespace Operators {
+
+template< typename Mesh,
+ typename Function,
+ typename Real = typename Mesh::RealType,
+ typename Index = typename Mesh::GlobalIndexType >
+class MomentumXBoundaryConditions
+{
+
+};
+
+/****
+ * Base
+ */
+template< typename Function >
+class MomentumXBoundaryConditionsBase
+{
+ public:
+
+ typedef Function FunctionType;
+
+ static void configSetup( const Config::ConfigDescription& config,
+ const String& prefix = "" )
+ {
+ Function::configSetup( config, prefix );
+ }
+
+ template< typename MeshPointer >
+ bool setup( const MeshPointer& meshPointer,
+ const Config::ParameterContainer& parameters,
+ const String& prefix = "" )
+ {
+ return Functions::FunctionAdapter< typename MeshPointer::ObjectType, FunctionType >::setup( this->function, meshPointer, parameters, prefix );
+ }
+
+ static void configSetup( Config::ConfigDescription& config,
+ const String& prefix = "" )
+ {
+ Function::configSetup( config, prefix );
+ };
+
+ bool setup( const Config::ParameterContainer& parameters,
+ const String& prefix = "" )
+ {
+ return this->function.setup( parameters, prefix );
+ };
+
+ void setFunction( const FunctionType& function )
+ {
+ this->function = function;
+ };
+
+ FunctionType& getFunction()
+ {
+ return this->function;
+ }
+
+ const FunctionType& getFunction() const
+ {
+ return this->function;
+ };
+
+ protected:
+
+ FunctionType function;
+
+};
+
+/****
+ * 1D grid
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Function,
+ typename Real,
+ typename Index >
+class MomentumXBoundaryConditions< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumXBoundaryConditionsBase< Function >,
+ public Operator< Meshes::Grid< 1, MeshReal, Device, MeshIndex >,
+ Functions::MeshBoundaryDomain,
+ 1, 1,
+ Real,
+ Index >
+{
+ public:
+
+ typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ typedef Function FunctionType;
+ typedef Functions::MeshFunction< MeshType > MeshFunctionType;
+ typedef Containers::Vector< RealType, DeviceType, IndexType> DofVectorType;
+ typedef Containers::StaticVector< 1, RealType > PointType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef MomentumXBoundaryConditions< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumXBoundaryConditionsBase< Function > BaseType;
+ typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
+ typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
+ typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
+
+ template< typename EntityType,
+ typename MeshFunction >
+ __cuda_callable__
+ const RealType operator()( const MeshFunction& u,
+ const EntityType& entity,
+ const RealType& time = 0 ) const
+ {
+ const MeshType& mesh = entity.getMesh();
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ if( entity.getCoordinates().x() == 0 )
+ return (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0 >()]
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0 >()]
+ + this->timestep
+ );
+ else
+ return u[ neighborEntities.template getEntityIndex< -1 >() ];
+
+ }
+
+
+ template< typename EntityType >
+ __cuda_callable__
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const EntityType& entity ) const
+ {
+ return 2;
+ }
+
+ template< typename PreimageFunction,
+ typename MeshEntity,
+ typename Matrix,
+ typename Vector >
+ __cuda_callable__
+ void setMatrixElements( const PreimageFunction& u,
+ const MeshEntity& entity,
+ const RealType& time,
+ const RealType& tau,
+ Matrix& matrix,
+ Vector& b ) const
+ {
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
+ if( entity.getCoordinates().x() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 1 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ else
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< -1 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ }
+
+ void setTimestep(const RealType timestep )
+ {
+ this->timestep = timestep;
+ }
+
+ void setGamma(const RealType gamma )
+ {
+ this->gamma = gamma;
+ }
+
+ void setCompressibleConservativeVariables(const CompressibleConservativeVariablesPointer& compressibleConservativeVariables)
+ {
+ this->compressibleConservativeVariables = compressibleConservativeVariables;
+ }
+
+ void setPressure(const MeshFunctionPointer& pressure)
+ {
+ this->pressure = pressure;
+ }
+
+ void setCavitySpeed(const RealType cavitySpeed)
+ {
+ this->cavitySpeed = cavitySpeed;
+ }
+
+ private:
+ CompressibleConservativeVariablesPointer compressibleConservativeVariables;
+ RealType timestep;
+ RealType cavitySpeed;
+ RealType gamma;
+ MeshFunctionPointer pressure;
+};
+
+/****
+ * 2D grid
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Function,
+ typename Real,
+ typename Index >
+class MomentumXBoundaryConditions< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumXBoundaryConditionsBase< Function >,
+ public Operator< Meshes::Grid< 2, MeshReal, Device, MeshIndex >,
+ Functions::MeshBoundaryDomain,
+ 2, 2,
+ Real,
+ Index >
+
+{
+ public:
+
+ typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ typedef Function FunctionType;
+ typedef Functions::MeshFunction< MeshType > MeshFunctionType;
+ typedef Containers::Vector< RealType, DeviceType, IndexType> DofVectorType;
+ typedef Containers::StaticVector< 2, RealType > PointType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef MomentumXBoundaryConditions< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumXBoundaryConditionsBase< Function > BaseType;
+ typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
+ typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
+ typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
+
+
+ template< typename EntityType,
+ typename MeshFunction >
+ __cuda_callable__
+ const RealType operator()( const MeshFunction& u,
+ const EntityType& entity,
+ const RealType& time = 0 ) const
+ {
+ const MeshType& mesh = entity.getMesh();
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ if( entity.getCoordinates().x() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< -0, 0 >() ];
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ }
+ // The following line is commented to avoid compiler warning
+ //if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ return (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ * (
+ this->cavitySpeed
+ );
+ }
+ }
+
+ template< typename EntityType >
+ __cuda_callable__
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const EntityType& entity ) const
+ {
+ return 2;
+ }
+
+ template< typename PreimageFunction,
+ typename MeshEntity,
+ typename Matrix,
+ typename Vector >
+ __cuda_callable__
+ void setMatrixElements( const PreimageFunction& u,
+ const MeshEntity& entity,
+ const RealType& time,
+ const RealType& tau,
+ Matrix& matrix,
+ Vector& b ) const
+ {
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
+ if( entity.getCoordinates().x() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 1, 0 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< -1, 0 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 0, 1 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< 0, -1 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ }
+
+ void setTimestep(const RealType timestep )
+ {
+ this->timestep = timestep;
+ }
+
+ void setGamma(const RealType gamma )
+ {
+ this->gamma = gamma;
+ }
+
+ void setCompressibleConservativeVariables(const CompressibleConservativeVariablesPointer& compressibleConservativeVariables)
+ {
+ this->compressibleConservativeVariables = compressibleConservativeVariables;
+ }
+
+ void setPressure(const MeshFunctionPointer& pressure)
+ {
+ this->pressure = pressure;
+ }
+
+ void setCavitySpeed(const RealType cavitySpeed)
+ {
+ this->cavitySpeed = cavitySpeed;
+ }
+
+ private:
+ CompressibleConservativeVariablesPointer compressibleConservativeVariables;
+ RealType timestep;
+ RealType cavitySpeed;
+ RealType gamma;
+ MeshFunctionPointer pressure;
+};
+
+/****
+ * 3D grid
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Function,
+ typename Real,
+ typename Index >
+class MomentumXBoundaryConditions< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumXBoundaryConditionsBase< Function >,
+ public Operator< Meshes::Grid< 3, MeshReal, Device, MeshIndex >,
+ Functions::MeshBoundaryDomain,
+ 3, 3,
+ Real,
+ Index >
+{
+ public:
+
+ typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ typedef Function FunctionType;
+ typedef Functions::MeshFunction< MeshType > MeshFunctionType;
+ typedef Containers::Vector< RealType, DeviceType, IndexType> DofVectorType;
+ typedef Containers::StaticVector< 3, RealType > PointType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef MomentumXBoundaryConditions< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumXBoundaryConditionsBase< Function > BaseType;
+ typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
+ typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
+ typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
+
+ template< typename EntityType,
+ typename MeshFunction >
+ __cuda_callable__
+ const RealType operator()( const MeshFunction& u,
+ const EntityType& entity,
+ const RealType& time = 0 ) const
+ {
+ const MeshType& mesh = entity.getMesh();
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ if( entity.getCoordinates().x() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 1, 0, 0 >() ] + entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< -1, 0, 0 >() ] + entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 1, 0 >() ] + entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, -1, 0 >() ] + entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().z() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 1 >() ] + entity.getMesh().getSpaceSteps().z() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ // The following line is commented to avoid compiler warning
+ //if( entity.getCoordinates().z() == entity.getMesh().getDimensions().z() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, -1 >() ] + entity.getMesh().getSpaceSteps().z() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ }
+
+
+ template< typename EntityType >
+ __cuda_callable__
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const EntityType& entity ) const
+ {
+ return 2;
+ }
+
+ template< typename PreimageFunction,
+ typename MeshEntity,
+ typename Matrix,
+ typename Vector >
+ __cuda_callable__
+ void setMatrixElements( const PreimageFunction& u,
+ const MeshEntity& entity,
+ const RealType& time,
+ const RealType& tau,
+ Matrix& matrix,
+ Vector& b ) const
+ {
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
+ if( entity.getCoordinates().x() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 1, 0, 0 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< -1, 0, 0 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 0, 1, 0 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< 0, -1, 0 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().z() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 0, 0, 1 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().z() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().z() == entity.getMesh().getDimensions().z() - 1 )
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< 0, 0, -1 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().z() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ }
+
+ void setTimestep(const RealType timestep )
+ {
+ this->timestep = timestep;
+ }
+
+ void setGamma(const RealType gamma )
+ {
+ this->gamma = gamma;
+ }
+
+ void setCompressibleConservativeVariables(const CompressibleConservativeVariablesPointer& compressibleConservativeVariables)
+ {
+ this->compressibleConservativeVariables = compressibleConservativeVariables;
+ }
+
+ void setPressure(const MeshFunctionPointer& pressure)
+ {
+ this->pressure = pressure;
+ }
+
+ void setCavitySpeed(const RealType cavitySpeed)
+ {
+ this->cavitySpeed = cavitySpeed;
+ }
+
+ private:
+ CompressibleConservativeVariablesPointer compressibleConservativeVariables;
+ RealType timestep;
+ RealType cavitySpeed;
+ RealType gamma;
+ MeshFunctionPointer pressure;
+};
+
+template< typename Mesh,
+ typename Function,
+ typename Real,
+ typename Index >
+std::ostream& operator << ( std::ostream& str, const MomentumXBoundaryConditions< Mesh, Function, Real, Index >& bc )
+{
+ str << "Neumann boundary conditions: function = " << bc.getFunction();
+ return str;
+}
+
+} // namespace Operators
+} // namespace TNL
+
diff --git a/examples/flow/1DMomentumYBoundaryCondition.h b/examples/flow/1DMomentumYBoundaryCondition.h
new file mode 100644
index 0000000000000000000000000000000000000000..1a8637bbac3a5781e83c7e7630c57c3a3835c713
--- /dev/null
+++ b/examples/flow/1DMomentumYBoundaryCondition.h
@@ -0,0 +1,560 @@
+/***************************************************************************
+ IdentityOperator.h - description
+ -------------------
+ begin : Nov 17, 2014
+ copyright : (C) 2014 by Tomas Oberhuber
+ email : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+
+#pragma once
+
+#include <TNL/Functions/FunctionAdapter.h>
+
+namespace TNL {
+namespace Operators {
+
+template< typename Mesh,
+ typename Function,
+ typename Real = typename Mesh::RealType,
+ typename Index = typename Mesh::GlobalIndexType >
+class MomentumYBoundaryConditions
+{
+
+};
+
+/****
+ * Base
+ */
+template< typename Function >
+class MomentumYBoundaryConditionsBase
+{
+ public:
+
+ typedef Function FunctionType;
+
+ static void configSetup( const Config::ConfigDescription& config,
+ const String& prefix = "" )
+ {
+ Function::configSetup( config, prefix );
+ }
+
+ template< typename MeshPointer >
+ bool setup( const MeshPointer& meshPointer,
+ const Config::ParameterContainer& parameters,
+ const String& prefix = "" )
+ {
+ return Functions::FunctionAdapter< typename MeshPointer::ObjectType, FunctionType >::setup( this->function, meshPointer, parameters, prefix );
+ }
+
+ static void configSetup( Config::ConfigDescription& config,
+ const String& prefix = "" )
+ {
+ Function::configSetup( config, prefix );
+ };
+
+ bool setup( const Config::ParameterContainer& parameters,
+ const String& prefix = "" )
+ {
+ return this->function.setup( parameters, prefix );
+ };
+
+ void setFunction( const FunctionType& function )
+ {
+ this->function = function;
+ };
+
+ FunctionType& getFunction()
+ {
+ return this->function;
+ }
+
+ const FunctionType& getFunction() const
+ {
+ return this->function;
+ };
+
+ protected:
+
+ FunctionType function;
+
+};
+
+/****
+ * 1D grid
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Function,
+ typename Real,
+ typename Index >
+class MomentumYBoundaryConditions< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumYBoundaryConditionsBase< Function >,
+ public Operator< Meshes::Grid< 1, MeshReal, Device, MeshIndex >,
+ Functions::MeshBoundaryDomain,
+ 1, 1,
+ Real,
+ Index >
+{
+ public:
+
+ typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ typedef Function FunctionType;
+ typedef Functions::MeshFunction< MeshType > MeshFunctionType;
+ typedef Containers::Vector< RealType, DeviceType, IndexType> DofVectorType;
+ typedef Containers::StaticVector< 1, RealType > PointType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef MomentumYBoundaryConditions< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumYBoundaryConditionsBase< Function > BaseType;
+ typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
+ typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
+ typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
+
+ template< typename EntityType,
+ typename MeshFunction >
+ __cuda_callable__
+ const RealType operator()( const MeshFunction& u,
+ const EntityType& entity,
+ const RealType& time = 0 ) const
+ {
+ const MeshType& mesh = entity.getMesh();
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ if( entity.getCoordinates().x() == 0 )
+ return (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0 >()]
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0 >()]
+ + this->timestep
+ );
+ else
+ return u[ neighborEntities.template getEntityIndex< -1 >() ];
+
+ }
+
+
+ template< typename EntityType >
+ __cuda_callable__
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const EntityType& entity ) const
+ {
+ return 2;
+ }
+
+ template< typename PreimageFunction,
+ typename MeshEntity,
+ typename Matrix,
+ typename Vector >
+ __cuda_callable__
+ void setMatrixElements( const PreimageFunction& u,
+ const MeshEntity& entity,
+ const RealType& time,
+ const RealType& tau,
+ Matrix& matrix,
+ Vector& b ) const
+ {
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
+ if( entity.getCoordinates().x() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 1 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ else
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< -1 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ }
+
+ void setTimestep(const RealType timestep )
+ {
+ this->timestep = timestep;
+ }
+
+ void setGamma(const RealType gamma )
+ {
+ this->gamma = gamma;
+ }
+
+ void setCompressibleConservativeVariables(const CompressibleConservativeVariablesPointer& compressibleConservativeVariables)
+ {
+ this->compressibleConservativeVariables = compressibleConservativeVariables;
+ }
+
+ void setPressure(const MeshFunctionPointer& pressure)
+ {
+ this->pressure = pressure;
+ }
+
+ void setCavitySpeed(const RealType cavitySpeed)
+ {
+ this->cavitySpeed = cavitySpeed;
+ }
+
+ private:
+ CompressibleConservativeVariablesPointer compressibleConservativeVariables;
+ RealType timestep;
+ RealType cavitySpeed;
+ RealType gamma;
+ MeshFunctionPointer pressure;
+};
+
+/****
+ * 2D grid
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Function,
+ typename Real,
+ typename Index >
+class MomentumYBoundaryConditions< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumYBoundaryConditionsBase< Function >,
+ public Operator< Meshes::Grid< 2, MeshReal, Device, MeshIndex >,
+ Functions::MeshBoundaryDomain,
+ 2, 2,
+ Real,
+ Index >
+
+{
+ public:
+
+ typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ typedef Function FunctionType;
+ typedef Functions::MeshFunction< MeshType > MeshFunctionType;
+ typedef Containers::Vector< RealType, DeviceType, IndexType> DofVectorType;
+ typedef Containers::StaticVector< 2, RealType > PointType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef MomentumYBoundaryConditions< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumYBoundaryConditionsBase< Function > BaseType;
+ typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
+ typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
+ typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
+
+
+ template< typename EntityType,
+ typename MeshFunction >
+ __cuda_callable__
+ const RealType operator()( const MeshFunction& u,
+ const EntityType& entity,
+ const RealType& time = 0 ) const
+ {
+ const MeshType& mesh = entity.getMesh();
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ if( entity.getCoordinates().x() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< -0, 0 >() ];
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ }
+ // The following line is commented to avoid compiler warning
+ //if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ }
+ }
+
+ template< typename EntityType >
+ __cuda_callable__
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const EntityType& entity ) const
+ {
+ return 2;
+ }
+
+ template< typename PreimageFunction,
+ typename MeshEntity,
+ typename Matrix,
+ typename Vector >
+ __cuda_callable__
+ void setMatrixElements( const PreimageFunction& u,
+ const MeshEntity& entity,
+ const RealType& time,
+ const RealType& tau,
+ Matrix& matrix,
+ Vector& b ) const
+ {
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
+ if( entity.getCoordinates().x() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 1, 0 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< -1, 0 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 0, 1 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< 0, -1 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ }
+
+ void setTimestep(const RealType timestep )
+ {
+ this->timestep = timestep;
+ }
+
+ void setGamma(const RealType gamma )
+ {
+ this->gamma = gamma;
+ }
+
+ void setCompressibleConservativeVariables(const CompressibleConservativeVariablesPointer& compressibleConservativeVariables)
+ {
+ this->compressibleConservativeVariables = compressibleConservativeVariables;
+ }
+
+ void setPressure(const MeshFunctionPointer& pressure)
+ {
+ this->pressure = pressure;
+ }
+
+ void setCavitySpeed(const RealType cavitySpeed)
+ {
+ this->cavitySpeed = cavitySpeed;
+ }
+
+ private:
+ CompressibleConservativeVariablesPointer compressibleConservativeVariables;
+ RealType timestep;
+ RealType cavitySpeed;
+ RealType gamma;
+ MeshFunctionPointer pressure;
+};
+
+/****
+ * 3D grid
+ */
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Function,
+ typename Real,
+ typename Index >
+class MomentumYBoundaryConditions< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumYBoundaryConditionsBase< Function >,
+ public Operator< Meshes::Grid< 3, MeshReal, Device, MeshIndex >,
+ Functions::MeshBoundaryDomain,
+ 3, 3,
+ Real,
+ Index >
+{
+ public:
+
+ typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ typedef Function FunctionType;
+ typedef Functions::MeshFunction< MeshType > MeshFunctionType;
+ typedef Containers::Vector< RealType, DeviceType, IndexType> DofVectorType;
+ typedef Containers::StaticVector< 3, RealType > PointType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef MomentumYBoundaryConditions< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumYBoundaryConditionsBase< Function > BaseType;
+ typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
+ typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
+ typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
+
+ template< typename EntityType,
+ typename MeshFunction >
+ __cuda_callable__
+ const RealType operator()( const MeshFunction& u,
+ const EntityType& entity,
+ const RealType& time = 0 ) const
+ {
+ const MeshType& mesh = entity.getMesh();
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ if( entity.getCoordinates().x() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 1, 0, 0 >() ] + entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< -1, 0, 0 >() ] + entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 1, 0 >() ] + entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, -1, 0 >() ] + entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().z() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 1 >() ] + entity.getMesh().getSpaceSteps().z() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ // The following line is commented to avoid compiler warning
+ //if( entity.getCoordinates().z() == entity.getMesh().getDimensions().z() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, -1 >() ] + entity.getMesh().getSpaceSteps().z() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ }
+
+
+ template< typename EntityType >
+ __cuda_callable__
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const EntityType& entity ) const
+ {
+ return 2;
+ }
+
+ template< typename PreimageFunction,
+ typename MeshEntity,
+ typename Matrix,
+ typename Vector >
+ __cuda_callable__
+ void setMatrixElements( const PreimageFunction& u,
+ const MeshEntity& entity,
+ const RealType& time,
+ const RealType& tau,
+ Matrix& matrix,
+ Vector& b ) const
+ {
+ const auto& neighborEntities = entity.getNeighborEntities();
+ const IndexType& index = entity.getIndex();
+ typename Matrix::MatrixRow matrixRow = matrix.getRow( index );
+ if( entity.getCoordinates().x() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 1, 0, 0 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< -1, 0, 0 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().x() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 0, 1, 0 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< 0, -1, 0 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().y() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().z() == 0 )
+ {
+ matrixRow.setElement( 0, index, 1.0 );
+ matrixRow.setElement( 1, neighborEntities.template getEntityIndex< 0, 0, 1 >(), -1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().z() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ if( entity.getCoordinates().z() == entity.getMesh().getDimensions().z() - 1 )
+ {
+ matrixRow.setElement( 0, neighborEntities.template getEntityIndex< 0, 0, -1 >(), -1.0 );
+ matrixRow.setElement( 1, index, 1.0 );
+ b[ index ] = entity.getMesh().getSpaceSteps().z() *
+ Functions::FunctionAdapter< MeshType, FunctionType >::getValue( this->function, entity, time );
+ }
+ }
+
+ void setTimestep(const RealType timestep )
+ {
+ this->timestep = timestep;
+ }
+
+ void setGamma(const RealType gamma )
+ {
+ this->gamma = gamma;
+ }
+
+ void setCompressibleConservativeVariables(const CompressibleConservativeVariablesPointer& compressibleConservativeVariables)
+ {
+ this->compressibleConservativeVariables = compressibleConservativeVariables;
+ }
+
+ void setPressure(const MeshFunctionPointer& pressure)
+ {
+ this->pressure = pressure;
+ }
+
+ void setCavitySpeed(const RealType cavitySpeed)
+ {
+ this->cavitySpeed = cavitySpeed;
+ }
+
+ private:
+ CompressibleConservativeVariablesPointer compressibleConservativeVariables;
+ RealType timestep;
+ RealType cavitySpeed;
+ RealType gamma;
+ MeshFunctionPointer pressure;
+};
+
+template< typename Mesh,
+ typename Function,
+ typename Real,
+ typename Index >
+std::ostream& operator << ( std::ostream& str, const MomentumYBoundaryConditions< Mesh, Function, Real, Index >& bc )
+{
+ str << "Neumann boundary conditions: function = " << bc.getFunction();
+ return str;
+}
+
+} // namespace Operators
+} // namespace TNL
+
diff --git a/examples/flow/LaxFridrichs.h b/examples/flow/LaxFridrichs.h
index 7c2230f33993b3991ff491230e63935286f1b964..34fbda09a259a7cbab5a6813737b63466ac23df0 100644
--- a/examples/flow/LaxFridrichs.h
+++ b/examples/flow/LaxFridrichs.h
@@ -68,10 +68,10 @@ class LaxFridrichs
const String& prefix = "" )
{
this->dynamicalViscosity = parameters.getParameter< double >( prefix + "dynamical-viscosity" );
- this->momentumXOperatorPointer->setDynamicalViscosity( artificialViscosity );
- this->momentumYOperatorPointer->setDynamicalViscosity( artificialViscosity );
- this->momentumZOperatorPointer->setDynamicalViscosity( artificialViscosity );
- this->energyOperatorPointer->setDynamicalViscosity( artificialViscosity );
+ this->momentumXOperatorPointer->setDynamicalViscosity( dynamicalViscosity );
+ this->momentumYOperatorPointer->setDynamicalViscosity( dynamicalViscosity );
+ this->momentumZOperatorPointer->setDynamicalViscosity( dynamicalViscosity );
+ this->energyOperatorPointer->setDynamicalViscosity( dynamicalViscosity );
this->artificialViscosity = parameters.getParameter< double >( prefix + "numerical-viscosity" );
this->continuityOperatorPointer->setArtificialViscosity( artificialViscosity );
this->momentumXOperatorPointer->setArtificialViscosity( artificialViscosity );
diff --git a/examples/flow/LaxFridrichsEnergy.h b/examples/flow/LaxFridrichsEnergy.h
index 546c41b1dc28a76adbec659077e84f20ac4a5175..7714edf0135d1f148e70765f3ce96ecd4fa64fe0 100644
--- a/examples/flow/LaxFridrichsEnergy.h
+++ b/examples/flow/LaxFridrichsEnergy.h
@@ -136,7 +136,7 @@ class LaxFridrichsEnergy< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real,
// 1D uT_11_x
- 4.0 / 3.0 * ( velocity_x_east * velocity_x_center - velocity_x_center * velocity_x_west
- velocity_x_center * velocity_x_center + velocity_x_west * velocity_x_west
- ) * hxSquareInverse
+ ) * hxSquareInverse * 4
* this->dynamicalViscosity;
}
@@ -235,34 +235,34 @@ class LaxFridrichsEnergy< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real,
// 2D uT_11_x
- ( 4.0 / 3.0 * ( velocity_x_east * velocity_x_center - velocity_x_center * velocity_x_west
- velocity_x_center * velocity_x_center + velocity_x_west * velocity_x_west
- ) * hxSquareInverse
+ ) * hxSquareInverse * 4
- 2.0 / 3.0 * ( velocity_y_northEast * velocity_x_east - velocity_y_southEast * velocity_x_east
- velocity_y_northWest * velocity_x_west + velocity_y_southWest * velocity_x_west
- ) * hxInverse * hyInverse
+ ) * hxInverse * hyInverse * 4
) * this->dynamicalViscosity
// vT_21_x
- ( ( velocity_y_northEast * velocity_y_east - velocity_y_southEast * velocity_y_east
- velocity_y_northWest * velocity_y_west + velocity_y_southWest * velocity_y_west
- ) * hxInverse * hyInverse
+ ) * hxInverse * hyInverse * 4
+ ( velocity_x_east * velocity_y_center - velocity_x_center * velocity_y_west
- velocity_x_center * velocity_y_center + velocity_x_west * velocity_y_west
- ) * hxSquareInverse
+ ) * hxSquareInverse * 4
) * this->dynamicalViscosity
// uT_12_y
- ( ( velocity_x_northEast * velocity_x_north - velocity_x_southEast * velocity_x_south
- velocity_x_northWest * velocity_x_north + velocity_x_southWest * velocity_x_south
- ) * hxInverse * hyInverse
+ ) * hxInverse * hyInverse * 4
+ ( velocity_y_north * velocity_x_center - velocity_y_center * velocity_x_south
- velocity_y_center * velocity_x_center + velocity_y_south * velocity_x_south
- ) * hySquareInverse
+ ) * hySquareInverse * 4
) * this->dynamicalViscosity
// 2D vT_22_y
- ( 4.0 / 3.0 * ( velocity_y_north * velocity_y_center - velocity_y_center * velocity_y_south
- velocity_y_center * velocity_y_center + velocity_y_south * velocity_y_south
- ) * hySquareInverse
+ ) * hySquareInverse * 4
- 2.0 / 3.0 * ( velocity_x_northEast * velocity_y_north - velocity_x_southEast * velocity_y_east
- velocity_x_northWest * velocity_y_north + velocity_x_southWest * velocity_y_west
- ) * hxInverse * hyInverse
+ ) * hxInverse * hyInverse * 4
) * this->dynamicalViscosity;
}
@@ -408,83 +408,83 @@ class LaxFridrichsEnergy< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real,
// 3D uT_11_x
- ( 4.0 / 3.0 * ( velocity_x_east * velocity_x_center - velocity_x_center * velocity_x_west
- velocity_x_center * velocity_x_center - velocity_x_west * velocity_x_west
- ) * hxSquareInverse
+ ) * hxSquareInverse * 4
- 2.0 / 3.0 * ( velocity_y_northEast * velocity_x_east - velocity_y_southEast * velocity_x_east
- velocity_y_northWest * velocity_x_west + velocity_y_southWest * velocity_x_west
- ) * hxInverse * hyInverse
+ ) * hxInverse * hyInverse * 4
- 2.0 / 3.0 * ( velocity_z_upEast * velocity_x_east - velocity_z_downEast * velocity_x_east
- velocity_z_upWest * velocity_x_west + velocity_z_downWest * velocity_x_west
- ) * hxInverse * hzInverse
+ ) * hxInverse * hzInverse * 4
) * this->dynamicalViscosity
// vT_21_x
- ( ( velocity_y_northEast * velocity_y_east - velocity_y_southEast * velocity_y_east
- velocity_y_northWest * velocity_y_west + velocity_y_southWest * velocity_y_west
- ) * hxInverse * hyInverse
+ ) * hxInverse * hyInverse * 4
+ ( velocity_x_east * velocity_y_center - velocity_x_center * velocity_y_west
- velocity_x_center * velocity_y_center - velocity_x_west * velocity_y_west
- ) * hxSquareInverse
+ ) * hxSquareInverse * 4
) * this->dynamicalViscosity
// wT_31_x
- ( ( velocity_z_upEast * velocity_z_east - velocity_z_downEast * velocity_z_east
- velocity_z_upWest * velocity_z_west + velocity_z_downWest * velocity_z_west
- ) * hxInverse * hzInverse
+ ) * hxInverse * hzInverse * 4
+ ( velocity_x_east * velocity_z_center - velocity_x_center * velocity_z_east
- velocity_x_center * velocity_z_center - velocity_x_west * velocity_z_west
- ) * hxSquareInverse
+ ) * hxSquareInverse * 4
) * this->dynamicalViscosity
// uT_12_y
- ( ( velocity_x_northEast * velocity_x_north - velocity_x_southEast * velocity_x_south
- velocity_x_northWest * velocity_x_north + velocity_x_southWest * velocity_x_south
- ) * hxInverse * hyInverse
+ ) * hxInverse * hyInverse * 4
+ ( velocity_y_north * velocity_x_center - velocity_y_center * velocity_x_south
+ velocity_y_center * velocity_x_center + velocity_y_south * velocity_x_south
- ) * hySquareInverse
+ ) * hySquareInverse * 4
) * this->dynamicalViscosity
// 3D vT_22_y
- ( 4.0 / 3.0 * ( velocity_y_north * velocity_y_center - velocity_y_center * velocity_y_south
- velocity_y_center * velocity_y_center + velocity_y_south * velocity_y_south
- ) * hySquareInverse
+ ) * hySquareInverse * 4
- 2.0 / 3.0 * ( velocity_x_northEast * velocity_y_north - velocity_x_southEast * velocity_y_south
- velocity_x_northWest * velocity_y_north + velocity_x_southWest * velocity_y_south
- ) * hxInverse * hyInverse
+ ) * hxInverse * hyInverse * 4
- 2.0 / 3.0 * ( velocity_z_upNorth * velocity_y_north - velocity_z_downNorth * velocity_y_north
- velocity_z_upSouth * velocity_y_south + velocity_z_downSouth * velocity_y_south
- ) * hyInverse * hzInverse
+ ) * hyInverse * hzInverse * 4
) * this->dynamicalViscosity
// wT_32_y
- ( ( velocity_z_upNorth * velocity_z_north - velocity_z_downNorth * velocity_y_north
- velocity_z_upSouth * velocity_z_south + velocity_z_downSouth * velocity_z_south
- ) * hyInverse * hzInverse
+ ) * hyInverse * hzInverse * 4
+ ( velocity_y_north * velocity_z_center - velocity_y_center * velocity_z_south
- velocity_y_center * velocity_z_center + velocity_y_south * velocity_z_south
- ) * hySquareInverse
+ ) * hySquareInverse * 4
) * this->dynamicalViscosity
// uT_13_z
- ( ( velocity_z_up * velocity_x_center - velocity_z_center * velocity_x_center
- velocity_z_center * velocity_x_down + velocity_z_down * velocity_x_down
- ) * hzSquareInverse
+ ) * hzSquareInverse * 4
+ ( velocity_x_upEast * velocity_x_up - velocity_x_downEast * velocity_x_down
- velocity_x_upWest * velocity_x_up + velocity_x_downWest * velocity_x_down
- ) * hxInverse * hzInverse
+ ) * hxInverse * hzInverse * 4
) * this->dynamicalViscosity
// T_23_z
- ( ( velocity_y_upNorth * velocity_y_up - velocity_y_downNorth * velocity_y_down
- velocity_y_upSouth * velocity_y_up + velocity_y_downSouth * velocity_y_down
- ) * hyInverse * hzInverse
+ ) * hyInverse * hzInverse * 4
+ ( velocity_z_up * velocity_y_center - velocity_z_center * velocity_y_down
- velocity_z_center * velocity_y_center + velocity_z_down * velocity_y_down
- ) * hzSquareInverse
+ ) * hzSquareInverse * 4
) * this->dynamicalViscosity
// 3D T_33_z
- ( 4.0 / 3.0 * ( velocity_z_up * velocity_z_center - velocity_z_center * velocity_z_center
- velocity_z_center * velocity_z_center + velocity_z_down * velocity_z_down
- ) * hzSquareInverse
+ ) * hzSquareInverse * 4
- 2.0 / 3.0 * ( velocity_y_upNorth * velocity_z_up - velocity_y_downNorth * velocity_z_down
- velocity_y_upSouth * velocity_z_up + velocity_y_downSouth * velocity_z_down
- ) * hyInverse * hzInverse
+ ) * hyInverse * hzInverse * 4
- 2.0 / 3.0 * ( velocity_x_upEast * velocity_z_up - velocity_x_downEast * velocity_z_down
- velocity_x_upWest * velocity_z_up + velocity_x_downWest * velocity_z_down
- ) * hxInverse * hzInverse
+ ) * hxInverse * hzInverse * 4
) * this->dynamicalViscosity;
}
diff --git a/examples/flow/LaxFridrichsMomentumX.h b/examples/flow/LaxFridrichsMomentumX.h
index d24b427f946798afb1a1f765194c71282111b218..d1c74e605f8637090901805c9beb3277232ec19b 100644
--- a/examples/flow/LaxFridrichsMomentumX.h
+++ b/examples/flow/LaxFridrichsMomentumX.h
@@ -82,7 +82,7 @@ class LaxFridrichsMomentumX< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Rea
-( rho_u[ west ] * velocity_x_west + pressure_west ) ) * hxInverse
// 1D T_11_x
- 4.0 / 3.0 *( velocity_x_east - 2 * velocity_x_center + velocity_x_west
- ) * hxSquareInverse
+ ) * hxSquareInverse * 4
* this->dynamicalViscosity;
}
@@ -183,15 +183,15 @@ class LaxFridrichsMomentumX< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Rea
- ( rho_u[ south ] * velocity_y_south ) ) * hyInverse )
// 2D T_11_x
- ( 4.0 / 3.0 * ( velocity_x_east - 2 * velocity_x_center + velocity_x_west
- ) * hxSquareInverse
+ ) * hxSquareInverse * 4
- 2.0 / 3.0 * ( velocity_y_northEast - velocity_y_southEast - velocity_y_northWest + velocity_y_southWest
- ) * hxInverse * hyInverse
+ ) * hxInverse * hyInverse * 4
) * this->dynamicalViscosity
// T_21_x
- ( ( velocity_y_northEast - velocity_y_southEast - velocity_y_northWest + velocity_y_southWest
- ) * hxInverse * hyInverse
+ ) * hxInverse * hyInverse * 4
+ ( velocity_x_east - 2 * velocity_x_center + velocity_x_west
- ) * hxInverse * hyInverse
+ ) * hxInverse * hyInverse * 4
) * this->dynamicalViscosity;
}
@@ -331,23 +331,23 @@ class LaxFridrichsMomentumX< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real
- ( rho_u[ down ] * velocity_z_down ) )* hzInverse )
// 3D T_11_x
- ( 4.0 / 3.0 * ( velocity_x_east - 2 * velocity_x_center + velocity_x_west
- ) * hxSquareInverse
+ ) * hxSquareInverse * 4
- 2.0 / 3.0 * ( velocity_y_northEast - velocity_y_southEast - velocity_y_northWest + velocity_y_southWest
- ) * hxInverse * hyInverse
+ ) * hxInverse * hyInverse * 4
- 2.0 / 3.0 * ( velocity_z_upEast - velocity_z_downEast - velocity_z_upWest + velocity_z_downWest
- ) * hxInverse * hzInverse
+ ) * hxInverse * hzInverse * 4
) * this->dynamicalViscosity
// T_21_x
- ( ( velocity_y_northEast - velocity_y_southEast - velocity_y_northWest + velocity_y_southWest
- ) * hxInverse * hyInverse
+ ) * hxInverse * hyInverse * 4
+ ( velocity_x_east - 2 * velocity_x_center + velocity_x_west
- ) * hxInverse * hyInverse
+ ) * hxInverse * hyInverse * 4
) * this->dynamicalViscosity
// T_31_x
- ( ( velocity_z_upEast - velocity_z_downEast - velocity_z_upWest + velocity_z_downWest
- ) * hxInverse * hzInverse
+ ) * hxInverse * hzInverse * 4
+ ( velocity_x_east - 2 * velocity_x_center + velocity_x_west
- ) * hxInverse * hyInverse
+ ) * hxInverse * hyInverse * 4
) * this->dynamicalViscosity;
}
diff --git a/examples/flow/LaxFridrichsMomentumY.h b/examples/flow/LaxFridrichsMomentumY.h
index 4c5c5f4b883c15147c380a33bbb7c9dd7f68618b..c0f2ead15f974cec2d9b0ddd1795e70bfd6031c4 100644
--- a/examples/flow/LaxFridrichsMomentumY.h
+++ b/examples/flow/LaxFridrichsMomentumY.h
@@ -167,15 +167,15 @@ class LaxFridrichsMomentumY< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Rea
- ( rho_v[ south ] * velocity_y_south + pressure_south ) )* hyInverse )
// 2D T_22_y
- ( 4.0 / 3.0 * ( velocity_y_north - 2 * velocity_y_center + velocity_y_south
- ) * hySquareInverse
+ ) * hySquareInverse * 4
- 2.0 / 3.0 * ( velocity_x_northEast - velocity_x_southEast - velocity_x_northWest + velocity_x_southWest
- ) * hxInverse * hyInverse
+ ) * hxInverse * hyInverse * 4
) * this->dynamicalViscosity
// T_12_y
- ( ( velocity_x_northEast - velocity_x_southEast - velocity_x_northWest + velocity_x_southWest
- ) * hxInverse * hyInverse
+ ) * hxInverse * hyInverse * 4
+ ( velocity_y_north - 2 * velocity_y_center + velocity_y_south
- ) * hySquareInverse
+ ) * hySquareInverse * 4
) * this->dynamicalViscosity;
}
@@ -315,23 +315,23 @@ class LaxFridrichsMomentumY< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real
- ( rho_v[ down ] * velocity_z_down ) ) * hzInverse )
// T_12_y
- ( ( velocity_x_northEast - velocity_x_southEast - velocity_x_northWest + velocity_x_southWest
- ) * hxInverse * hyInverse
+ ) * hxInverse * hyInverse * 4
+ ( velocity_y_north - 2 * velocity_y_center + velocity_y_south
- ) * hySquareInverse
+ ) * hySquareInverse * 4
) * this->dynamicalViscosity
// 3D T_22_y
- ( 4.0 / 3.0 * ( velocity_y_north - 2 * velocity_y_center + velocity_y_south
- ) * hySquareInverse
+ ) * hySquareInverse * 4
- 2.0 / 3.0 * ( velocity_x_northEast - velocity_x_southEast - velocity_x_northWest + velocity_x_southWest
- ) * hxInverse * hyInverse
+ ) * hxInverse * hyInverse * 4
- 2.0 / 3.0 * ( velocity_z_upNorth - velocity_z_downNorth - velocity_z_upSouth + velocity_z_downSouth
- ) * hyInverse * hzInverse
+ ) * hyInverse * hzInverse * 4
) * this->dynamicalViscosity
// T_32_y
- ( ( velocity_z_upNorth - velocity_z_downNorth - velocity_z_upSouth + velocity_z_downSouth
- ) * hyInverse * hzInverse
+ ) * hyInverse * hzInverse * 4
+ ( velocity_y_north - 2 * velocity_y_center + velocity_y_south
- ) * hySquareInverse
+ ) * hySquareInverse * 4
) * this->dynamicalViscosity;
}
diff --git a/examples/flow/LaxFridrichsMomentumZ.h b/examples/flow/LaxFridrichsMomentumZ.h
index 0ae1e6f357cf74d6a43a7df4b79a962077063538..be3bae5c3f4cb703fc56cda2b4cdad9fed2e2c89 100644
--- a/examples/flow/LaxFridrichsMomentumZ.h
+++ b/examples/flow/LaxFridrichsMomentumZ.h
@@ -266,17 +266,17 @@ class LaxFridrichsMomentumZ< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real
+ ( ( rho_w[ up ] * velocity_z_up + pressure_up )
- ( rho_w[ down ] * velocity_z_down + pressure_down ) )* hzInverse )
// T_13_z
- - ( ( velocity_z_up - 2 * velocity_z_center + velocity_z_down ) * hzSquareInverse
- + ( velocity_x_upEast - velocity_x_downEast - velocity_x_upWest + velocity_x_downWest ) * hxInverse * hzInverse )
+ - ( ( velocity_z_up - 2 * velocity_z_center + velocity_z_down ) * hzSquareInverse * 4
+ + ( velocity_x_upEast - velocity_x_downEast - velocity_x_upWest + velocity_x_downWest ) * hxInverse * hzInverse * 4 )
* this->dynamicalViscosity
// T_23_z
- - ( ( velocity_y_upNorth - velocity_y_downNorth - velocity_y_upSouth + velocity_y_downSouth ) * hyInverse * hzInverse
- + ( velocity_z_up - 2 * velocity_z_center + velocity_z_down ) * hzSquareInverse )
+ - ( ( velocity_y_upNorth - velocity_y_downNorth - velocity_y_upSouth + velocity_y_downSouth ) * hyInverse * hzInverse * 4
+ + ( velocity_z_up - 2 * velocity_z_center + velocity_z_down ) * hzSquareInverse * 4 )
* this->dynamicalViscosity
// 3D T_33_z
- - ( 4.0 / 3.0 * ( velocity_z_up - 2 * velocity_z_center + velocity_z_down ) * hzSquareInverse
- - 2.0 / 3.0 * ( velocity_y_upNorth - velocity_y_downNorth - velocity_y_upSouth + velocity_y_downSouth ) * hyInverse * hzInverse
- - 2.0 / 3.0 * ( velocity_x_upEast - velocity_x_downEast - velocity_x_upWest + velocity_x_downWest ) * hxInverse * hzInverse )
+ - ( 4.0 / 3.0 * ( velocity_z_up - 2 * velocity_z_center + velocity_z_down ) * hzSquareInverse * 4
+ - 2.0 / 3.0 * ( velocity_y_upNorth - velocity_y_downNorth - velocity_y_upSouth + velocity_y_downSouth ) * hyInverse * hzInverse * 4
+ - 2.0 / 3.0 * ( velocity_x_upEast - velocity_x_downEast - velocity_x_upWest + velocity_x_downWest ) * hxInverse * hzInverse * 4 )
* this->dynamicalViscosity;
}
diff --git a/examples/flow/PhysicalVariablesGetter.h b/examples/flow/PhysicalVariablesGetter.h
index 58fd5df7a5b467a85506aff1c2b28841443269f3..f1ba6bd1222b8653faeaac041606c101a071e188 100644
--- a/examples/flow/PhysicalVariablesGetter.h
+++ b/examples/flow/PhysicalVariablesGetter.h
@@ -50,8 +50,11 @@ class PhysicalVariablesGetter
RealType operator()( const EntityType& meshEntity,
const RealType& time = 0.0 ) const
{
- return momentum.template getData< DeviceType >()( meshEntity ) /
- density.template getData< DeviceType >()( meshEntity );
+ if( density.template getData< DeviceType >()( meshEntity ) == 0.0 )
+ return 0;
+ else
+ return momentum.template getData< DeviceType >()( meshEntity ) /
+ density.template getData< DeviceType >()( meshEntity );
}
protected:
@@ -77,7 +80,10 @@ class PhysicalVariablesGetter
const RealType e = energy.template getData< DeviceType >()( meshEntity );
const RealType rho = density.template getData< DeviceType >()( meshEntity );
const RealType momentumNorm = momentum.template getData< DeviceType >().getVector( meshEntity ).lpNorm( 2.0 );
- return ( gamma - 1.0 ) * ( e - 0.5 * momentumNorm * momentumNorm / rho );
+ if( rho == 0.0 )
+ return 0;
+ else
+ return ( gamma - 1.0 ) * ( e - 0.5 * momentumNorm * momentumNorm / rho );
}
protected:
diff --git a/examples/flow/RiemannProblemInitialCondition.h b/examples/flow/RiemannProblemInitialCondition.h
index 664d93ea326250f47218704868003b8128084845..a712934143c551899e2e005a98e98eea6f1adf7b 100644
--- a/examples/flow/RiemannProblemInitialCondition.h
+++ b/examples/flow/RiemannProblemInitialCondition.h
@@ -887,7 +887,7 @@ class RiemannProblemInitialCondition
this->NWDDensity = parameters.getParameter< RealType >( prefix + "NWD-density" );
this->NWDVelocity.setup( parameters, prefix + "NWD-velocity-" );
this->NWDPressure = parameters.getParameter< RealType >( prefix + "NWD-pressure" );
- this->SWUEnergy = Energy( NWDDensity, NWDPressure, gamma, NWDVelocity);
+ this->NWDEnergy = Energy( NWDDensity, NWDPressure, gamma, NWDVelocity);
this->NWDMomentum = NWDVelocity * NWDDensity;
this->SWDDensity = parameters.getParameter< RealType >( prefix + "SWD-density" );
@@ -1257,11 +1257,7 @@ class RiemannProblemInitialCondition
this->SEDVelocity = PointLoad(preSEDVelocityX, preSEDVelocityY, preSEDVelocityZ);
this->SEDPressure = preSEDPressure;
this->SEDEnergy = Energy( SEDDensity, SEDPressure, gamma, SEDVelocity);
- this->SEDMomentum = SEDVelocity * SEDDensity;
-
- std::cout << this->SEDEnergy;
- std::cout << this->SWDEnergy;
-
+ this->SEDMomentum = SEDVelocity * SEDDensity;
}
PointType PointLoad( RealType ValueX, RealType ValueY, RealType ValueZ)
diff --git a/examples/flow/navierStokes.h b/examples/flow/navierStokes.h
index 38a7affbee41fb4de4c7ff0be3e6d4d8a89a4416..eb257e8cf07fc087180b89a55921bd460a61ab21 100644
--- a/examples/flow/navierStokes.h
+++ b/examples/flow/navierStokes.h
@@ -10,6 +10,7 @@
#include "navierStokesBuildConfigTag.h"
#include "RiemannProblemInitialCondition.h"
+#include "1DBoundaryConditions.h"
using namespace TNL;
@@ -35,6 +36,9 @@ template< typename ConfigTag >class navierStokesConfig
config.addEntryEnum< String >( "dirichlet" );
config.addEntryEnum< String >( "neumann" );
config.addEntry< double >( "boundary-conditions-constant", "This sets a value in case of the constant boundary conditions." );
+ config.addEntry< double >( "speed-increment", "This sets increment of input speed.", 0.0 );
+ config.addEntry< double >( "speed-increment-until", "This sets time until input speed will rose", -0.1 );
+ config.addEntry< double >( "cavity-speed", "This sets speed parameter of cavity", 0.0 );
typedef Meshes::Grid< 3 > Mesh;
LaxFridrichs< Mesh >::configSetup( config, "inviscid-operators-" );
RiemannProblemInitialCondition< Mesh >::configSetup( config );
@@ -72,6 +76,12 @@ class navierStokesSetter
* The following code is for the Dirichlet and the Neumann boundary conditions.
* Both can be constant or defined as descrete values of Vector.
*/
+
+ typedef Functions::Analytic::Constant< Dimension, Real > Constant;
+ typedef BoundaryConditions< MeshType, Constant, Real, Index > BoundaryConditions;
+ typedef navierStokesProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
+ SolverStarter solverStarter;
+ return solverStarter.template run< Problem >( parameters );/*
String boundaryConditionsType = parameters.getParameter< String >( "boundary-conditions-type" );
if( parameters.checkParameter( "boundary-conditions-constant" ) )
{
@@ -102,7 +112,7 @@ class navierStokesSetter
typedef navierStokesProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
SolverStarter solverStarter;
return solverStarter.template run< Problem >( parameters );
- }
+ }*/
return true;}
@@ -114,4 +124,4 @@ int main( int argc, char* argv[] )
if( ! solver. run( argc, argv ) )
return EXIT_FAILURE;
return EXIT_SUCCESS;
-}
+};
diff --git a/examples/flow/navierStokesProblem.h b/examples/flow/navierStokesProblem.h
index d587ffef8fc91d4c26f2d5807fc25b8b81518f79..ce4473d55746ecf7a3ed3f35f9b81b2fe4de5304 100644
--- a/examples/flow/navierStokesProblem.h
+++ b/examples/flow/navierStokesProblem.h
@@ -121,7 +121,10 @@ class navierStokesProblem:
VelocityFieldPointer velocity;
MeshFunctionPointer pressure;
- RealType gamma;
+ RealType gamma;
+ RealType speedIncrement;
+ RealType cavitySpeed;
+ RealType speedIncrementUntil;
};
} // namespace TNL
diff --git a/examples/flow/navierStokesProblem_impl.h b/examples/flow/navierStokesProblem_impl.h
index 478b4b206cefe7a6f0b2703c52633037c6eb3605..1da32e7c9cbb054ea497a8308009f1626a737977 100644
--- a/examples/flow/navierStokesProblem_impl.h
+++ b/examples/flow/navierStokesProblem_impl.h
@@ -28,6 +28,11 @@
#include "LaxFridrichsMomentumY.h"
#include "LaxFridrichsMomentumZ.h"
+#include "1DDensityBoundaryCondition.h"
+#include "1DMomentumXBoundaryCondition.h"
+#include "1DMomentumYBoundaryCondition.h"
+#include "1DEnergyBoundaryCondition.h"
+
namespace TNL {
template< typename Mesh,
@@ -127,6 +132,9 @@ setInitialCondition( const Config::ParameterContainer& parameters,
CompressibleConservativeVariables< MeshType > conservativeVariables;
conservativeVariables.bind( mesh, dofs );
const String& initialConditionType = parameters.getParameter< String >( "initial-condition" );
+ this->speedIncrementUntil = parameters.getParameter< RealType >( "speed-increment-until" );
+ this->speedIncrement = parameters.getParameter< RealType >( "speed-increment" );
+ this->cavitySpeed = parameters.getParameter< RealType >( "cavity-speed" );
if( initialConditionType == "riemann-problem" )
{
RiemannProblemInitialCondition< MeshType > initialCondition;
@@ -202,6 +210,9 @@ makeSnapshot( const RealType& time,
fileName.setFileNameBase( "energy-" );
if( ! this->conservativeVariables->getEnergy()->save( fileName.getFileName() ) )
return false;
+ fileName.setFileNameBase( "momentum-" );
+ if( ! this->conservativeVariables->getMomentum()->save( fileName.getFileName() ) )
+ return false;
return true;
}
@@ -249,23 +260,50 @@ getExplicitUpdate( const RealType& time,
this->inviscidOperatorsPointer->setVelocity( this->velocity );
this->inviscidOperatorsPointer->setPressure( this->pressure );
+ /****
+ * Set Up Boundary Conditions
+ */
+
+ typedef typename BoundaryCondition::DensityBoundaryConditionsType DensityBoundaryConditionsType;
+ typedef typename BoundaryCondition::MomentumXBoundaryConditionsType MomentumXBoundaryConditionsType;
+ typedef typename BoundaryCondition::MomentumYBoundaryConditionsType MomentumYBoundaryConditionsType;
+ typedef typename BoundaryCondition::EnergyBoundaryConditionsType EnergyBoundaryConditionsType;
+
+ /****
+ * Update Boundary Conditions
+ */
+ if(this->speedIncrementUntil > time )
+ {
+ this->boundaryConditionPointer->setTimestep(this->speedIncrement);
+ }
+ else
+ {
+ this->boundaryConditionPointer->setTimestep(0);
+ }
+ std::cout<<tau<< std::endl;
+ std::cout<<time<< std::endl;
+ getchar();
+ this->boundaryConditionPointer->setCavitySpeed(this->cavitySpeed);
+ this->boundaryConditionPointer->setCompressibleConservativeVariables(this->conservativeVariables);
+ this->boundaryConditionPointer->setGamma(this->gamma);
+ this->boundaryConditionPointer->setPressure(this->pressure);
+
/****
* Continuity equation
*/
- Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, ContinuityOperatorType, BoundaryCondition, RightHandSide > explicitUpdaterContinuity;
+ Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, ContinuityOperatorType, DensityBoundaryConditionsType, RightHandSide > explicitUpdaterContinuity;
explicitUpdaterContinuity.setDifferentialOperator( this->inviscidOperatorsPointer->getContinuityOperator() );
- explicitUpdaterContinuity.setBoundaryConditions( this->boundaryConditionPointer );
+ explicitUpdaterContinuity.setBoundaryConditions( this->boundaryConditionPointer->getDensityBoundaryCondition() );
explicitUpdaterContinuity.setRightHandSide( this->rightHandSidePointer );
explicitUpdaterContinuity.template update< typename Mesh::Cell >( time, tau, mesh,
this->conservativeVariables->getDensity(),
this->conservativeVariablesRHS->getDensity() );
-
/****
* Momentum equations
*/
- Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, MomentumXOperatorType, BoundaryCondition, RightHandSide > explicitUpdaterMomentumX;
+ Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, MomentumXOperatorType, MomentumXBoundaryConditionsType, RightHandSide > explicitUpdaterMomentumX;
explicitUpdaterMomentumX.setDifferentialOperator( this->inviscidOperatorsPointer->getMomentumXOperator() );
- explicitUpdaterMomentumX.setBoundaryConditions( this->boundaryConditionPointer );
+ explicitUpdaterMomentumX.setBoundaryConditions( this->boundaryConditionPointer->getMomentumXBoundaryCondition() );
explicitUpdaterMomentumX.setRightHandSide( this->rightHandSidePointer );
explicitUpdaterMomentumX.template update< typename Mesh::Cell >( time, tau, mesh,
( *this->conservativeVariables->getMomentum() )[ 0 ], // uRhoVelocityX,
@@ -273,9 +311,9 @@ getExplicitUpdate( const RealType& time,
if( Dimensions > 1 )
{
- Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, MomentumYOperatorType, BoundaryCondition, RightHandSide > explicitUpdaterMomentumY;
+ Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, MomentumYOperatorType, MomentumYBoundaryConditionsType, RightHandSide > explicitUpdaterMomentumY;
explicitUpdaterMomentumY.setDifferentialOperator( this->inviscidOperatorsPointer->getMomentumYOperator() );
- explicitUpdaterMomentumY.setBoundaryConditions( this->boundaryConditionPointer );
+ explicitUpdaterMomentumY.setBoundaryConditions( this->boundaryConditionPointer->getMomentumYBoundaryCondition() );
explicitUpdaterMomentumY.setRightHandSide( this->rightHandSidePointer );
explicitUpdaterMomentumY.template update< typename Mesh::Cell >( time, tau, mesh,
( *this->conservativeVariables->getMomentum() )[ 1 ], // uRhoVelocityX,
@@ -284,22 +322,21 @@ getExplicitUpdate( const RealType& time,
if( Dimensions > 2 )
{
- Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, MomentumZOperatorType, BoundaryCondition, RightHandSide > explicitUpdaterMomentumZ;
+ Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, MomentumZOperatorType, MomentumXBoundaryConditionsType, RightHandSide > explicitUpdaterMomentumZ;
explicitUpdaterMomentumZ.setDifferentialOperator( this->inviscidOperatorsPointer->getMomentumZOperator() );
- explicitUpdaterMomentumZ.setBoundaryConditions( this->boundaryConditionPointer );
+ explicitUpdaterMomentumZ.setBoundaryConditions( this->boundaryConditionPointer->getMomentumXBoundaryCondition() );
explicitUpdaterMomentumZ.setRightHandSide( this->rightHandSidePointer );
explicitUpdaterMomentumZ.template update< typename Mesh::Cell >( time, tau, mesh,
( *this->conservativeVariables->getMomentum() )[ 2 ], // uRhoVelocityX,
( *this->conservativeVariablesRHS->getMomentum() )[ 2 ] ); //, fuRhoVelocityX );
}
-
/****
* Energy equation
*/
- Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, EnergyOperatorType, BoundaryCondition, RightHandSide > explicitUpdaterEnergy;
+ Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, EnergyOperatorType, EnergyBoundaryConditionsType, RightHandSide > explicitUpdaterEnergy;
explicitUpdaterEnergy.setDifferentialOperator( this->inviscidOperatorsPointer->getEnergyOperator() );
- explicitUpdaterEnergy.setBoundaryConditions( this->boundaryConditionPointer );
+ explicitUpdaterEnergy.setBoundaryConditions( this->boundaryConditionPointer->getEnergyBoundaryCondition() );
explicitUpdaterEnergy.setRightHandSide( this->rightHandSidePointer );
explicitUpdaterEnergy.template update< typename Mesh::Cell >( time, tau, mesh,
this->conservativeVariables->getEnergy(), // uRhoVelocityX,
diff --git a/examples/flow/run-euler b/examples/flow/run-navier-stokes
similarity index 100%
rename from examples/flow/run-euler
rename to examples/flow/run-navier-stokes
diff --git a/examples/inviscid-flow/LaxFridrichsContinuity.h b/examples/inviscid-flow/LaxFridrichsContinuity.h
index c9980874656bdcb7f8c1eab04286995dcfea61c8..45ad4d52b12d402365a40cac043d5525e230cecb 100644
--- a/examples/inviscid-flow/LaxFridrichsContinuity.h
+++ b/examples/inviscid-flow/LaxFridrichsContinuity.h
@@ -119,11 +119,8 @@ class LaxFridrichsContinuity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Re
const IndexType& west = neighborEntities.template getEntityIndex< -1 >();
const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
const RealType& velocity_x_east = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ];
-// return 1.0 / ( 2.0 * this->tau ) * this->artificialViscosity * ( u[ west ] - 2.0 * u[ center ] + u[ east ] )
-// - 0.5 * ( u[ west ] * velocity_x_west - u[ east ] * velocity_x_east ) * hxInverse;
return 1.0 / ( 2.0 * this->tau ) * this->artificialViscosity * ( u[ west ] - 2.0 * u[ center ] + u[ east ] )
- 0.5 * ( u[ east ] * velocity_x_east - u[ west ] * velocity_x_west ) * hxInverse;
-
}
/*template< typename MeshEntity >
diff --git a/examples/inviscid-flow/LaxFridrichsEnergy.h b/examples/inviscid-flow/LaxFridrichsEnergy.h
index 3549e85f35265f68816f58cc6dd65e1a0db45e1c..18c824762b8c677253dbd4e494be7ad3aea7e769 100644
--- a/examples/inviscid-flow/LaxFridrichsEnergy.h
+++ b/examples/inviscid-flow/LaxFridrichsEnergy.h
@@ -125,7 +125,7 @@ class LaxFridrichsEnergy< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real,
const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
return 1.0 / ( 2.0 * this->tau ) * this->artificialViscosity * ( e[ west ] - 2.0 * e[ center ] + e[ east ] )
- 0.5 * ( ( e[ east ] + pressure_east ) * velocity_x_east
- - ( e[ west ] + pressure_west ) * velocity_x_west ) * hxInverse;
+ - ( e[ west ] + pressure_west ) * velocity_x_west ) * hxInverse;
}
diff --git a/examples/inviscid-flow/PhysicalVariablesGetter.h b/examples/inviscid-flow/PhysicalVariablesGetter.h
index 58fd5df7a5b467a85506aff1c2b28841443269f3..f1ba6bd1222b8653faeaac041606c101a071e188 100644
--- a/examples/inviscid-flow/PhysicalVariablesGetter.h
+++ b/examples/inviscid-flow/PhysicalVariablesGetter.h
@@ -50,8 +50,11 @@ class PhysicalVariablesGetter
RealType operator()( const EntityType& meshEntity,
const RealType& time = 0.0 ) const
{
- return momentum.template getData< DeviceType >()( meshEntity ) /
- density.template getData< DeviceType >()( meshEntity );
+ if( density.template getData< DeviceType >()( meshEntity ) == 0.0 )
+ return 0;
+ else
+ return momentum.template getData< DeviceType >()( meshEntity ) /
+ density.template getData< DeviceType >()( meshEntity );
}
protected:
@@ -77,7 +80,10 @@ class PhysicalVariablesGetter
const RealType e = energy.template getData< DeviceType >()( meshEntity );
const RealType rho = density.template getData< DeviceType >()( meshEntity );
const RealType momentumNorm = momentum.template getData< DeviceType >().getVector( meshEntity ).lpNorm( 2.0 );
- return ( gamma - 1.0 ) * ( e - 0.5 * momentumNorm * momentumNorm / rho );
+ if( rho == 0.0 )
+ return 0;
+ else
+ return ( gamma - 1.0 ) * ( e - 0.5 * momentumNorm * momentumNorm / rho );
}
protected:
diff --git a/examples/inviscid-flow/RiemannProblemInitialCondition.h b/examples/inviscid-flow/RiemannProblemInitialCondition.h
index 664d93ea326250f47218704868003b8128084845..a036747f255fd1477c02a55366d2f57ef2f7f990 100644
--- a/examples/inviscid-flow/RiemannProblemInitialCondition.h
+++ b/examples/inviscid-flow/RiemannProblemInitialCondition.h
@@ -119,7 +119,7 @@ class RiemannProblemInitialConditionSetter< Meshes::Grid< 1,MeshReal, Device, Me
typedef typename MeshType::CoordinatesType CoordinatesType;
MeshType mesh = (* conservativeVariables.getDensity()).getMesh();
for( int i = 0; i < mesh.getDimensions().x(); i++)
- if ( i <= this->discontinuityPlacement[ 0 ] * mesh.getDimensions().x() )
+ if ( i < this->discontinuityPlacement[ 0 ] * mesh.getDimensions().x() )
{
CellType cell(mesh, CoordinatesType(i));
cell.refresh();
@@ -139,7 +139,7 @@ class RiemannProblemInitialConditionSetter< Meshes::Grid< 1,MeshReal, Device, Me
typedef typename MeshType::CoordinatesType CoordinatesType;
MeshType mesh = (* conservativeVariables.getDensity()).getMesh();
for( int i = 0; i < mesh.getDimensions().x(); i++)
- if ( i <= this->discontinuityPlacement[ 0 ] * mesh.getDimensions().x() )
+ if ( i < this->discontinuityPlacement[ 0 ] * mesh.getDimensions().x() )
{
CellType cell(mesh, CoordinatesType(i));
cell.refresh();
@@ -159,7 +159,7 @@ class RiemannProblemInitialConditionSetter< Meshes::Grid< 1,MeshReal, Device, Me
typedef typename MeshType::CoordinatesType CoordinatesType;
MeshType mesh = (* conservativeVariables.getDensity()).getMesh();
for( int i = 0; i < mesh.getDimensions().x(); i++)
- if ( i <= this->discontinuityPlacement[ 0 ] * mesh.getDimensions().x() )
+ if ( i < this->discontinuityPlacement[ 0 ] * mesh.getDimensions().x() )
{
CellType cell(mesh, CoordinatesType(i));
cell.refresh();