From 61cb4617aa352db6b2ddc2d0dac1be4f35a5dd43 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Vladim=C3=ADr=20Klement?= <wlada@post.cz>
Date: Mon, 16 Feb 2015 12:15:53 +0100
Subject: [PATCH] To co jsme tady udelali
---
.../tnlLinearDiffusion.h | 168 +++++++++++
.../tnlLinearDiffusion_impl.h | 272 ++++++++++++++++++
.../tnlNSProblem.h | 107 +++++++
.../tnlNSProblem_impl.h | 247 ++++++++++++++++
4 files changed, 794 insertions(+)
create mode 100644 examples/incompressible-navier-stokes/tnlLinearDiffusion.h
create mode 100644 examples/incompressible-navier-stokes/tnlLinearDiffusion_impl.h
create mode 100644 examples/incompressible-navier-stokes/tnlNSProblem.h
create mode 100644 examples/incompressible-navier-stokes/tnlNSProblem_impl.h
diff --git a/examples/incompressible-navier-stokes/tnlLinearDiffusion.h b/examples/incompressible-navier-stokes/tnlLinearDiffusion.h
new file mode 100644
index 0000000000..35eebca181
--- /dev/null
+++ b/examples/incompressible-navier-stokes/tnlLinearDiffusion.h
@@ -0,0 +1,168 @@
+#ifndef TNLLINEARDIFFUSION_H
+#define TNLLINEARDIFFUSION_H
+
+#include <core/vectors/tnlVector.h>
+#include <mesh/tnlGrid.h>
+
+template< typename Mesh,
+ typename Real = typename Mesh::RealType,
+ typename Index = typename Mesh::IndexType >
+class tnlLinearDiffusion
+{
+
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class tnlLinearDiffusion< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index >
+{
+ public:
+
+ typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const RealType& time ) const;
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const;
+
+ template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const;
+
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >
+{
+ public:
+
+ typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const;
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const;
+
+ template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const;
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >
+{
+ public:
+
+ typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+
+ static tnlString getType();
+
+ template< typename Vector >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Real getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const;
+
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ Index getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const;
+
+ template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+ void updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const;
+
+};
+
+
+#include <operators/diffusion/tnlLinearDiffusion_impl.h>
+
+
+#endif /* TNLLINEARDIFFUSION_H */
diff --git a/examples/incompressible-navier-stokes/tnlLinearDiffusion_impl.h b/examples/incompressible-navier-stokes/tnlLinearDiffusion_impl.h
new file mode 100644
index 0000000000..4838978478
--- /dev/null
+++ b/examples/incompressible-navier-stokes/tnlLinearDiffusion_impl.h
@@ -0,0 +1,272 @@
+
+#ifndef TNLLINEARDIFFUSION_IMP_H
+#define TNLLINEARDIFFUSION_IMP_H
+
+#include <operators/diffusion/tnlLinearDiffusion.h>
+#include <mesh/tnlGrid.h>
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString
+tnlLinearDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+ return tnlString( "tnlLinearDiffusion< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlLinearDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return ( u[ mesh.getCellXPredecessor( cellIndex ) ]
+ - 2.0 * u[ cellIndex ]
+ + u[ mesh.getCellXSuccessor( cellIndex ) ] ) * mesh.getHxSquareInverse();
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+Index
+tnlLinearDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
+getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const
+{
+ return 3;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+ template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlLinearDiffusion< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::
+updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const
+{
+ const RealType lambdaX = tau * mesh.getHxSquareInverse();
+ //printf( "tau = %f lambda = %f dx_sqr = %f dx = %f, \n", tau, lambdaX, mesh.getHxSquareInverse(), mesh.getHx() );
+ matrixRow.setElement( 0, mesh.getCellXPredecessor( index ), - lambdaX );
+ matrixRow.setElement( 1, index, 2.0 * lambdaX );
+ matrixRow.setElement( 2, mesh.getCellXSuccessor( index ), - lambdaX );
+ //printf( "Linear diffusion index %d columns %d %d %d \n", index, columns[ 0 ], columns[ 1 ], columns[ 2 ] );
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString
+tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+ return tnlString( "tnlLinearDiffusion< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+Index
+tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const
+{
+ return 5;
+}
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return ( u[ mesh.getCellXPredecessor( cellIndex ) ]
+ - 2.0 * u[ cellIndex ]
+ + u[ mesh.getCellXSuccessor( cellIndex ) ] ) * mesh.getHxSquareInverse() +
+ ( u[ mesh.getCellYPredecessor( cellIndex ) ]
+ - 2.0 * u[ cellIndex ]
+ + u[ mesh.getCellYSuccessor( cellIndex ) ] ) * mesh.getHySquareInverse();
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+ template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlLinearDiffusion< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::
+updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const
+{
+ const RealType lambdaX = tau * mesh.getHxSquareInverse();
+ const RealType lambdaY = tau * mesh.getHySquareInverse();
+ matrixRow.setElement( 0, mesh.getCellYPredecessor( index ), -lambdaY );
+ matrixRow.setElement( 1, mesh.getCellXPredecessor( index ), -lambdaX );
+ matrixRow.setElement( 2, index, 2.0 * ( lambdaX + lambdaY ) );
+ matrixRow.setElement( 3, mesh.getCellXSuccessor( index ), -lambdaX );
+ matrixRow.setElement( 4, mesh.getCellYSuccessor( index ), -lambdaY );
+}
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString
+tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
+getType()
+{
+ return tnlString( "tnlLinearDiffusion< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + " >";
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real
+tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
+getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time ) const
+{
+ return ( u[ mesh.getCellXPredecessor( cellIndex ) ]
+ - 2.0 * u[ cellIndex ]
+ + u[ mesh.getCellXSuccessor( cellIndex ) ] ) * mesh.getHxSquareInverse() +
+ ( u[ mesh.getCellYPredecessor( cellIndex ) ]
+ - 2.0 * u[ cellIndex ]
+ + u[ mesh.getCellYSuccessor( cellIndex ) ] ) * mesh.getHySquareInverse() +
+ ( u[ mesh.getCellZPredecessor( cellIndex ) ]
+ - 2.0 * u[ cellIndex ]
+ + u[ mesh.getCellZSuccessor( cellIndex ) ] ) * mesh.getHzSquareInverse();
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+#ifdef HAVE_CUDA
+ __device__ __host__
+#endif
+Index
+tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
+getLinearSystemRowLength( const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates ) const
+{
+ return 7;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+ template< typename Vector, typename MatrixRow >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+void
+tnlLinearDiffusion< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::
+updateLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ const IndexType& index,
+ const CoordinatesType& coordinates,
+ Vector& u,
+ Vector& b,
+ MatrixRow& matrixRow ) const
+{
+ const RealType lambdaX = tau * mesh.getHxSquareInverse();
+ const RealType lambdaY = tau * mesh.getHySquareInverse();
+ const RealType lambdaZ = tau * mesh.getHzSquareInverse();
+ matrixRow.setElement( 0, mesh.getCellZPredecessor( index ), -lambdaZ );
+ matrixRow.setElement( 1, mesh.getCellYPredecessor( index ), -lambdaY );
+ matrixRow.setElement( 2, mesh.getCellXPredecessor( index ), -lambdaX );
+ matrixRow.setElement( 3, index, 2.0 * ( lambdaX + lambdaY + lambdaZ ) );
+ matrixRow.setElement( 4, mesh.getCellXSuccessor( index ), -lambdaX );
+ matrixRow.setElement( 5, mesh.getCellYSuccessor( index ), -lambdaY );
+ matrixRow.setElement( 6, mesh.getCellZSuccessor( index ), -lambdaZ );
+}
+
+
+
+#endif /* TNLLINEARDIFFUSION_IMP_H */
diff --git a/examples/incompressible-navier-stokes/tnlNSProblem.h b/examples/incompressible-navier-stokes/tnlNSProblem.h
new file mode 100644
index 0000000000..6014a754d7
--- /dev/null
+++ b/examples/incompressible-navier-stokes/tnlNSProblem.h
@@ -0,0 +1,107 @@
+/***************************************************************************
+ tnlHeatEquationProblem.h - description
+ -------------------
+ begin : Feb 23, 2013
+ copyright : (C) 2013 by Tomas Oberhuber
+ email : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ ***************************************************************************/
+
+#ifndef TNLHEATEQUATIONPROBLEM_H_
+#define TNLHEATEQUATIONPROBLEM_H_
+
+#include <problems/tnlPDEProblem.h>
+#include <operators/diffusion/tnlLinearDiffusion.h>
+#include <array/tnlStaticArray.h>
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator = tnlLinearDiffusion< Mesh,
+ typename BoundaryCondition::RealType > >
+class tnlNSProblem : public tnlPDEProblem< Mesh,
+ typename DifferentialOperator::RealType,
+ typename Mesh::DeviceType,
+ typename DifferentialOperator::IndexType >
+{
+ public:
+
+ typedef typename DifferentialOperator::RealType RealType;
+ typedef typename Mesh::DeviceType DeviceType;
+ typedef typename DifferentialOperator::IndexType IndexType;
+ typedef tnlPDEProblem< Mesh, RealType, DeviceType, IndexType > BaseType;
+
+ using typename BaseType::MeshType;
+ using typename BaseType::DofVectorType;
+
+ enum { Dimensions = Mesh::Dimensions };
+
+ static tnlString getTypeStatic() {return tnlString( "tnlNSProblem< " ) + Mesh :: getTypeStatic() + " >";}
+
+ tnlString getPrologHeader() const{return tnlString( "NS equation" );}
+
+ void writeProlog( tnlLogger& logger,
+ const tnlParameterContainer& parameters ) const {}
+
+ bool setup( const tnlParameterContainer& parameters );
+
+ bool setInitialCondition( const tnlParameterContainer& parameters,
+ const MeshType& mesh,
+ DofVectorType& dofs,
+ DofVectorType& auxDofs );
+
+ template< typename MatrixType >
+ bool setupLinearSystem( const MeshType& mesh,
+ MatrixType& matrix );
+
+ bool makeSnapshot( const RealType& time,
+ const IndexType& step,
+ const MeshType& mesh,
+ DofVectorType& dofs,
+ DofVectorType& auxDofs );
+
+ IndexType getDofs( const MeshType& mesh ) const;
+
+ void bindDofs( const MeshType& mesh,
+ DofVectorType& dofs );
+
+ void getExplicitRHS( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ DofVectorType& _u,
+ DofVectorType& _fu );
+
+ template< typename MatrixType >
+ void assemblyLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ DofVectorType& dofs,
+ DofVectorType& auxDofs,
+ MatrixType& matrix,
+ DofVectorType& rightHandSide );
+
+
+ protected:
+
+ tnlSharedVector< RealType, DeviceType, IndexType > p;
+
+ tnlStaticArray< Dimensions, tnlSharedVector< RealType, DeviceType, IndexType > > u;
+
+ DifferentialOperator differentialOperator;
+
+ BoundaryCondition boundaryCondition;
+
+ RightHandSide rightHandSide;
+};
+
+#include <problems/tnlHeatEquationProblem_impl.h>
+
+#endif /* TNLHEATEQUATIONPROBLEM_H_ */
diff --git a/examples/incompressible-navier-stokes/tnlNSProblem_impl.h b/examples/incompressible-navier-stokes/tnlNSProblem_impl.h
new file mode 100644
index 0000000000..1050c0ffd8
--- /dev/null
+++ b/examples/incompressible-navier-stokes/tnlNSProblem_impl.h
@@ -0,0 +1,247 @@
+/***************************************************************************
+ tnlHeatEquationProblem_impl.h - description
+ -------------------
+ begin : Mar 10, 2013
+ copyright : (C) 2013 by Tomas Oberhuber
+ email : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ ***************************************************************************/
+
+#ifndef TNLHEATEQUATIONPROBLEM_IMPL_H_
+#define TNLHEATEQUATIONPROBLEM_IMPL_H_
+
+#include <core/mfilename.h>
+#include <matrices/tnlMatrixSetter.h>
+#include <matrices/tnlMultidiagonalMatrixSetter.h>
+#include <core/tnlLogger.h>
+#include <solvers/pde/tnlExplicitUpdater.h>
+#include <solvers/pde/tnlLinearSystemAssembler.h>
+
+
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+bool
+tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+setup( const tnlParameterContainer& parameters )
+{
+ if( ! this->boundaryCondition.setup( parameters, "boundary-conditions-" ) ||
+ ! this->rightHandSide.setup( parameters, "right-hand-side-" ) )
+ return false;
+ return true;
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+typename tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::IndexType
+tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+getDofs( const MeshType& mesh ) const
+{
+ /****
+ * Set-up DOFs and supporting grid functions
+ */
+ return mesh.getNumberOfFaces();
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+typename tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::IndexType
+tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+getAuxiliaryDofs( const MeshType& mesh ) const
+{
+ /****
+ * Set-up DOFs and supporting grid functions
+ */
+ return mesh.getNumberOfCells();
+}
+
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+void
+tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+bindDofs( const MeshType& mesh,
+ DofVectorType& dofVector )
+{
+ const IndexType uDofs = mesh.getNumberOfFaces() / 2;
+ this->u[ 0 ].bind( dofVector.getData(), uDofs );
+ this->u[ 1 ].bind( &dofVector.getData()[ uDofs ], uDofs );
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+void
+tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+bindAuxiliaryDofs( const MeshType& mesh,
+ DofVectorType& auxDofs )
+{
+ const IndexType pDofs = mesh.getNumberOfCells();
+ this->p.bind( auxDofs.getData(), pDofs );
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+bool
+tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+setInitialCondition( const tnlParameterContainer& parameters,
+ const MeshType& mesh,
+ DofVectorType& dofs,
+ DofVectorType& auxiliaryDofs )
+{
+ this->bindDofs( mesh, dofs );
+ this->bindAuxiliaryDofs( mesh, dofs );
+ this->p.setValue( 0.0 );
+ this->u[ 0 ].setValue( 0.0 );
+ this->u[ 1 ].setValue( 0.0 );
+
+ /*const tnlString& initialConditionFile = parameters.getParameter< tnlString >( "initial-condition" );
+ if( ! this->solution.load( initialConditionFile ) )
+ {
+ cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." << endl;
+ return false;
+ }*/
+ return true;
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+ template< typename MatrixType >
+bool
+tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+setupLinearSystem( const MeshType& mesh,
+ MatrixType& matrix )
+{
+ const IndexType dofs = this->getDofs( mesh );
+ typedef typename MatrixType::RowLengthsVector RowLengthsVectorType;
+ RowLengthsVectorType rowLengths;
+ if( ! rowLengths.setSize( dofs ) )
+ return false;
+ tnlMatrixSetter< MeshType, DifferentialOperator, BoundaryCondition, RowLengthsVectorType > matrixSetter;
+ matrixSetter.template getRowLengths< Mesh::Dimensions >( mesh,
+ differentialOperator,
+ boundaryCondition,
+ rowLengths );
+ matrix.setDimensions( dofs, dofs );
+ if( ! matrix.setRowLengths( rowLengths ) )
+ return false;
+ return true;
+ //return tnlMultidiagonalMatrixSetter< Mesh >::setupMatrix( mesh, matrix );
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+bool
+tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+makeSnapshot( const RealType& time,
+ const IndexType& step,
+ const MeshType& mesh,
+ DofVectorType& dofs,
+ DofVectorType& auxiliaryDofs )
+{
+ cout << endl << "Writing output at time " << time << " step " << step << "." << endl;
+
+ this->bindDofs( mesh, dofs );
+ this->bindAuxiliaryDofs( mesh, auxiliaryDofs );
+ //cout << "dofs = " << dofs << endl;
+ tnlString fileName;
+ FileNameBaseNumberEnding( "u-", step, 5, ".tnl", fileName );
+ if( ! this->solution.save( fileName ) )
+ return false;
+ return true;
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+void
+tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+getExplicitRHS( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ DofVectorType& u,
+ DofVectorType& fu )
+{
+ /****
+ * If you use an explicit solver like tnlEulerSolver or tnlMersonSolver, you
+ * need to implement this method. Compute the right-hand side of
+ *
+ * d/dt u(x) = fu( x, u )
+ *
+ * You may use supporting vectors again if you need.
+ */
+
+ //cout << "u = " << u << endl;
+ this->bindDofs( mesh, u );
+
+ tnlExplicitUpdater< Mesh, DofVectorType, DifferentialOperator, BoundaryCondition, RightHandSide > explicitUpdater;
+ explicitUpdater.template update< Mesh::Dimensions >( time,
+ mesh,
+ this->differentialOperator,
+ this->boundaryCondition,
+ this->rightHandSide,
+ u,
+ fu );
+ //cout << "u = " << u << endl;
+ //cout << "fu = " << fu << endl;
+ //_u.save( "u.tnl" );
+ //_fu.save( "fu.tnl" );
+ //getchar();
+}
+
+template< typename Mesh,
+ typename BoundaryCondition,
+ typename RightHandSide,
+ typename DifferentialOperator >
+ template< typename MatrixType >
+void
+tnlNSProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
+assemblyLinearSystem( const RealType& time,
+ const RealType& tau,
+ const MeshType& mesh,
+ DofVectorType& u,
+ DofVectorType& auxDofs,
+ MatrixType& matrix,
+ DofVectorType& b )
+{
+ tnlLinearSystemAssembler< Mesh, DofVectorType, DifferentialOperator, BoundaryCondition, RightHandSide, MatrixType > systemAssembler;
+ systemAssembler.template assembly< Mesh::Dimensions >( time,
+ tau,
+ mesh,
+ this->differentialOperator,
+ this->boundaryCondition,
+ this->rightHandSide,
+ u,
+ matrix,
+ b );
+ /*matrix.print( cout );
+ cout << endl << b << endl;
+ cout << endl << u << endl;
+ abort();*/
+}
+
+#endif /* TNLHEATEQUATIONPROBLEM_IMPL_H_ */
--
GitLab