Commit 8110508d authored by Tomáš Oberhuber's avatar Tomáš Oberhuber
Browse files

The transport equation is refactored.

parent edfebfb0

examples/advection/advection.h

+54 −48
Original line number Diff line number Diff line
@@ -3,12 +3,11 @@ 
#include <TNL/Solvers/BuildConfigTags.h>

#include <TNL/Operators/DirichletBoundaryConditions.h>

#include <TNL/Operators/NeumannBoundaryConditions.h>

#include <TNL/Operators/Advection/LaxFridrichs.h>

#include <TNL/Functions/Analytic/Constant.h>

#include <TNL/Functions/VectorField.h>

#include <TNL/Meshes/Grid.h>

#include "advectionProblem.h"

#include "LaxFridrichs.h"

#include "advectionRhs.h"

#include "advectionBuildConfigTag.h"



using namespace TNL;
@@ -16,11 +15,11 @@ using namespace TNL; 
typedef advectionBuildConfigTag BuildConfig;



/****

 * Uncoment the following (and comment the previous line) for the complete build.

 * Uncomment the following (and comment the previous line) for the complete build.

 * This will include support for all floating point precisions, all indexing types

 * and more solvers. You may then choose between them from the command line.

 * The compile time may, however, take tens of minutes or even several hours,

 * esppecially if CUDA is enabled. Use this, if you want, only for the final build,

 * especially if CUDA is enabled. Use this, if you want, only for the final build,

 * not in the development phase.

 */

//typedef tnlDefaultConfigTag BuildConfig;
@@ -33,11 +32,10 @@ template< typename ConfigTag >class advectionConfig 
         config.addDelimiter( "Advection settings:" );

         config.addEntry< String >( "velocity-field", "Type of velocity field.", "constant" );

            config.addEntryEnum< String >( "constant" );

            //config.addEntryEnum< String >( "file" );

         Functions::VectorField< 3, Functions::Analytic::Constant< 3 > >::configSetup( config, "velocity-field-" );

         

         typedef Meshes::Grid< 3 > MeshType;

         LaxFridrichs< MeshType >::configSetup( config, "lax-fridrichs" );

         Operators::Advection::LaxFridrichs< MeshType >::configSetup( config );

         

         config.addEntry< String >( "boundary-conditions-type", "Choose the boundary conditions type.", "dirichlet");

            config.addEntryEnum< String >( "dirichlet" );
@@ -60,50 +58,58 @@ class advectionSetter 
      typedef Device DeviceType;

      typedef Index IndexType;

      

      static bool run( const Config::ParameterContainer & parameters )

      {

      static const int Dimensions = MeshType::getMeshDimensions();

          typedef Functions::Analytic::Constant< Dimensions, RealType > ConstantFunctionType;

          typedef Functions::VectorField< Dimensions, ConstantFunctionType > VelocityFieldType;

          typedef LaxFridrichs< MeshType, Real, Index, VelocityFieldType > ApproximateOperator;

          typedef advectionRhs< MeshType, Real > RightHandSide;

          typedef Containers::StaticVector < MeshType::getMeshDimensions(), Real > Vertex;

      

         /****

          * Resolve the template arguments of your solver here.

          * The following code is for the Dirichlet and the Neumann boundary conditions.

          * Both can be constant or defined as descrete values of Vector.

          */

          String boundaryConditionsType = parameters.getParameter< String >( "boundary-conditions-type" );

          if( parameters.checkParameter( "boundary-conditions-constant" ) )

          {

             typedef Functions::Analytic::Constant< Dimensions, Real > Constant;

             if( boundaryConditionsType == "dirichlet" )

      template< typename Problem >

      static bool callSolverStarter( const Config::ParameterContainer& parameters )

      {

                typedef Operators::DirichletBoundaryConditions< MeshType, Constant, MeshType::getMeshDimensions(), Real, Index > BoundaryConditions;

                typedef advectionProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;

         SolverStarter solverStarter;

         return solverStarter.template run< Problem >( parameters );

      }

             typedef Operators::NeumannBoundaryConditions< MeshType, Constant, Real, Index > BoundaryConditions;

             typedef advectionProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;

             SolverStarter solverStarter;

             return solverStarter.template run< Problem >( parameters );

          }

          typedef Functions::MeshFunction< MeshType > MeshFunction;

      

      template< typename DifferentialOperatorType >

      static bool setBoundaryConditionsType( const Config::ParameterContainer& parameters )

      {

         typedef Functions::Analytic::Constant< Dimensions, Real > ConstantFunctionType;

         String boundaryConditionsType = parameters.getParameter< String >( "boundary-conditions-type" );

         if( boundaryConditionsType == "dirichlet" )

         {

             typedef Operators::DirichletBoundaryConditions< MeshType, MeshFunction, MeshType::getMeshDimensions(), Real, Index > BoundaryConditions;

             typedef advectionProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;

             SolverStarter solverStarter;

             return solverStarter.template run< Problem >( parameters );

            typedef Operators::DirichletBoundaryConditions< MeshType, ConstantFunctionType, MeshType::getMeshDimensions(), Real, Index > BoundaryConditions;

            typedef advectionProblem< MeshType, BoundaryConditions, ConstantFunctionType, DifferentialOperatorType > Problem;

            return callSolverStarter< Problem >( parameters );

         }

          typedef Operators::NeumannBoundaryConditions< MeshType, MeshFunction, Real, Index > BoundaryConditions;

          typedef advectionProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;

          SolverStarter solverStarter;

          return solverStarter.template run< Problem >( parameters );

         if( boundaryConditionsType == "neumann" )

         {

            typedef Operators::DirichletBoundaryConditions< MeshType, ConstantFunctionType, MeshType::getMeshDimensions(), Real, Index > BoundaryConditions;

            typedef advectionProblem< MeshType, BoundaryConditions, ConstantFunctionType, DifferentialOperatorType > Problem;

            return callSolverStarter< Problem >( parameters );

         }

         std::cerr << "Unknown boundary conditions type: " << boundaryConditionsType << "." << std::endl;

         return false;

      }

      

      template< typename VelocityFieldType >

      static bool setDifferentialOperatorType( const Config::ParameterContainer& parameters )

      {

         typedef Operators::Advection::LaxFridrichs< MeshType, Real, Index, VelocityFieldType > DifferentialOperatorType;

         return setBoundaryConditionsType< DifferentialOperatorType >( parameters );

      }

      

      static bool setVelocityFieldType( const Config::ParameterContainer& parameters )

      {

         String velocityFieldType = parameters.getParameter< String >( "velocity-field" );

         if( velocityFieldType == "constant" )

         {

            typedef Functions::Analytic::Constant< Dimensions, RealType > VelocityFieldType;

            return setDifferentialOperatorType< VelocityFieldType >( parameters );

         }

         return false;

      }



      static bool run( const Config::ParameterContainer& parameters )

      {

         return setVelocityFieldType( parameters );

      }      

};



int main( int argc, char* argv[] )

examples/advection/advectionProblem.h

+2 −3
Original line number Diff line number Diff line
@@ -12,8 +12,7 @@ namespace TNL { 
template< typename Mesh,

          typename BoundaryCondition,

          typename RightHandSide,

          typename DifferentialOperator,

          typename VelocityField = Functions::MeshFunction< Mesh > >

          typename DifferentialOperator >

class advectionProblem:

public PDEProblem< Mesh,

                   typename DifferentialOperator::RealType,
@@ -31,7 +30,7 @@ public PDEProblem< Mesh, 
      typedef SharedPointer< DifferentialOperator > DifferentialOperatorPointer;

      typedef SharedPointer< BoundaryCondition > BoundaryConditionPointer;

      typedef SharedPointer< RightHandSide, DeviceType > RightHandSidePointer;

      typedef VelocityField VelocityFieldType;

      typedef typename DifferentialOperator::VelocityFieldType VelocityFieldType;

      typedef SharedPointer< VelocityFieldType, DeviceType > VelocityFieldPointer;



      using typename BaseType::MeshType;

examples/advection/advectionProblem_impl.h

+28 −45
Original line number Diff line number Diff line
@@ -7,15 +7,16 @@ 
#include <TNL/Solvers/PDE/LinearSystemAssembler.h>

#include <TNL/Solvers/PDE/BackwardTimeDiscretisation.h>



#include "advectionProblem.h"



namespace TNL {



template< typename Mesh,

          typename BoundaryCondition,

          typename RightHandSide,

          typename DifferentialOperator,

          typename VelocityField >

          typename DifferentialOperator >

String

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator, VelocityField >::

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::

getTypeStatic()

{

   return String( "advectionProblem< " ) + Mesh :: getTypeStatic() + " >";
@@ -24,10 +25,9 @@ getTypeStatic() 
template< typename Mesh,

          typename BoundaryCondition,

          typename RightHandSide,

          typename DifferentialOperator,

          typename VelocityField >

          typename DifferentialOperator >

String

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator, VelocityField >::

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::

getPrologHeader() const

{

   return String( "advection" );
@@ -36,10 +36,9 @@ getPrologHeader() const 
template< typename Mesh,

          typename BoundaryCondition,

          typename RightHandSide,

          typename DifferentialOperator,

          typename VelocityField >

          typename DifferentialOperator >

void

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator, VelocityField >::

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::

writeProlog( Logger& logger, const Config::ParameterContainer& parameters ) const

{

   /****
@@ -51,26 +50,16 @@ writeProlog( Logger& logger, const Config::ParameterContainer& parameters ) cons 
template< typename Mesh,

          typename BoundaryCondition,

          typename RightHandSide,

          typename DifferentialOperator,

          typename VelocityField >

          typename DifferentialOperator >

bool

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator, VelocityField >::

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::

setup( const MeshPointer& meshPointer,

       const Config::ParameterContainer& parameters,

       const String& prefix )

{

   //this->velocityType = parameters.getParameter< String >( "move" );

   /*const double artificalViscosity = parameters.getParameter< double >( "artifical-viscosity" );

   differentialOperatorPointer->setViscosity(artificalViscosity);

   const double advectionSpeedX = parameters.getParameter< double >( "advection-speedX" );

   differentialOperatorPointer->setAdvectionSpeedX(advectionSpeedX);

   const double advectionSpeedY = parameters.getParameter< double >( "advection-speedY" );

   differentialOperatorPointer->setAdvectionSpeedY(advectionSpeedY);*/

   

   if( ! this->velocityField->setup( meshPointer, parameters, prefix + "velocity-field-" ) ||

       ! this->differentialOperatorPointer->setup( meshPointer, parameters, prefix + "advection-" ) ||

       ! this->boundaryConditionPointer->setup( meshPointer, parameters, prefix + "boundary-conditions-" ) ||

       ! this->rightHandSidePointer->setup( parameters, prefix + "right-hand-side-" ) )

       ! this->differentialOperatorPointer->setup( meshPointer, parameters, prefix ) ||

       ! this->boundaryConditionPointer->setup( meshPointer, parameters, prefix + "boundary-conditions-" ) )

      return false;

   return true;

}
@@ -78,10 +67,9 @@ setup( const MeshPointer& meshPointer, 
template< typename Mesh,

          typename BoundaryCondition,

          typename RightHandSide,

          typename DifferentialOperator,

          typename VelocityField >

typename advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator, VelocityField >::IndexType

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator, VelocityField >::

          typename DifferentialOperator >

typename advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::IndexType

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::

getDofs( const MeshPointer& mesh ) const

{

   /****
@@ -94,10 +82,9 @@ getDofs( const MeshPointer& mesh ) const 
template< typename Mesh,

          typename BoundaryCondition,

          typename RightHandSide,

          typename DifferentialOperator,

          typename VelocityField >

          typename DifferentialOperator >

void

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator, VelocityField >::

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::

bindDofs( const MeshPointer& meshPointer,

          DofVectorPointer& dofVector )

{
@@ -108,10 +95,9 @@ bindDofs( const MeshPointer& meshPointer, 
template< typename Mesh,

          typename BoundaryCondition,

          typename RightHandSide,

          typename DifferentialOperator,

          typename VelocityField >

          typename DifferentialOperator >

bool

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator, VelocityField >::

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::

setInitialCondition( const Config::ParameterContainer& parameters,

                     const MeshPointer& meshPointer,

                     DofVectorPointer& dofs,
@@ -130,11 +116,10 @@ setInitialCondition( const Config::ParameterContainer& parameters, 
template< typename Mesh,

          typename BoundaryCondition,

          typename RightHandSide,

          typename DifferentialOperator,

          typename VelocityField >

          typename DifferentialOperator >

   template< typename Matrix >

bool

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator, VelocityField >::

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::

setupLinearSystem( const MeshPointer& mesh,

                   Matrix& matrix )

{
@@ -157,10 +142,9 @@ setupLinearSystem( const MeshPointer& mesh, 
template< typename Mesh,

          typename BoundaryCondition,

          typename RightHandSide,

          typename DifferentialOperator,

          typename VelocityField >

          typename DifferentialOperator >

bool

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator, VelocityField >::

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::

makeSnapshot( const RealType& time,

              const IndexType& step,

              const MeshPointer& mesh,
@@ -181,10 +165,9 @@ makeSnapshot( const RealType& time, 
template< typename Mesh,

          typename BoundaryCondition,

          typename RightHandSide,

          typename DifferentialOperator,

          typename VelocityField >

          typename DifferentialOperator >

void

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator, VelocityField >::

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::

getExplicitRHS( const RealType& time,

                const RealType& tau,

                const MeshPointer& mesh,
@@ -207,6 +190,7 @@ getExplicitRHS( const RealType& time, 
   SharedPointer< MeshFunctionType > u( mesh, _u ); 

   SharedPointer< MeshFunctionType > fu( mesh, _fu );

   differentialOperatorPointer->setTau(tau); 

   differentialOperatorPointer->setVelocityField( this->velocityField );

   explicitUpdater.template update< typename Mesh::Cell >( time,

                                                           mesh,

                                                           this->differentialOperatorPointer,
@@ -223,11 +207,10 @@ getExplicitRHS( const RealType& time, 
template< typename Mesh,

          typename BoundaryCondition,

          typename RightHandSide,

          typename DifferentialOperator,

          typename VelocityField >

          typename DifferentialOperator >

   template< typename Matrix >

void

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator, VelocityField >::

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::

assemblyLinearSystem( const RealType& time,

                      const RealType& tau,

                      const MeshPointer& mesh,

examples/advection/advectionRhs.h

deleted100644 → 0
+0 −33
Original line number Diff line number Diff line 
#pragma once



#include <TNL/Functions/Domain.h>



namespace TNL {



template< typename Mesh, typename Real >class advectionRhs

  : public Functions::Domain< Mesh::meshDimensions, Functions::MeshDomain > 

{

   public:



      typedef Mesh MeshType;

      typedef Real RealType;



      bool setup( const Config::ParameterContainer& parameters,

                  const String& prefix = "" )

      {

         return true;

      }



      template< typename MeshEntity >

      __cuda_callable__

      Real operator()( const MeshEntity& entity,

                       const Real& time = 0.0 ) const

      {

         typedef typename MeshEntity::MeshType::VertexType VertexType;

         VertexType v = entity.getCenter();

         return 0.0;

      }

};



} // namespace TNL

src/TNL/Config/ParameterContainer.h

+4 −0
Original line number Diff line number Diff line
@@ -14,6 +14,7 @@ 
#include <TNL/Config/ConfigDescription.h>

#include <TNL/mpi-supp.h>

#include <TNL/param-types.h>

//#include <TNL/Debugging/StackBacktrace.h>



namespace TNL {

namespace Config {
@@ -73,7 +74,10 @@ class ParameterContainer 
            return true;

         }

      if( verbose )

      {

         std::cerr << "Missing parameter '" << name << "'." << std::endl;

         throw(0); //PrintStackBacktrace;

      }

      return false;

   }