From 7453b698ba8b54d803a58b8ba7c97525d0c5b805 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jan=20Sch=C3=A4fer?= <schafjan@fjfi.cvut.cz>
Date: Wed, 18 Apr 2018 21:31:11 +0200
Subject: [PATCH] navier stokes + tenzory 2D fixed
---
...onditions.h => BoundaryConditionsBoiler.h} | 35 +-
examples/flow/BoundaryConditionsCavity.h | 137 +++
.../flow/DensityBoundaryConditionBoiler.h | 542 +++++++++++
...ion.h => DensityBoundaryConditionCavity.h} | 50 +-
examples/flow/EnergyBoundaryConditionBoiler.h | 854 ++++++++++++++++++
examples/flow/EnergyBoundaryConditionCavity.h | 673 ++++++++++++++
examples/flow/LaxFridrichsEnergy.h | 24 +-
examples/flow/LaxFridrichsMomentumX.h | 18 +-
examples/flow/LaxFridrichsMomentumY.h | 18 +-
.../flow/MomentumXBoundaryConditionBoiler.h | 594 ++++++++++++
...n.h => MomentumXBoundaryConditionCavity.h} | 67 +-
.../flow/MomentumYBoundaryConditionBoiler.h | 588 ++++++++++++
...n.h => MomentumYBoundaryConditionCavity.h} | 102 +--
.../flow/MomentumZBoundaryConditionBoiler.h | 563 ++++++++++++
...n.h => MomentumZBoundaryConditionCavity.h} | 50 +-
examples/flow/navierStokes.h | 47 +-
examples/flow/navierStokesProblem_impl.h | 55 +-
.../tnl-transport-equation.h | 9 +
src/TNL/Operators/Advection/CMakeLists.txt | 3 +-
19 files changed, 4186 insertions(+), 243 deletions(-)
rename examples/flow/{1DBoundaryConditions.h => BoundaryConditionsBoiler.h} (72%)
create mode 100644 examples/flow/BoundaryConditionsCavity.h
create mode 100644 examples/flow/DensityBoundaryConditionBoiler.h
rename examples/flow/{1DDensityBoundaryCondition.h => DensityBoundaryConditionCavity.h} (90%)
create mode 100644 examples/flow/EnergyBoundaryConditionBoiler.h
create mode 100644 examples/flow/EnergyBoundaryConditionCavity.h
create mode 100644 examples/flow/MomentumXBoundaryConditionBoiler.h
rename examples/flow/{1DMomentumXBoundaryCondition.h => MomentumXBoundaryConditionCavity.h} (90%)
create mode 100644 examples/flow/MomentumYBoundaryConditionBoiler.h
rename examples/flow/{1DEnergyBoundaryCondition.h => MomentumYBoundaryConditionCavity.h} (83%)
create mode 100644 examples/flow/MomentumZBoundaryConditionBoiler.h
rename examples/flow/{1DMomentumYBoundaryCondition.h => MomentumZBoundaryConditionCavity.h} (90%)
diff --git a/examples/flow/1DBoundaryConditions.h b/examples/flow/BoundaryConditionsBoiler.h
similarity index 72%
rename from examples/flow/1DBoundaryConditions.h
rename to examples/flow/BoundaryConditionsBoiler.h
index 0a07664dcb..0cba68d7fa 100644
--- a/examples/flow/1DBoundaryConditions.h
+++ b/examples/flow/BoundaryConditionsBoiler.h
@@ -1,9 +1,10 @@
#include <TNL/Functions/FunctionAdapter.h>
-#include "1DDensityBoundaryCondition.h"
-#include "1DMomentumXBoundaryCondition.h"
-#include "1DMomentumYBoundaryCondition.h"
-#include "1DEnergyBoundaryCondition.h"
+#include "DensityBoundaryConditionBoiler.h"
+#include "MomentumXBoundaryConditionBoiler.h"
+#include "MomentumYBoundaryConditionBoiler.h"
+#include "MomentumZBoundaryConditionBoiler.h"
+#include "EnergyBoundaryConditionBoiler.h"
namespace TNL {
@@ -11,7 +12,7 @@ template< typename Mesh,
typename Function,
typename Real = typename Mesh::RealType,
typename Index = typename Mesh::IndexType >
-class BoundaryConditions
+class BoundaryConditionsBoiler
{
public:
typedef Mesh MeshType;
@@ -21,15 +22,17 @@ class BoundaryConditions
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 TNL::Operators::DensityBoundaryConditionsBoiler< MeshType, FunctionType, RealType, IndexType > DensityBoundaryConditionsType;
+ typedef TNL::Operators::MomentumXBoundaryConditionsBoiler< MeshType, FunctionType, RealType, IndexType > MomentumXBoundaryConditionsType;
+ typedef TNL::Operators::MomentumYBoundaryConditionsBoiler< MeshType, FunctionType, RealType, IndexType > MomentumYBoundaryConditionsType;
+ typedef TNL::Operators::MomentumZBoundaryConditionsBoiler< MeshType, FunctionType, RealType, IndexType > MomentumZBoundaryConditionsType;
+ typedef TNL::Operators::EnergyBoundaryConditionsBoiler< MeshType, FunctionType, RealType, IndexType > EnergyBoundaryConditionsType;
typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
typedef SharedPointer< DensityBoundaryConditionsType > DensityBoundaryConditionsTypePointer;
typedef SharedPointer< MomentumXBoundaryConditionsType > MomentumXBoundaryConditionsTypePointer;
typedef SharedPointer< MomentumYBoundaryConditionsType > MomentumYBoundaryConditionsTypePointer;
+ typedef SharedPointer< MomentumZBoundaryConditionsType > MomentumZBoundaryConditionsTypePointer;
typedef SharedPointer< EnergyBoundaryConditionsType > EnergyBoundaryConditionsTypePointer;
typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
typedef SharedPointer< MeshType > MeshPointer;
@@ -47,6 +50,7 @@ class BoundaryConditions
this->densityBoundaryConditionsPointer->setup( meshPointer, parameters, prefix);
this->momentumXBoundaryConditionsPointer->setup( meshPointer, parameters, prefix);
this->momentumYBoundaryConditionsPointer->setup( meshPointer, parameters, prefix);
+ this->momentumZBoundaryConditionsPointer->setup( meshPointer, parameters, prefix);
this->energyBoundaryConditionsPointer->setup( meshPointer, parameters, prefix);
return true;
}
@@ -56,6 +60,7 @@ class BoundaryConditions
this->densityBoundaryConditionsPointer->setCompressibleConservativeVariables(compressibleConservativeVariables);
this->momentumXBoundaryConditionsPointer->setCompressibleConservativeVariables(compressibleConservativeVariables);
this->momentumYBoundaryConditionsPointer->setCompressibleConservativeVariables(compressibleConservativeVariables);
+ this->momentumZBoundaryConditionsPointer->setCompressibleConservativeVariables(compressibleConservativeVariables);
this->energyBoundaryConditionsPointer->setCompressibleConservativeVariables(compressibleConservativeVariables);
}
@@ -64,6 +69,7 @@ class BoundaryConditions
this->densityBoundaryConditionsPointer->setTimestep(timestep);
this->momentumXBoundaryConditionsPointer->setTimestep(timestep);
this->momentumYBoundaryConditionsPointer->setTimestep(timestep);
+ this->momentumZBoundaryConditionsPointer->setTimestep(timestep);
this->energyBoundaryConditionsPointer->setTimestep(timestep);
}
@@ -72,6 +78,7 @@ class BoundaryConditions
this->densityBoundaryConditionsPointer->setGamma(gamma);
this->momentumXBoundaryConditionsPointer->setGamma(gamma);
this->momentumYBoundaryConditionsPointer->setGamma(gamma);
+ this->momentumZBoundaryConditionsPointer->setGamma(gamma);
this->energyBoundaryConditionsPointer->setGamma(gamma);
}
@@ -80,13 +87,15 @@ class BoundaryConditions
this->densityBoundaryConditionsPointer->setPressure(pressure);
this->momentumXBoundaryConditionsPointer->setPressure(pressure);
this->momentumYBoundaryConditionsPointer->setPressure(pressure);
+ this->momentumZBoundaryConditionsPointer->setPressure(pressure);
this->energyBoundaryConditionsPointer->setPressure(pressure);
}
- void setCavitySpeed(const RealType cavitySpeed)
+ void setSpeed(const RealType cavitySpeed)
{
this->momentumXBoundaryConditionsPointer->setCavitySpeed(cavitySpeed);
this->momentumYBoundaryConditionsPointer->setCavitySpeed(cavitySpeed);
+ this->momentumZBoundaryConditionsPointer->setCavitySpeed(cavitySpeed);
this->energyBoundaryConditionsPointer->setCavitySpeed(cavitySpeed);
}
@@ -105,6 +114,11 @@ class BoundaryConditions
return this->momentumYBoundaryConditionsPointer;
}
+ MomentumZBoundaryConditionsTypePointer& getMomentumZBoundaryCondition()
+ {
+ return this->momentumZBoundaryConditionsPointer;
+ }
+
EnergyBoundaryConditionsTypePointer& getEnergyBoundaryCondition()
{
return this->energyBoundaryConditionsPointer;
@@ -115,6 +129,7 @@ class BoundaryConditions
DensityBoundaryConditionsTypePointer densityBoundaryConditionsPointer;
MomentumXBoundaryConditionsTypePointer momentumXBoundaryConditionsPointer;
MomentumYBoundaryConditionsTypePointer momentumYBoundaryConditionsPointer;
+ MomentumZBoundaryConditionsTypePointer momentumZBoundaryConditionsPointer;
EnergyBoundaryConditionsTypePointer energyBoundaryConditionsPointer;
};
diff --git a/examples/flow/BoundaryConditionsCavity.h b/examples/flow/BoundaryConditionsCavity.h
new file mode 100644
index 0000000000..8a42faea17
--- /dev/null
+++ b/examples/flow/BoundaryConditionsCavity.h
@@ -0,0 +1,137 @@
+#include <TNL/Functions/FunctionAdapter.h>
+
+#include "DensityBoundaryConditionCavity.h"
+#include "MomentumXBoundaryConditionCavity.h"
+#include "MomentumYBoundaryConditionCavity.h"
+#include "MomentumZBoundaryConditionCavity.h"
+#include "EnergyBoundaryConditionCavity.h"
+
+namespace TNL {
+
+template< typename Mesh,
+ typename Function,
+ typename Real = typename Mesh::RealType,
+ typename Index = typename Mesh::IndexType >
+class BoundaryConditionsCavity
+{
+ 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::DensityBoundaryConditionsCavity< MeshType, FunctionType, RealType, IndexType > DensityBoundaryConditionsType;
+ typedef TNL::Operators::MomentumXBoundaryConditionsCavity< MeshType, FunctionType, RealType, IndexType > MomentumXBoundaryConditionsType;
+ typedef TNL::Operators::MomentumYBoundaryConditionsCavity< MeshType, FunctionType, RealType, IndexType > MomentumYBoundaryConditionsType;
+ typedef TNL::Operators::MomentumZBoundaryConditionsCavity< MeshType, FunctionType, RealType, IndexType > MomentumZBoundaryConditionsType;
+ typedef TNL::Operators::EnergyBoundaryConditionsCavity< MeshType, FunctionType, RealType, IndexType > EnergyBoundaryConditionsType;
+ typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
+
+ typedef SharedPointer< DensityBoundaryConditionsType > DensityBoundaryConditionsTypePointer;
+ typedef SharedPointer< MomentumXBoundaryConditionsType > MomentumXBoundaryConditionsTypePointer;
+ typedef SharedPointer< MomentumYBoundaryConditionsType > MomentumYBoundaryConditionsTypePointer;
+ typedef SharedPointer< MomentumZBoundaryConditionsType > MomentumZBoundaryConditionsTypePointer;
+ 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->momentumZBoundaryConditionsPointer->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->momentumZBoundaryConditionsPointer->setCompressibleConservativeVariables(compressibleConservativeVariables);
+ this->energyBoundaryConditionsPointer->setCompressibleConservativeVariables(compressibleConservativeVariables);
+ }
+
+ void setTimestep(const RealType timestep)
+ {
+ this->densityBoundaryConditionsPointer->setTimestep(timestep);
+ this->momentumXBoundaryConditionsPointer->setTimestep(timestep);
+ this->momentumYBoundaryConditionsPointer->setTimestep(timestep);
+ this->momentumZBoundaryConditionsPointer->setTimestep(timestep);
+ this->energyBoundaryConditionsPointer->setTimestep(timestep);
+ }
+
+ void setGamma(const RealType gamma)
+ {
+ this->densityBoundaryConditionsPointer->setGamma(gamma);
+ this->momentumXBoundaryConditionsPointer->setGamma(gamma);
+ this->momentumYBoundaryConditionsPointer->setGamma(gamma);
+ this->momentumZBoundaryConditionsPointer->setGamma(gamma);
+ this->energyBoundaryConditionsPointer->setGamma(gamma);
+ }
+
+ void setPressure(const MeshFunctionPointer& pressure)
+ {
+ this->densityBoundaryConditionsPointer->setPressure(pressure);
+ this->momentumXBoundaryConditionsPointer->setPressure(pressure);
+ this->momentumYBoundaryConditionsPointer->setPressure(pressure);
+ this->momentumZBoundaryConditionsPointer->setPressure(pressure);
+ this->energyBoundaryConditionsPointer->setPressure(pressure);
+ }
+
+ void setSpeed(const RealType cavitySpeed)
+ {
+ this->momentumXBoundaryConditionsPointer->setCavitySpeed(cavitySpeed);
+ this->momentumYBoundaryConditionsPointer->setCavitySpeed(cavitySpeed);
+ this->momentumZBoundaryConditionsPointer->setCavitySpeed(cavitySpeed);
+ this->energyBoundaryConditionsPointer->setCavitySpeed(cavitySpeed);
+ }
+
+ DensityBoundaryConditionsTypePointer& getDensityBoundaryCondition()
+ {
+ return this->densityBoundaryConditionsPointer;
+ }
+
+ MomentumXBoundaryConditionsTypePointer& getMomentumXBoundaryCondition()
+ {
+ return this->momentumXBoundaryConditionsPointer;
+ }
+
+ MomentumYBoundaryConditionsTypePointer& getMomentumYBoundaryCondition()
+ {
+ return this->momentumYBoundaryConditionsPointer;
+ }
+
+ MomentumZBoundaryConditionsTypePointer& getMomentumZBoundaryCondition()
+ {
+ return this->momentumZBoundaryConditionsPointer;
+ }
+
+ EnergyBoundaryConditionsTypePointer& getEnergyBoundaryCondition()
+ {
+ return this->energyBoundaryConditionsPointer;
+ }
+
+
+ protected:
+ DensityBoundaryConditionsTypePointer densityBoundaryConditionsPointer;
+ MomentumXBoundaryConditionsTypePointer momentumXBoundaryConditionsPointer;
+ MomentumYBoundaryConditionsTypePointer momentumYBoundaryConditionsPointer;
+ MomentumZBoundaryConditionsTypePointer momentumZBoundaryConditionsPointer;
+ EnergyBoundaryConditionsTypePointer energyBoundaryConditionsPointer;
+
+};
+
+} //namespace TNL
diff --git a/examples/flow/DensityBoundaryConditionBoiler.h b/examples/flow/DensityBoundaryConditionBoiler.h
new file mode 100644
index 0000000000..c3bae7e3d9
--- /dev/null
+++ b/examples/flow/DensityBoundaryConditionBoiler.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 DensityBoundaryConditionsBoiler
+{
+
+};
+
+/****
+ * Base
+ */
+template< typename Function >
+class DensityBoundaryConditionsBoilerBase
+{
+ 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 DensityBoundaryConditionsBoiler< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public DensityBoundaryConditionsBoilerBase< 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 DensityBoundaryConditionsBoiler< MeshType, Function, Real, Index > ThisType;
+ typedef DensityBoundaryConditionsBoilerBase< 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 DensityBoundaryConditionsBoiler< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public DensityBoundaryConditionsBoilerBase< 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 DensityBoundaryConditionsBoiler< MeshType, Function, Real, Index > ThisType;
+ typedef DensityBoundaryConditionsBoilerBase< 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 )
+ {
+ if (entity.getCoordinates().y() < 0.8 * ( entity.getMesh().getDimensions().y() - 1 ) && false)
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ else
+ 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, -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, 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 DensityBoundaryConditionsBoiler< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public DensityBoundaryConditionsBoilerBase< 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 DensityBoundaryConditionsBoiler< MeshType, Function, Real, Index > ThisType;
+ typedef DensityBoundaryConditionsBoilerBase< 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, 0 >() ];
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ if (entity.getCoordinates().z() < 0.8 * ( entity.getMesh().getDimensions().z() - 1 ) )
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ else
+ return u[ neighborEntities.template getEntityIndex< -1, 0, 0 >() ];
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ if( entity.getCoordinates().z() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ // 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, 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, 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 DensityBoundaryConditionsBoiler< Mesh, Function, Real, Index >& bc )
+{
+ str << "Neumann boundary ConditionBoilers: function = " << bc.getFunction();
+ return str;
+}
+
+} // namespace Operators
+} // namespace TNL
+
diff --git a/examples/flow/1DDensityBoundaryCondition.h b/examples/flow/DensityBoundaryConditionCavity.h
similarity index 90%
rename from examples/flow/1DDensityBoundaryCondition.h
rename to examples/flow/DensityBoundaryConditionCavity.h
index b45afdd6d3..a2d34ce540 100644
--- a/examples/flow/1DDensityBoundaryCondition.h
+++ b/examples/flow/DensityBoundaryConditionCavity.h
@@ -20,7 +20,7 @@ template< typename Mesh,
typename Function,
typename Real = typename Mesh::RealType,
typename Index = typename Mesh::GlobalIndexType >
-class DensityBoundaryConditions
+class DensityBoundaryConditionsCavity
{
};
@@ -29,7 +29,7 @@ class DensityBoundaryConditions
* Base
*/
template< typename Function >
-class DensityBoundaryConditionsBase
+class DensityBoundaryConditionsCavityBase
{
public:
@@ -91,8 +91,8 @@ template< typename MeshReal,
typename Function,
typename Real,
typename Index >
-class DensityBoundaryConditions< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >
- : public DensityBoundaryConditionsBase< Function >,
+class DensityBoundaryConditionsCavity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public DensityBoundaryConditionsCavityBase< Function >,
public Operator< Meshes::Grid< 1, MeshReal, Device, MeshIndex >,
Functions::MeshBoundaryDomain,
1, 1,
@@ -111,8 +111,8 @@ class DensityBoundaryConditions< Meshes::Grid< 1, MeshReal, Device, MeshIndex >,
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 DensityBoundaryConditionsCavity< MeshType, Function, Real, Index > ThisType;
+ typedef DensityBoundaryConditionsCavityBase< Function > BaseType;
typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
@@ -213,8 +213,8 @@ template< typename MeshReal,
typename Function,
typename Real,
typename Index >
-class DensityBoundaryConditions< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >
- : public DensityBoundaryConditionsBase< Function >,
+class DensityBoundaryConditionsCavity< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public DensityBoundaryConditionsCavityBase< Function >,
public Operator< Meshes::Grid< 2, MeshReal, Device, MeshIndex >,
Functions::MeshBoundaryDomain,
2, 2,
@@ -234,8 +234,8 @@ class DensityBoundaryConditions< Meshes::Grid< 2, MeshReal, Device, MeshIndex >,
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 DensityBoundaryConditionsCavity< MeshType, Function, Real, Index > ThisType;
+ typedef DensityBoundaryConditionsCavityBase< Function > BaseType;
typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
@@ -361,8 +361,8 @@ template< typename MeshReal,
typename Function,
typename Real,
typename Index >
-class DensityBoundaryConditions< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >
- : public DensityBoundaryConditionsBase< Function >,
+class DensityBoundaryConditionsCavity< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public DensityBoundaryConditionsCavityBase< Function >,
public Operator< Meshes::Grid< 3, MeshReal, Device, MeshIndex >,
Functions::MeshBoundaryDomain,
3, 3,
@@ -381,8 +381,8 @@ class DensityBoundaryConditions< Meshes::Grid< 3, MeshReal, Device, MeshIndex >,
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 DensityBoundaryConditionsCavity< MeshType, Function, Real, Index > ThisType;
+ typedef DensityBoundaryConditionsCavityBase< Function > BaseType;
typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
@@ -399,34 +399,28 @@ class DensityBoundaryConditions< Meshes::Grid< 3, MeshReal, Device, MeshIndex >,
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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
// 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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
}
@@ -531,9 +525,9 @@ template< typename Mesh,
typename Function,
typename Real,
typename Index >
-std::ostream& operator << ( std::ostream& str, const DensityBoundaryConditions< Mesh, Function, Real, Index >& bc )
+std::ostream& operator << ( std::ostream& str, const DensityBoundaryConditionsCavity< Mesh, Function, Real, Index >& bc )
{
- str << "Neumann boundary conditions: function = " << bc.getFunction();
+ str << "Neumann boundary ConditionsCavity: function = " << bc.getFunction();
return str;
}
diff --git a/examples/flow/EnergyBoundaryConditionBoiler.h b/examples/flow/EnergyBoundaryConditionBoiler.h
new file mode 100644
index 0000000000..fe227d68f8
--- /dev/null
+++ b/examples/flow/EnergyBoundaryConditionBoiler.h
@@ -0,0 +1,854 @@
+/***************************************************************************
+ 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 EnergyBoundaryConditionsBoiler
+{
+
+};
+
+/****
+ * Base
+ */
+template< typename Function>
+class EnergyBoundaryConditionsBoilerBase
+{
+ 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 EnergyBoundaryConditionsBoiler< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public EnergyBoundaryConditionsBoilerBase< 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 EnergyBoundaryConditionsBoiler< MeshType, Function, Real, Index > ThisType;
+ typedef EnergyBoundaryConditionsBoilerBase< 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 EnergyBoundaryConditionsBoiler< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public EnergyBoundaryConditionsBoilerBase< 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 EnergyBoundaryConditionsBoiler< MeshType, Function, Real, Index > ThisType;
+ typedef EnergyBoundaryConditionsBoilerBase< 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 )
+ {
+ if( ( entity.getCoordinates().y() < 0.20 * ( entity.getMesh().getDimensions().y() - 1 ) ) && ( entity.getCoordinates().y() > 0.19 * ( 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
+ )
+ );
+ else
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ /*( (* this->pressure)[ neighborEntities.template getEntityIndex< 0, 0 >() ]
+ / ( this->gamma - 1 )
+ )
+ + 0.5
+ * (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ * (
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< 1, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 1, 0 >()]
+ + 0
+ )
+ *
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< 1, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 1, 0 >()]
+ + 0
+ )
+ );*/
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ if( ( entity.getCoordinates().y() < 0.20 * ( entity.getMesh().getDimensions().y() - 1 ) ) && ( entity.getCoordinates().y() > 0.19 * ( 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 * ( -1 )
+ *
+ this->cavitySpeed * ( -1 )
+ +
+ ( (* (* 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
+ )
+ );
+ else if( entity.getCoordinates().y() > 0.8 * ( entity.getMesh().getDimensions().y() - 1 ) && false )
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ else
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ /*( (* this->pressure)[ neighborEntities.template getEntityIndex< 0, 0 >() ]
+ / ( this->gamma - 1 )
+ )
+ + 0.5
+ * (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ * (
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< -1, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< -1, 0 >()]
+ + 0
+ )
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< -1, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< -1, 0 >()]
+ + 0
+ )
+ );*/
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ if( ( entity.getCoordinates().x() < 0.6 * ( entity.getMesh().getDimensions().x() - 1 ) ) && ( entity.getCoordinates().x() > 0.4 * ( entity.getMesh().getDimensions().x() - 1 ) ) )
+ return ( (* this->pressure)[ neighborEntities.template getEntityIndex< 0, 0 >() ]
+ / ( this->gamma - 1 )
+ )
+ + 0.5
+ * (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ * (
+ (
+ (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ + 0
+ )
+ * (
+ (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ + 0
+ )
+ +
+ this->cavitySpeed
+ *
+ this->cavitySpeed
+ );
+ else
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ /*( (* this->pressure)[ neighborEntities.template getEntityIndex< 0, 0 >() ]
+ / ( this->gamma - 1 )
+ )
+ + 0.5
+ * (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ * (
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0, 1 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 1 >()]
+ + 0
+ )
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0, 1 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 1 >()]
+ + 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, -1 >() ];
+ /*return ( (* this->pressure)[ neighborEntities.template getEntityIndex< 0, 0 >() ]
+ / ( this->gamma - 1 )
+ )
+ + 0.5
+ * (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ * (
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0, -1 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, -1 >()]
+ + 0
+ )
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0, -1 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, -1 >()]
+ + 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 EnergyBoundaryConditionsBoiler< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public EnergyBoundaryConditionsBoilerBase< 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 EnergyBoundaryConditionsBoiler< MeshType, Function, Real, Index > ThisType;
+ typedef EnergyBoundaryConditionsBoilerBase< 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 )
+ {
+ if( ( entity.getCoordinates().y() < 0.59 * ( entity.getMesh().getDimensions().y() - 1 ) ) && ( entity.getCoordinates().y() > 0.39 * ( entity.getMesh().getDimensions().y() - 1 ) )
+ &&( entity.getCoordinates().z() < 0.59 * ( entity.getMesh().getDimensions().z() - 1 ) ) && ( entity.getCoordinates().z() > 0.39 * ( entity.getMesh().getDimensions().z() - 1 ) ) )
+ return ( (* this->pressure)[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ]
+ / ( this->gamma - 1 )
+ )
+ + 0.5
+ * (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ * (
+ this->cavitySpeed
+ *
+ this->cavitySpeed
+ +
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ +
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 2 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 2 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ );
+ else
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ if( ( entity.getCoordinates().y() < 0.59 * ( entity.getMesh().getDimensions().y() - 1 ) ) && ( entity.getCoordinates().y() > 0.39 * ( entity.getMesh().getDimensions().y() - 1 ) )
+ &&( entity.getCoordinates().z() < 0.59 * ( entity.getMesh().getDimensions().z() - 1 ) ) && ( entity.getCoordinates().z() > 0.39 * ( entity.getMesh().getDimensions().z() - 1 ) ) )
+ return ( (* this->pressure)[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ]
+ / ( this->gamma - 1 )
+ )
+ + 0.5
+ * (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ * (
+ this->cavitySpeed * ( -1 )
+ *
+ this->cavitySpeed * ( -1 )
+ +
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ +
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 2 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 2 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ );
+ else if( entity.getCoordinates().y() >0.8 * ( entity.getMesh().getDimensions().y() - 1 ) )
+ return u[ neighborEntities.template getEntityIndex< -1, 0, 0 >() ];
+ else
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ if( ( entity.getCoordinates().x() < 0.59 * ( entity.getMesh().getDimensions().x() - 1 ) ) && ( entity.getCoordinates().x() > 0.39 * ( entity.getMesh().getDimensions().x() - 1 ) )
+ &&( entity.getCoordinates().z() < 0.59 * ( entity.getMesh().getDimensions().z() - 1 ) ) && ( entity.getCoordinates().z() > 0.39 * ( entity.getMesh().getDimensions().z() - 1 ) ) )
+ return ( (* this->pressure)[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ]
+ / ( this->gamma - 1 )
+ )
+ + 0.5
+ * (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ * (
+ this->cavitySpeed
+ *
+ this->cavitySpeed
+ +
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ +
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 2 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 2 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ );
+ else
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ if( ( entity.getCoordinates().x() < 0.59 * ( entity.getMesh().getDimensions().x() - 1 ) ) && ( entity.getCoordinates().x() > 0.39 * ( entity.getMesh().getDimensions().x() - 1 ) )
+ &&( entity.getCoordinates().z() < 0.59 * ( entity.getMesh().getDimensions().z() - 1 ) ) && ( entity.getCoordinates().z() > 0.39 * ( entity.getMesh().getDimensions().z() - 1 ) ) )
+ return ( (* this->pressure)[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ]
+ / ( this->gamma - 1 )
+ )
+ + 0.5
+ * (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ * (
+ this->cavitySpeed * ( -1 )
+ *
+ this->cavitySpeed * ( -1 )
+ +
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ +
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 2 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 2 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ );
+ else
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ if( entity.getCoordinates().z() == 0 )
+{
+ if( ( entity.getCoordinates().x() < 0.6 * ( entity.getMesh().getDimensions().y() - 1 ) ) && ( entity.getCoordinates().y() > 0.4 * ( entity.getMesh().getDimensions().y() - 1 ) )
+ &&( entity.getCoordinates().y() < 0.6 * ( entity.getMesh().getDimensions().z() - 1 ) ) && ( entity.getCoordinates().z() > 0.4 * ( entity.getMesh().getDimensions().z() - 1 ) ) )
+ return ( (* this->pressure)[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ]
+ / ( this->gamma - 1 )
+ )
+ + 0.5
+ * (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ *
+ (
+ (
+ (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ +
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ +
+ (
+ this->cavitySpeed
+ *
+ this->cavitySpeed
+ )
+ );
+
+ else
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ // 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, 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, 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 EnergyBoundaryConditionsBoiler< Mesh, Function, Real, Index >& bc )
+{
+ str << "Neumann boundary ConditionsBoiler: function = " << bc.getFunction();
+ return str;
+}
+
+} // namespace Operators
+} // namespace TNL
+
diff --git a/examples/flow/EnergyBoundaryConditionCavity.h b/examples/flow/EnergyBoundaryConditionCavity.h
new file mode 100644
index 0000000000..9d1a55176e
--- /dev/null
+++ b/examples/flow/EnergyBoundaryConditionCavity.h
@@ -0,0 +1,673 @@
+/***************************************************************************
+ 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 EnergyBoundaryConditionsCavity
+{
+
+};
+
+/****
+ * Base
+ */
+template< typename Function>
+class EnergyBoundaryConditionsCavityBase
+{
+ 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 EnergyBoundaryConditionsCavity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public EnergyBoundaryConditionsCavityBase< 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 EnergyBoundaryConditionsCavity< MeshType, Function, Real, Index > ThisType;
+ typedef EnergyBoundaryConditionsCavityBase< 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 EnergyBoundaryConditionsCavity< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public EnergyBoundaryConditionsCavityBase< 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 EnergyBoundaryConditionsCavity< MeshType, Function, Real, Index > ThisType;
+ typedef EnergyBoundaryConditionsCavityBase< 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 >() ];
+ /*( (* this->pressure)[ neighborEntities.template getEntityIndex< 0, 0 >() ]
+ / ( this->gamma - 1 )
+ )
+ + 0.5
+ * (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ * (
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< 1, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 1, 0 >()]
+ + 0
+ )
+ *
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< 1, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 1, 0 >()]
+ + 0
+ )
+ );*/
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ /*( (* this->pressure)[ neighborEntities.template getEntityIndex< 0, 0 >() ]
+ / ( this->gamma - 1 )
+ )
+ + 0.5
+ * (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ * (
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< -1, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< -1, 0 >()]
+ + 0
+ )
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< -1, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< -1, 0 >()]
+ + 0
+ )
+ );*/
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ /*( (* this->pressure)[ neighborEntities.template getEntityIndex< 0, 0 >() ]
+ / ( this->gamma - 1 )
+ )
+ + 0.5
+ * (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ * (
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0, 1 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 1 >()]
+ + 0
+ )
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0, 1 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 1 >()]
+ + 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/*
+ * (
+ entity.getMesh().getDimensions().x() / 2 - std::abs( (entity.getCoordinates().x() - entity.getMesh().getDimensions().x() / 2 ) )
+ )
+ / ( entity.getMesh().getDimensions().x() / 2 )*/
+ *
+ this->cavitySpeed/*
+ * (
+ entity.getMesh().getDimensions().x() / 2 - std::abs( (entity.getCoordinates().x() - entity.getMesh().getDimensions().x() / 2 ) )
+ )
+ / ( entity.getMesh().getDimensions().x() / 2 )*/
+ +
+ ( (* (* 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 EnergyBoundaryConditionsCavity< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public EnergyBoundaryConditionsCavityBase< 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 EnergyBoundaryConditionsCavity< MeshType, Function, Real, Index > ThisType;
+ typedef EnergyBoundaryConditionsCavityBase< 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, 0 >() ];
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ if( entity.getCoordinates().z() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ // The following line is commented to avoid compiler warning
+ //if( entity.getCoordinates().z() == entity.getMesh().getDimensions().z() - 1 )
+ {
+ return ( (* this->pressure)[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ]
+ / ( this->gamma - 1 )
+ )
+ + 0.5
+ * (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ * (
+ ( this->cavitySpeed
+ * (
+ entity.getMesh().getDimensions().x() / 2 - std::abs( (entity.getCoordinates().x() - entity.getMesh().getDimensions().x() / 2 ) )
+ )
+ / ( entity.getMesh().getDimensions().x() / 2 )
+ )
+ *
+ ( this->cavitySpeed
+ * (
+ entity.getMesh().getDimensions().x() / 2 - std::abs( (entity.getCoordinates().x() - entity.getMesh().getDimensions().x() / 2 ) )
+ )
+ / ( entity.getMesh().getDimensions().x() / 2 )
+ )
+ +
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ +
+ ( (* (* this->compressibleConservativeVariables->getMomentum())[ 2 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ + 0
+ )
+ * ( (* (* this->compressibleConservativeVariables->getMomentum())[ 2 ])[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 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, 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 EnergyBoundaryConditionsCavity< Mesh, Function, Real, Index >& bc )
+{
+ str << "Neumann boundary ConditionsCavity: function = " << bc.getFunction();
+ return str;
+}
+
+} // namespace Operators
+} // namespace TNL
+
diff --git a/examples/flow/LaxFridrichsEnergy.h b/examples/flow/LaxFridrichsEnergy.h
index 7714edf013..fc4f7a4a83 100644
--- a/examples/flow/LaxFridrichsEnergy.h
+++ b/examples/flow/LaxFridrichsEnergy.h
@@ -233,36 +233,36 @@ class LaxFridrichsEnergy< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real,
+ ( ( ( e[ north ] + pressure_north ) * velocity_y_north )
-( ( e[ south ] + pressure_south ) * velocity_y_south ) ) * hyInverse )
// 2D uT_11_x
- - ( 4.0 / 3.0 * ( velocity_x_east * velocity_x_center - velocity_x_center * velocity_x_west
+ + ( 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 * 4
+ ) * hxSquareInverse
- 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 * 4
+ ) * hxInverse * hyInverse / 4
) * this->dynamicalViscosity
// vT_21_x
- - ( ( velocity_y_northEast * velocity_y_east - velocity_y_southEast * velocity_y_east
+ + ( ( 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 * 4
+ ) * 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 * 4
+ ) * hxSquareInverse
) * this->dynamicalViscosity
// uT_12_y
- - ( ( velocity_x_northEast * velocity_x_north - velocity_x_southEast * velocity_x_south
+ + ( ( 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 * 4
+ ) * 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 * 4
+ ) * hySquareInverse
) * this->dynamicalViscosity
// 2D vT_22_y
- - ( 4.0 / 3.0 * ( velocity_y_north * velocity_y_center - velocity_y_center * velocity_y_south
+ + ( 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 * 4
+ ) * hySquareInverse
- 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 * 4
+ ) * hxInverse * hyInverse / 4
) * this->dynamicalViscosity;
}
diff --git a/examples/flow/LaxFridrichsMomentumX.h b/examples/flow/LaxFridrichsMomentumX.h
index d1c74e605f..c2d1886a19 100644
--- a/examples/flow/LaxFridrichsMomentumX.h
+++ b/examples/flow/LaxFridrichsMomentumX.h
@@ -163,6 +163,8 @@ class LaxFridrichsMomentumX< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Rea
const RealType& pressure_east = this->pressure.template getData< DeviceType >()[ east ];
const RealType& velocity_x_east = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ];
const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
+ const RealType& velocity_x_north = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ north ];
+ const RealType& velocity_x_south = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ south ];
const RealType& velocity_x_center = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ center ];
const RealType& velocity_x_southEast = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ southEast ];
const RealType& velocity_x_southWest = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ southWest ];
@@ -182,16 +184,16 @@ class LaxFridrichsMomentumX< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Rea
+ ( ( rho_u[ north ] * velocity_y_north )
- ( 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 * 4
+ + ( 4.0 / 3.0 * ( velocity_x_east - 2 * velocity_x_center + velocity_x_west
+ ) * hxSquareInverse
- 2.0 / 3.0 * ( velocity_y_northEast - velocity_y_southEast - velocity_y_northWest + velocity_y_southWest
- ) * hxInverse * hyInverse * 4
+ ) * hxInverse * hyInverse / 4
) * this->dynamicalViscosity
-// T_21_x
- - ( ( velocity_y_northEast - velocity_y_southEast - velocity_y_northWest + velocity_y_southWest
- ) * hxInverse * hyInverse * 4
- + ( velocity_x_east - 2 * velocity_x_center + velocity_x_west
- ) * hxInverse * hyInverse * 4
+// T_21_y
+ + ( ( velocity_y_northEast - velocity_y_southEast - velocity_y_northWest + velocity_y_southWest
+ ) * hxInverse * hyInverse / 4
+ + ( velocity_x_north - 2 * velocity_x_center + velocity_x_south
+ ) * hxInverse * hyInverse
) * this->dynamicalViscosity;
}
diff --git a/examples/flow/LaxFridrichsMomentumY.h b/examples/flow/LaxFridrichsMomentumY.h
index c0f2ead15f..17e7133724 100644
--- a/examples/flow/LaxFridrichsMomentumY.h
+++ b/examples/flow/LaxFridrichsMomentumY.h
@@ -153,6 +153,8 @@ class LaxFridrichsMomentumY< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Rea
const RealType& velocity_x_northWest = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ northWest ];
const RealType& velocity_y_north = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ north ];
const RealType& velocity_y_south = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ];
+ const RealType& velocity_y_east = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ east ];
+ const RealType& velocity_y_west = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ west ];
const RealType& velocity_y_center = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ center ];
const RealType& velocity_y_southEast = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ southEast ];
const RealType& velocity_y_southWest = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ southWest ];
@@ -166,16 +168,16 @@ class LaxFridrichsMomentumY< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Rea
+ ( ( rho_v[ north ] * velocity_y_north + pressure_north )
- ( 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 * 4
+ + ( 4.0 / 3.0 * ( velocity_y_north - 2 * velocity_y_center + velocity_y_south
+ ) * hySquareInverse
- 2.0 / 3.0 * ( velocity_x_northEast - velocity_x_southEast - velocity_x_northWest + velocity_x_southWest
- ) * hxInverse * hyInverse * 4
+ ) * hxInverse * hyInverse / 4
) * this->dynamicalViscosity
-// T_12_y
- - ( ( velocity_x_northEast - velocity_x_southEast - velocity_x_northWest + velocity_x_southWest
- ) * hxInverse * hyInverse * 4
- + ( velocity_y_north - 2 * velocity_y_center + velocity_y_south
- ) * hySquareInverse * 4
+// T_12_x
+ + ( ( velocity_x_northEast - velocity_x_southEast - velocity_x_northWest + velocity_x_southWest
+ ) * hxInverse * hyInverse / 4
+ + ( velocity_y_west - 2 * velocity_y_center + velocity_y_east
+ ) * hxSquareInverse
) * this->dynamicalViscosity;
}
diff --git a/examples/flow/MomentumXBoundaryConditionBoiler.h b/examples/flow/MomentumXBoundaryConditionBoiler.h
new file mode 100644
index 0000000000..823ec475a5
--- /dev/null
+++ b/examples/flow/MomentumXBoundaryConditionBoiler.h
@@ -0,0 +1,594 @@
+/***************************************************************************
+ 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 MomentumXBoundaryConditionsBoiler
+{
+
+};
+
+/****
+ * Base
+ */
+template< typename Function >
+class MomentumXBoundaryConditionsBoilerBase
+{
+ 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 MomentumXBoundaryConditionsBoiler< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumXBoundaryConditionsBoilerBase< 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 MomentumXBoundaryConditionsBoiler< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumXBoundaryConditionsBoilerBase< 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 MomentumXBoundaryConditionsBoiler< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumXBoundaryConditionsBoilerBase< 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 MomentumXBoundaryConditionsBoiler< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumXBoundaryConditionsBoilerBase< 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 )
+ {
+ if( ( entity.getCoordinates().y() < 0.20 * ( entity.getMesh().getDimensions().y() - 1 ) ) && ( entity.getCoordinates().y() > 0.19 * ( entity.getMesh().getDimensions().y() - 1 ) ) )
+ return (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ * (
+ this->cavitySpeed
+ );
+ else
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ if( ( entity.getCoordinates().y() < 0.20 * ( entity.getMesh().getDimensions().y() - 1 ) ) && ( entity.getCoordinates().y() > 0.19 * ( entity.getMesh().getDimensions().y() - 1 ) ) )
+ return (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ * (
+ this->cavitySpeed * ( -1 )
+ );
+ else if( entity.getCoordinates().y() > 0.8 * ( entity.getMesh().getDimensions().y() - 1 ) && false )
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ else
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ /*(* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ * (
+ (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0, 1 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 1 >()]
+ );*/
+ }
+ // The following line is commented to avoid compiler warning
+ //if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, -1 >() ];
+ /*return (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ * (
+ (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0, -1 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, -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, 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 MomentumXBoundaryConditionsBoiler< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumXBoundaryConditionsBoilerBase< 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 MomentumXBoundaryConditionsBoiler< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumXBoundaryConditionsBoilerBase< 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 )
+ {
+ if( ( entity.getCoordinates().y() < 0.59 * ( entity.getMesh().getDimensions().y() - 1 ) ) && ( entity.getCoordinates().y() > 0.39 * ( entity.getMesh().getDimensions().y() - 1 ) )
+ &&( entity.getCoordinates().z() < 0.59 * ( entity.getMesh().getDimensions().z() - 1 ) ) && ( entity.getCoordinates().z() > 0.39 * ( entity.getMesh().getDimensions().z() - 1 ) ) )
+ return (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ * (
+ this->cavitySpeed
+ );
+ else
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ if( ( entity.getCoordinates().y() < 0.59 * ( entity.getMesh().getDimensions().y() - 1 ) ) && ( entity.getCoordinates().y() > 0.39 * ( entity.getMesh().getDimensions().y() - 1 ) )
+ &&( entity.getCoordinates().z() < 0.59 * ( entity.getMesh().getDimensions().z() - 1 ) ) && ( entity.getCoordinates().z() > 0.39 * ( entity.getMesh().getDimensions().z() - 1 ) ) )
+ return (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ * (
+ this->cavitySpeed * ( -1 )
+ );
+ else if( entity.getCoordinates().y() >0.8 * ( entity.getMesh().getDimensions().y() - 1 ) )
+ return u[ neighborEntities.template getEntityIndex< -1, 0, 0 >() ];
+ else
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ if( entity.getCoordinates().z() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ // 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, 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, 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 MomentumXBoundaryConditionsBoiler< Mesh, Function, Real, Index >& bc )
+{
+ str << "Neumann boundary ConditionsBoiler: function = " << bc.getFunction();
+ return str;
+}
+
+} // namespace Operators
+} // namespace TNL
+
diff --git a/examples/flow/1DMomentumXBoundaryCondition.h b/examples/flow/MomentumXBoundaryConditionCavity.h
similarity index 90%
rename from examples/flow/1DMomentumXBoundaryCondition.h
rename to examples/flow/MomentumXBoundaryConditionCavity.h
index 101050bc09..5d3be8ba0e 100644
--- a/examples/flow/1DMomentumXBoundaryCondition.h
+++ b/examples/flow/MomentumXBoundaryConditionCavity.h
@@ -20,7 +20,7 @@ template< typename Mesh,
typename Function,
typename Real = typename Mesh::RealType,
typename Index = typename Mesh::GlobalIndexType >
-class MomentumXBoundaryConditions
+class MomentumXBoundaryConditionsCavity
{
};
@@ -29,7 +29,7 @@ class MomentumXBoundaryConditions
* Base
*/
template< typename Function >
-class MomentumXBoundaryConditionsBase
+class MomentumXBoundaryConditionsCavityBase
{
public:
@@ -91,8 +91,8 @@ template< typename MeshReal,
typename Function,
typename Real,
typename Index >
-class MomentumXBoundaryConditions< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >
- : public MomentumXBoundaryConditionsBase< Function >,
+class MomentumXBoundaryConditionsCavity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumXBoundaryConditionsCavityBase< Function >,
public Operator< Meshes::Grid< 1, MeshReal, Device, MeshIndex >,
Functions::MeshBoundaryDomain,
1, 1,
@@ -111,8 +111,8 @@ class MomentumXBoundaryConditions< Meshes::Grid< 1, MeshReal, Device, MeshIndex
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 MomentumXBoundaryConditionsCavity< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumXBoundaryConditionsCavityBase< Function > BaseType;
typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
@@ -221,8 +221,8 @@ template< typename MeshReal,
typename Function,
typename Real,
typename Index >
-class MomentumXBoundaryConditions< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >
- : public MomentumXBoundaryConditionsBase< Function >,
+class MomentumXBoundaryConditionsCavity< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumXBoundaryConditionsCavityBase< Function >,
public Operator< Meshes::Grid< 2, MeshReal, Device, MeshIndex >,
Functions::MeshBoundaryDomain,
2, 2,
@@ -242,8 +242,8 @@ class MomentumXBoundaryConditions< Meshes::Grid< 2, MeshReal, Device, MeshIndex
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 MomentumXBoundaryConditionsCavity< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumXBoundaryConditionsCavityBase< Function > BaseType;
typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
@@ -265,18 +265,28 @@ class MomentumXBoundaryConditions< Meshes::Grid< 2, MeshReal, Device, MeshIndex
}
if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
{
- return u[ neighborEntities.template getEntityIndex< -0, 0 >() ];
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
}
if( entity.getCoordinates().y() == 0 )
{
return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ /*(* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ * (
+ (* (* this->compressibleConservativeVariables->getMomentum())[ 0 ])[neighborEntities.template getEntityIndex< 0, 1 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 1 >()]
+ );*/
}
// 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
+ ( this->cavitySpeed/*
+ * (
+ entity.getMesh().getDimensions().x() / 2 - std::abs( (entity.getCoordinates().x() - entity.getMesh().getDimensions().x() / 2 ) )
+ )
+ / ( entity.getMesh().getDimensions().x() / 2 )*/
+ )
);
}
}
@@ -377,8 +387,8 @@ template< typename MeshReal,
typename Function,
typename Real,
typename Index >
-class MomentumXBoundaryConditions< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >
- : public MomentumXBoundaryConditionsBase< Function >,
+class MomentumXBoundaryConditionsCavity< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumXBoundaryConditionsCavityBase< Function >,
public Operator< Meshes::Grid< 3, MeshReal, Device, MeshIndex >,
Functions::MeshBoundaryDomain,
3, 3,
@@ -397,8 +407,8 @@ class MomentumXBoundaryConditions< Meshes::Grid< 3, MeshReal, Device, MeshIndex
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 MomentumXBoundaryConditionsCavity< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumXBoundaryConditionsCavityBase< Function > BaseType;
typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
@@ -415,34 +425,31 @@ class MomentumXBoundaryConditions< Meshes::Grid< 3, MeshReal, Device, MeshIndex
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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
// 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 );
+ return (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ * (
+ this->cavitySpeed
+ );
}
}
@@ -552,9 +559,9 @@ template< typename Mesh,
typename Function,
typename Real,
typename Index >
-std::ostream& operator << ( std::ostream& str, const MomentumXBoundaryConditions< Mesh, Function, Real, Index >& bc )
+std::ostream& operator << ( std::ostream& str, const MomentumXBoundaryConditionsCavity< Mesh, Function, Real, Index >& bc )
{
- str << "Neumann boundary conditions: function = " << bc.getFunction();
+ str << "Neumann boundary ConditionsCavity: function = " << bc.getFunction();
return str;
}
diff --git a/examples/flow/MomentumYBoundaryConditionBoiler.h b/examples/flow/MomentumYBoundaryConditionBoiler.h
new file mode 100644
index 0000000000..76f3ff0573
--- /dev/null
+++ b/examples/flow/MomentumYBoundaryConditionBoiler.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>
+
+namespace TNL {
+namespace Operators {
+
+template< typename Mesh,
+ typename Function,
+ typename Real = typename Mesh::RealType,
+ typename Index = typename Mesh::GlobalIndexType >
+class MomentumYBoundaryConditionsBoiler
+{
+
+};
+
+/****
+ * Base
+ */
+template< typename Function >
+class MomentumYBoundaryConditionsBoilerBase
+{
+ 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 MomentumYBoundaryConditionsBoiler< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumYBoundaryConditionsBoilerBase< 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 MomentumYBoundaryConditionsBoiler< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumYBoundaryConditionsBoilerBase< 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 MomentumYBoundaryConditionsBoiler< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumYBoundaryConditionsBoilerBase< 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 MomentumYBoundaryConditionsBoiler< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumYBoundaryConditionsBoilerBase< 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 >() ];
+ /*(* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ * (
+ (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< 1, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 1, 0 >()]
+ );*/
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ if( entity.getCoordinates().y() > 0.8 * ( entity.getMesh().getDimensions().y() - 1 ) && false)
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ else
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ /*(* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ * (
+ (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< -1, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< -1, 0 >()]
+ );*/
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ if( ( entity.getCoordinates().x() < 0.6 * ( entity.getMesh().getDimensions().x() - 1 ) ) && ( entity.getCoordinates().x() > 0.4 * ( entity.getMesh().getDimensions().x() - 1 ) ) )
+ return (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ * (
+ this->cavitySpeed
+ );
+ else 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, -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 MomentumYBoundaryConditionsBoiler< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumYBoundaryConditionsBoilerBase< 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 MomentumYBoundaryConditionsBoiler< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumYBoundaryConditionsBoilerBase< 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, 0 >() ];
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ if( entity.getCoordinates().y() >0.8 * ( entity.getMesh().getDimensions().y() - 1 ) )
+ return u[ neighborEntities.template getEntityIndex< -1, 0, 0 >() ];
+ else
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ if( ( entity.getCoordinates().x() < 0.59 * ( entity.getMesh().getDimensions().x() - 1 ) ) && ( entity.getCoordinates().x() > 0.39 * ( entity.getMesh().getDimensions().x() - 1 ) )
+ &&( entity.getCoordinates().z() < 0.59 * ( entity.getMesh().getDimensions().z() - 1 ) ) && ( entity.getCoordinates().z() > 0.39 * ( entity.getMesh().getDimensions().z() - 1 ) ) )
+ return (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ * (
+ this->cavitySpeed
+ );
+ else return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ if( ( entity.getCoordinates().x() < 0.59 * ( entity.getMesh().getDimensions().x() - 1 ) ) && ( entity.getCoordinates().x() > 0.39 * ( entity.getMesh().getDimensions().x() - 1 ) )
+ &&( entity.getCoordinates().z() < 0.59 * ( entity.getMesh().getDimensions().z() - 1 ) ) && ( entity.getCoordinates().z() > 0.39 * ( entity.getMesh().getDimensions().z() - 1 ) ) )
+ return (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ * (
+ this->cavitySpeed * ( -1 )
+ );
+ else return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ if( entity.getCoordinates().z() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ // 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, 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, 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 MomentumYBoundaryConditionsBoiler< Mesh, Function, Real, Index >& bc )
+{
+ str << "Neumann boundary ConditionsBoiler: function = " << bc.getFunction();
+ return str;
+}
+
+} // namespace Operators
+} // namespace TNL
+
diff --git a/examples/flow/1DEnergyBoundaryCondition.h b/examples/flow/MomentumYBoundaryConditionCavity.h
similarity index 83%
rename from examples/flow/1DEnergyBoundaryCondition.h
rename to examples/flow/MomentumYBoundaryConditionCavity.h
index bce86634a9..afce8239b7 100644
--- a/examples/flow/1DEnergyBoundaryCondition.h
+++ b/examples/flow/MomentumYBoundaryConditionCavity.h
@@ -12,7 +12,6 @@
#pragma once
#include <TNL/Functions/FunctionAdapter.h>
-#include "CompressibleConservativeVariables.h"
namespace TNL {
namespace Operators {
@@ -21,7 +20,7 @@ template< typename Mesh,
typename Function,
typename Real = typename Mesh::RealType,
typename Index = typename Mesh::GlobalIndexType >
-class EnergyBoundaryConditions
+class MomentumYBoundaryConditionsCavity
{
};
@@ -29,8 +28,8 @@ class EnergyBoundaryConditions
/****
* Base
*/
-template< typename Function>
-class EnergyBoundaryConditionsBase
+template< typename Function >
+class MomentumYBoundaryConditionsCavityBase
{
public:
@@ -81,7 +80,6 @@ class EnergyBoundaryConditionsBase
FunctionType function;
-
};
/****
@@ -93,8 +91,8 @@ template< typename MeshReal,
typename Function,
typename Real,
typename Index >
-class EnergyBoundaryConditions< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >
- : public EnergyBoundaryConditionsBase< Function >,
+class MomentumYBoundaryConditionsCavity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumYBoundaryConditionsCavityBase< Function >,
public Operator< Meshes::Grid< 1, MeshReal, Device, MeshIndex >,
Functions::MeshBoundaryDomain,
1, 1,
@@ -113,8 +111,8 @@ class EnergyBoundaryConditions< Meshes::Grid< 1, MeshReal, Device, MeshIndex >,
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 MomentumYBoundaryConditionsCavity< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumYBoundaryConditionsCavityBase< Function > BaseType;
typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
@@ -130,19 +128,11 @@ class EnergyBoundaryConditions< Meshes::Grid< 1, MeshReal, Device, MeshIndex >,
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
- );
+ 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 >() ];
@@ -231,8 +221,8 @@ template< typename MeshReal,
typename Function,
typename Real,
typename Index >
-class EnergyBoundaryConditions< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >
- : public EnergyBoundaryConditionsBase< Function >,
+class MomentumYBoundaryConditionsCavity< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumYBoundaryConditionsCavityBase< Function >,
public Operator< Meshes::Grid< 2, MeshReal, Device, MeshIndex >,
Functions::MeshBoundaryDomain,
2, 2,
@@ -252,8 +242,8 @@ class EnergyBoundaryConditions< Meshes::Grid< 2, MeshReal, Device, MeshIndex >,
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 MomentumYBoundaryConditionsCavity< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumYBoundaryConditionsCavityBase< Function > BaseType;
typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
@@ -272,10 +262,20 @@ class EnergyBoundaryConditions< Meshes::Grid< 2, MeshReal, Device, MeshIndex >,
if( entity.getCoordinates().x() == 0 )
{
return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ /*(* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ * (
+ (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< 1, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 1, 0 >()]
+ );*/
}
if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
{
return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
+ /*(* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0 >()]
+ * (
+ (* (* this->compressibleConservativeVariables->getMomentum())[ 1 ])[neighborEntities.template getEntityIndex< -1, 0 >()]
+ / (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< -1, 0 >()]
+ );*/
}
if( entity.getCoordinates().y() == 0 )
{
@@ -284,25 +284,7 @@ class EnergyBoundaryConditions< Meshes::Grid< 2, MeshReal, Device, MeshIndex >,
// 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
- )
- );
+ return u[ neighborEntities.template getEntityIndex< 0, 0 >() ];
}
}
@@ -402,8 +384,8 @@ template< typename MeshReal,
typename Function,
typename Real,
typename Index >
-class EnergyBoundaryConditions< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >
- : public EnergyBoundaryConditionsBase< Function >,
+class MomentumYBoundaryConditionsCavity< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumYBoundaryConditionsCavityBase< Function >,
public Operator< Meshes::Grid< 3, MeshReal, Device, MeshIndex >,
Functions::MeshBoundaryDomain,
3, 3,
@@ -422,8 +404,8 @@ class EnergyBoundaryConditions< Meshes::Grid< 3, MeshReal, Device, MeshIndex >,
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 MomentumYBoundaryConditionsCavity< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumYBoundaryConditionsCavityBase< Function > BaseType;
typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
@@ -440,34 +422,28 @@ class EnergyBoundaryConditions< Meshes::Grid< 3, MeshReal, Device, MeshIndex >,
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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
// 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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
}
@@ -577,9 +553,9 @@ template< typename Mesh,
typename Function,
typename Real,
typename Index >
-std::ostream& operator << ( std::ostream& str, const EnergyBoundaryConditions< Mesh, Function, Real, Index >& bc )
+std::ostream& operator << ( std::ostream& str, const MomentumYBoundaryConditionsCavity< Mesh, Function, Real, Index >& bc )
{
- str << "Neumann boundary conditions: function = " << bc.getFunction();
+ str << "Neumann boundary ConditionsCavity: function = " << bc.getFunction();
return str;
}
diff --git a/examples/flow/MomentumZBoundaryConditionBoiler.h b/examples/flow/MomentumZBoundaryConditionBoiler.h
new file mode 100644
index 0000000000..188aaa9851
--- /dev/null
+++ b/examples/flow/MomentumZBoundaryConditionBoiler.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 MomentumZBoundaryConditionsBoiler
+{
+
+};
+
+/****
+ * Base
+ */
+template< typename Function >
+class MomentumZBoundaryConditionsBoilerBase
+{
+ 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 MomentumZBoundaryConditionsBoiler< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumZBoundaryConditionsBoilerBase< 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 MomentumZBoundaryConditionsBoiler< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumZBoundaryConditionsBoilerBase< 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 MomentumZBoundaryConditionsBoiler< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumZBoundaryConditionsBoilerBase< 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 MomentumZBoundaryConditionsBoiler< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumZBoundaryConditionsBoilerBase< 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 MomentumZBoundaryConditionsBoiler< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumZBoundaryConditionsBoilerBase< 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 MomentumZBoundaryConditionsBoiler< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumZBoundaryConditionsBoilerBase< 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, 0 >() ];
+ }
+ if( entity.getCoordinates().x() == entity.getMesh().getDimensions().x() - 1 )
+ {
+ if( entity.getCoordinates().y() >0.8 * ( entity.getMesh().getDimensions().y() - 1 ) )
+ return u[ neighborEntities.template getEntityIndex< -1, 0, 0 >() ];
+ else
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ if( entity.getCoordinates().y() == 0 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ if( entity.getCoordinates().y() == entity.getMesh().getDimensions().y() - 1 )
+ {
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ if( entity.getCoordinates().z() == 0 )
+ {
+ if( ( entity.getCoordinates().x() < 0.6 * ( entity.getMesh().getDimensions().y() - 1 ) ) && ( entity.getCoordinates().y() > 0.4 * ( entity.getMesh().getDimensions().y() - 1 ) )
+ &&( entity.getCoordinates().y() < 0.6 * ( entity.getMesh().getDimensions().z() - 1 ) ) && ( entity.getCoordinates().z() > 0.4 * ( entity.getMesh().getDimensions().z() - 1 ) ) )
+ return (* this->compressibleConservativeVariables->getDensity())[neighborEntities.template getEntityIndex< 0, 0, 0 >()]
+ * (
+ this->cavitySpeed
+ );
+ else return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
+ }
+ // 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, 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, 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 MomentumZBoundaryConditionsBoiler< Mesh, Function, Real, Index >& bc )
+{
+ str << "Neumann boundary ConditionsBoiler: function = " << bc.getFunction();
+ return str;
+}
+
+} // namespace Operators
+} // namespace TNL
+
diff --git a/examples/flow/1DMomentumYBoundaryCondition.h b/examples/flow/MomentumZBoundaryConditionCavity.h
similarity index 90%
rename from examples/flow/1DMomentumYBoundaryCondition.h
rename to examples/flow/MomentumZBoundaryConditionCavity.h
index 1a8637bbac..1942cd5893 100644
--- a/examples/flow/1DMomentumYBoundaryCondition.h
+++ b/examples/flow/MomentumZBoundaryConditionCavity.h
@@ -20,7 +20,7 @@ template< typename Mesh,
typename Function,
typename Real = typename Mesh::RealType,
typename Index = typename Mesh::GlobalIndexType >
-class MomentumYBoundaryConditions
+class MomentumZBoundaryConditionsCavity
{
};
@@ -29,7 +29,7 @@ class MomentumYBoundaryConditions
* Base
*/
template< typename Function >
-class MomentumYBoundaryConditionsBase
+class MomentumZBoundaryConditionsCavityBase
{
public:
@@ -91,8 +91,8 @@ template< typename MeshReal,
typename Function,
typename Real,
typename Index >
-class MomentumYBoundaryConditions< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >
- : public MomentumYBoundaryConditionsBase< Function >,
+class MomentumZBoundaryConditionsCavity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumZBoundaryConditionsCavityBase< Function >,
public Operator< Meshes::Grid< 1, MeshReal, Device, MeshIndex >,
Functions::MeshBoundaryDomain,
1, 1,
@@ -111,8 +111,8 @@ class MomentumYBoundaryConditions< Meshes::Grid< 1, MeshReal, Device, MeshIndex
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 MomentumZBoundaryConditionsCavity< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumZBoundaryConditionsCavityBase< Function > BaseType;
typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
@@ -221,8 +221,8 @@ template< typename MeshReal,
typename Function,
typename Real,
typename Index >
-class MomentumYBoundaryConditions< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >
- : public MomentumYBoundaryConditionsBase< Function >,
+class MomentumZBoundaryConditionsCavity< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumZBoundaryConditionsCavityBase< Function >,
public Operator< Meshes::Grid< 2, MeshReal, Device, MeshIndex >,
Functions::MeshBoundaryDomain,
2, 2,
@@ -242,8 +242,8 @@ class MomentumYBoundaryConditions< Meshes::Grid< 2, MeshReal, Device, MeshIndex
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 MomentumZBoundaryConditionsCavity< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumZBoundaryConditionsCavityBase< Function > BaseType;
typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
@@ -374,8 +374,8 @@ template< typename MeshReal,
typename Function,
typename Real,
typename Index >
-class MomentumYBoundaryConditions< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >
- : public MomentumYBoundaryConditionsBase< Function >,
+class MomentumZBoundaryConditionsCavity< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Function, Real, Index >
+ : public MomentumZBoundaryConditionsCavityBase< Function >,
public Operator< Meshes::Grid< 3, MeshReal, Device, MeshIndex >,
Functions::MeshBoundaryDomain,
3, 3,
@@ -394,8 +394,8 @@ class MomentumYBoundaryConditions< Meshes::Grid< 3, MeshReal, Device, MeshIndex
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 MomentumZBoundaryConditionsCavity< MeshType, Function, Real, Index > ThisType;
+ typedef MomentumZBoundaryConditionsCavityBase< Function > BaseType;
typedef CompressibleConservativeVariables< MeshType > CompressibleConservativeVariablesType;
typedef SharedPointer< CompressibleConservativeVariablesType > CompressibleConservativeVariablesPointer;
typedef SharedPointer< MeshFunctionType, DeviceType > MeshFunctionPointer;
@@ -412,34 +412,28 @@ class MomentumYBoundaryConditions< Meshes::Grid< 3, MeshReal, Device, MeshIndex
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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
// 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 );
+ return u[ neighborEntities.template getEntityIndex< 0, 0, 0 >() ];
}
}
@@ -549,9 +543,9 @@ template< typename Mesh,
typename Function,
typename Real,
typename Index >
-std::ostream& operator << ( std::ostream& str, const MomentumYBoundaryConditions< Mesh, Function, Real, Index >& bc )
+std::ostream& operator << ( std::ostream& str, const MomentumZBoundaryConditionsCavity< Mesh, Function, Real, Index >& bc )
{
- str << "Neumann boundary conditions: function = " << bc.getFunction();
+ str << "Neumann boundary ConditionsCavity: function = " << bc.getFunction();
return str;
}
diff --git a/examples/flow/navierStokes.h b/examples/flow/navierStokes.h
index eb257e8cf0..b0ea90d8bf 100644
--- a/examples/flow/navierStokes.h
+++ b/examples/flow/navierStokes.h
@@ -10,7 +10,8 @@
#include "navierStokesBuildConfigTag.h"
#include "RiemannProblemInitialCondition.h"
-#include "1DBoundaryConditions.h"
+#include "BoundaryConditionsCavity.h"
+#include "BoundaryConditionsBoiler.h"
using namespace TNL;
@@ -32,9 +33,9 @@ template< typename ConfigTag >class navierStokesConfig
static void configSetup( Config::ConfigDescription & config )
{
config.addDelimiter( "Inviscid flow settings:" );
- config.addEntry< String >( "boundary-conditions-type", "Choose the boundary conditions type.", "dirichlet");
- config.addEntryEnum< String >( "dirichlet" );
- config.addEntryEnum< String >( "neumann" );
+ config.addEntry< String >( "boundary-conditions-type", "Choose the boundary conditions type.", "cavity");
+ config.addEntryEnum< String >( "boiler" );
+ config.addEntryEnum< String >( "cavity" );
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 );
@@ -78,41 +79,21 @@ class navierStokesSetter
*/
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" ) )
- {
- typedef Functions::Analytic::Constant< Dimension, Real > Constant;
- if( boundaryConditionsType == "dirichlet" )
+ if( boundaryConditionsType == "cavity" )
{
- typedef Operators::DirichletBoundaryConditions< MeshType, Constant, MeshType::getMeshDimension(), Real, Index > BoundaryConditions;
+ typedef BoundaryConditionsCavity< MeshType, Constant, Real, Index > BoundaryConditions;
typedef navierStokesProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
SolverStarter solverStarter;
return solverStarter.template run< Problem >( parameters );
}
- typedef Operators::NeumannBoundaryConditions< MeshType, Constant, Real, Index > BoundaryConditions;
- typedef navierStokesProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
- SolverStarter solverStarter;
- return solverStarter.template run< Problem >( parameters );
- }
- typedef Functions::MeshFunction< MeshType > MeshFunction;
- if( boundaryConditionsType == "dirichlet" )
- {
- typedef Operators::DirichletBoundaryConditions< MeshType, MeshFunction, MeshType::getMeshDimension(), Real, Index > BoundaryConditions;
- typedef navierStokesProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
- SolverStarter solverStarter;
- return solverStarter.template run< Problem >( parameters );
- }
- if( boundaryConditionsType == "neumann" )
- {
- typedef Operators::NeumannBoundaryConditions< MeshType, MeshFunction, Real, Index > BoundaryConditions;
- typedef navierStokesProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
- SolverStarter solverStarter;
- return solverStarter.template run< Problem >( parameters );
- }*/
+ if( boundaryConditionsType == "boiler" )
+ {
+ typedef BoundaryConditionsBoiler< MeshType, Constant, Real, Index > BoundaryConditions;
+ typedef navierStokesProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;
+ SolverStarter solverStarter;
+ return solverStarter.template run< Problem >( parameters );
+ }
return true;}
diff --git a/examples/flow/navierStokesProblem_impl.h b/examples/flow/navierStokesProblem_impl.h
index 0b421ba16e..3c7e97dd20 100644
--- a/examples/flow/navierStokesProblem_impl.h
+++ b/examples/flow/navierStokesProblem_impl.h
@@ -27,12 +27,19 @@
#include "LaxFridrichsMomentumX.h"
#include "LaxFridrichsMomentumY.h"
#include "LaxFridrichsMomentumZ.h"
-
-#include "1DDensityBoundaryCondition.h"
-#include "1DMomentumXBoundaryCondition.h"
-#include "1DMomentumYBoundaryCondition.h"
-#include "1DEnergyBoundaryCondition.h"
-
+/*
+#include "DensityBoundaryConditionCavity.h"
+#include "MomentumXBoundaryConditionCavity.h"
+#include "MomentumYBoundaryConditionCavity.h"
+#include "MomentumZBoundaryConditionCavity.h"
+#include "EnergyBoundaryConditionCavity.h"
+
+#include "DensityBoundaryConditionBoiler.h"
+#include "MomentumXBoundaryConditionBoiler.h"
+#include "MomentumYBoundaryConditionBoiler.h"
+#include "MomentumZBoundaryConditionBoiler.h"
+#include "EnergyBoundaryConditionBoiler.h"
+*/
namespace TNL {
template< typename Mesh,
@@ -196,23 +203,24 @@ makeSnapshot( const RealType& time,
fileName.setExtension( "tnl" );
fileName.setIndex( step );
fileName.setFileNameBase( "density-" );
- if( ! this->conservativeVariables->getDensity()->save( fileName.getFileName() ) )
- return false;
+// if( ! this->conservativeVariables->getDensity()->save( fileName.getFileName() ) )
+// return false;
fileName.setFileNameBase( "velocity-" );
if( ! this->velocity->save( fileName.getFileName() ) )
return false;
- fileName.setFileNameBase( "pressure-" );
- if( ! this->pressure->save( fileName.getFileName() ) )
- return false;
+// fileName.setFileNameBase( "pressure-" );
+// if( ! this->pressure->save( fileName.getFileName() ) )
+// return false;
- fileName.setFileNameBase( "energy-" );
- if( ! this->conservativeVariables->getEnergy()->save( fileName.getFileName() ) )
- return false;
- fileName.setFileNameBase( "momentum-" );
- if( ! this->conservativeVariables->getMomentum()->save( fileName.getFileName() ) )
- return false;
+// 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;
}
@@ -267,6 +275,7 @@ getExplicitUpdate( const RealType& time,
typedef typename BoundaryCondition::DensityBoundaryConditionsType DensityBoundaryConditionsType;
typedef typename BoundaryCondition::MomentumXBoundaryConditionsType MomentumXBoundaryConditionsType;
typedef typename BoundaryCondition::MomentumYBoundaryConditionsType MomentumYBoundaryConditionsType;
+ typedef typename BoundaryCondition::MomentumZBoundaryConditionsType MomentumZBoundaryConditionsType;
typedef typename BoundaryCondition::EnergyBoundaryConditionsType EnergyBoundaryConditionsType;
/****
@@ -280,7 +289,7 @@ getExplicitUpdate( const RealType& time,
{
this->boundaryConditionPointer->setTimestep(0);
}
- this->boundaryConditionPointer->setCavitySpeed(this->cavitySpeed);
+ this->boundaryConditionPointer->setSpeed(this->cavitySpeed);
this->boundaryConditionPointer->setCompressibleConservativeVariables(this->conservativeVariables);
this->boundaryConditionPointer->setGamma(this->gamma);
this->boundaryConditionPointer->setPressure(this->pressure);
@@ -319,9 +328,9 @@ getExplicitUpdate( const RealType& time,
if( Dimensions > 2 )
{
- Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, MomentumZOperatorType, MomentumXBoundaryConditionsType, RightHandSide > explicitUpdaterMomentumZ;
+ Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, MomentumZOperatorType, MomentumZBoundaryConditionsType, RightHandSide > explicitUpdaterMomentumZ;
explicitUpdaterMomentumZ.setDifferentialOperator( this->inviscidOperatorsPointer->getMomentumZOperator() );
- explicitUpdaterMomentumZ.setBoundaryConditions( this->boundaryConditionPointer->getMomentumXBoundaryCondition() );
+ explicitUpdaterMomentumZ.setBoundaryConditions( this->boundaryConditionPointer->getMomentumZBoundaryCondition() );
explicitUpdaterMomentumZ.setRightHandSide( this->rightHandSidePointer );
explicitUpdaterMomentumZ.template update< typename Mesh::Cell >( time, tau, mesh,
( *this->conservativeVariables->getMomentum() )[ 2 ], // uRhoVelocityX,
@@ -339,10 +348,12 @@ getExplicitUpdate( const RealType& time,
this->conservativeVariables->getEnergy(), // uRhoVelocityX,
this->conservativeVariablesRHS->getEnergy() ); //, fuRhoVelocityX );
- /*this->conservativeVariablesRHS->getDensity()->write( "density", "gnuplot" );
+ /*
+ this->conservativeVariablesRHS->getDensity()->write( "density", "gnuplot" );
this->conservativeVariablesRHS->getEnergy()->write( "energy", "gnuplot" );
this->conservativeVariablesRHS->getMomentum()->write( "momentum", "gnuplot", 0.05 );
- getchar();*/
+ getchar();
+ */
}
diff --git a/examples/transport-equation/tnl-transport-equation.h b/examples/transport-equation/tnl-transport-equation.h
index be0083153b..381a1691ae 100644
--- a/examples/transport-equation/tnl-transport-equation.h
+++ b/examples/transport-equation/tnl-transport-equation.h
@@ -14,6 +14,7 @@
#include <TNL/Operators/DirichletBoundaryConditions.h>
#include <TNL/Operators/NeumannBoundaryConditions.h>
#include <TNL/Operators/Advection/LaxFridrichs.h>
+#include <TNL/Operators/Advection/Upwind.h>
#include <TNL/Functions/Analytic/Constant.h>
#include <TNL/Functions/VectorField.h>
#include <TNL/Meshes/Grid.h>
@@ -51,6 +52,9 @@ template< typename ConfigTag >class advectionConfig
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< String >( "differential-operator-type", "Choose the differential operator type.", "lax-friedrichs");
+ config.addEntryEnum< String >( "lax-friedrichs" );
+ config.addEntryEnum< String >( "upwind" );
}
};
@@ -101,6 +105,11 @@ class advectionSetter
template< typename VelocityFieldType >
static bool setDifferentialOperatorType( const Config::ParameterContainer& parameters )
{
+ String differentialOperatorType = parameters.getParameter< String >( "differential-operator-type" );
+ if( differentialOperatorType == "upwind" )
+ {
+ typedef Operators::Advection::Upwind< MeshType, Real, Index, VelocityFieldType > DifferentialOperatorType;
+ }
typedef Operators::Advection::LaxFridrichs< MeshType, Real, Index, VelocityFieldType > DifferentialOperatorType;
return setBoundaryConditionsType< DifferentialOperatorType >( parameters );
}
diff --git a/src/TNL/Operators/Advection/CMakeLists.txt b/src/TNL/Operators/Advection/CMakeLists.txt
index 5dff0bfc3d..ef489bf38e 100644
--- a/src/TNL/Operators/Advection/CMakeLists.txt
+++ b/src/TNL/Operators/Advection/CMakeLists.txt
@@ -1,4 +1,5 @@
-SET( headers LaxFridrichs.h )
+SET( headers LaxFridrichs.h
+ Upwind.h )
SET( CURRENT_DIR ${CMAKE_SOURCE_DIR}/src/TNL/Operators/Advection )
--
GitLab