From 1c0871189526291fab7b48467b2942c6b77f6727 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Vladim=C3=ADr=20Klement?= <wlada@post.cz>
Date: Thu, 3 Sep 2015 14:25:52 +0200
Subject: [PATCH] Committing explicist NS solver
---
examples/heat-equation/tnl-run-heat-equation | 0
.../CMakeLists.txt | 13 +-
.../tnl-incompressible-navier-stokes.h | 6 +-
.../tnlExplicitINSTimeStepper.h | 2 +-
.../tnlExplicitINSTimeStepper_impl.h | 37 ++-
.../tnlIncompressibleNavierStokesProblem.h | 193 ++++++++++++-
...nlIncompressibleNavierStokesProblem_impl.h | 62 +++-
.../tnlLinearDiffusion.h | 168 -----------
.../tnlLinearDiffusion_impl.h | 272 ------------------
install | 22 +-
.../linear/tnlLinearResidueGetter_impl.h | 3 +-
11 files changed, 294 insertions(+), 484 deletions(-)
mode change 100644 => 100755 examples/heat-equation/tnl-run-heat-equation
delete mode 100644 examples/incompressible-navier-stokes/tnlLinearDiffusion.h
delete mode 100644 examples/incompressible-navier-stokes/tnlLinearDiffusion_impl.h
diff --git a/examples/heat-equation/tnl-run-heat-equation b/examples/heat-equation/tnl-run-heat-equation
old mode 100644
new mode 100755
diff --git a/examples/incompressible-navier-stokes/CMakeLists.txt b/examples/incompressible-navier-stokes/CMakeLists.txt
index c960e25453..3a37066441 100755
--- a/examples/incompressible-navier-stokes/CMakeLists.txt
+++ b/examples/incompressible-navier-stokes/CMakeLists.txt
@@ -1,12 +1,23 @@
set( tnl_incompressible_navier_stokes_SOURCES
tnl-incompressible-navier-stokes.cpp
+tnlExplicitINSTimeStepper_impl.h
+tnlExplicitINSTimeStepper.h
+tnlPoissonProblem_impl.h
+tnlPoissonProblem.h
+tnlINSRightHandSide.h
+tnlINSBoundaryConditions.h
+tnlIncompressibleNavierStokesProblem_impl.h
+tnlIncompressibleNavierStokesProblem.h
+tnlNSFastBuildConfig.h
+visit_writer.h
+visit_writer.cpp
)
IF( BUILD_CUDA )
CUDA_ADD_EXECUTABLE(tnl-incompressible-navier-stokes${debugExt} tnl-incompressible-navier-stokes.cu)
ELSE( BUILD_CUDA )
- ADD_EXECUTABLE(tnl-incompressible-navier-stokes${debugExt} tnl-incompressible-navier-stokes.cpp)
+ ADD_EXECUTABLE(tnl-incompressible-navier-stokes${debugExt} ${tnl_incompressible_navier_stokes_SOURCES})
ENDIF( BUILD_CUDA )
diff --git a/examples/incompressible-navier-stokes/tnl-incompressible-navier-stokes.h b/examples/incompressible-navier-stokes/tnl-incompressible-navier-stokes.h
index 1526acdb95..6dcea9e87c 100644
--- a/examples/incompressible-navier-stokes/tnl-incompressible-navier-stokes.h
+++ b/examples/incompressible-navier-stokes/tnl-incompressible-navier-stokes.h
@@ -19,16 +19,16 @@
#define TNL_INCOMPRESSIBLE_NAVIER_STOKES_H_
#include <solvers/tnlSolver.h>
-#include <solvers/tnlFastBuildConfig.h>
#include <solvers/tnlConfigTags.h>
#include <operators/diffusion/tnlLinearDiffusion.h>
#include "tnlINSBoundaryConditions.h"
#include "tnlINSRightHandSide.h"
#include "tnlIncompressibleNavierStokes.h"
#include "tnlIncompressibleNavierStokesProblem.h"
+#include "tnlNSFastBuildConfig.h"
//typedef tnlDefaultConfigTag BuildConfig;
-typedef tnlFastBuildConfig BuildConfig;
+typedef tnlNSFastBuildConfig BuildConfig;
template< typename ConfigTag >
class tnlIncompressibleNavierStokesConfig
@@ -37,6 +37,8 @@ class tnlIncompressibleNavierStokesConfig
static void configSetup( tnlConfigDescription& config )
{
config.addDelimiter( "Incompressible Navier-Stokes solver settings:" );
+ config.addEntry< double >( "viscosity", "Viscosity of the diffusion." );
+ config.addEntry< double >( "inletVelocity", "Maximal X velocity on the inlet." );
/*config.addEntry< tnlString >( "boundary-conditions-type", "Choose the boundary conditions type.", "dirichlet");
config.addEntryEnum< tnlString >( "dirichlet" );
diff --git a/examples/incompressible-navier-stokes/tnlExplicitINSTimeStepper.h b/examples/incompressible-navier-stokes/tnlExplicitINSTimeStepper.h
index 8319d0bf58..f0e8593ebc 100644
--- a/examples/incompressible-navier-stokes/tnlExplicitINSTimeStepper.h
+++ b/examples/incompressible-navier-stokes/tnlExplicitINSTimeStepper.h
@@ -58,7 +58,7 @@ class tnlExplicitINSTimeStepper
protected:
Problem* problem;
-
+ //LinearSolver _matSolver;
RealType timeStep;
};
diff --git a/examples/incompressible-navier-stokes/tnlExplicitINSTimeStepper_impl.h b/examples/incompressible-navier-stokes/tnlExplicitINSTimeStepper_impl.h
index 3aeaefa51b..90d2998b1d 100644
--- a/examples/incompressible-navier-stokes/tnlExplicitINSTimeStepper_impl.h
+++ b/examples/incompressible-navier-stokes/tnlExplicitINSTimeStepper_impl.h
@@ -18,6 +18,7 @@
#ifndef EXAMPLES_INCOMPRESSIBLE_NAVIER_STOKES_TNLEXPLICITINSTIMESTEPPER_IMPL_H_
#define EXAMPLES_INCOMPRESSIBLE_NAVIER_STOKES_TNLEXPLICITINSTIMESTEPPER_IMPL_H_
+#include "tnlExplicitINSTimeStepper.h"
template< typename Problem,
typename LinearSystemSolver >
tnlExplicitINSTimeStepper< Problem, LinearSystemSolver >::
@@ -32,7 +33,7 @@ template< typename Problem,
void
tnlExplicitINSTimeStepper< Problem, LinearSystemSolver >::
configSetup( tnlConfigDescription& config,
- const tnlString& prefix )
+ const tnlString& prefix ) //Pridavani parametru prikazove radky
{
config.addEntry< bool >( "verbose", "Verbose mode.", true );
}
@@ -42,7 +43,7 @@ template< typename Problem,
bool
tnlExplicitINSTimeStepper< Problem, LinearSystemSolver >::
setup( const tnlParameterContainer& parameters,
- const tnlString& prefix )
+ const tnlString& prefix ) //Nacteni parametru prikazove radky
{
//this->verbose = parameters.getParameter< bool >( "verbose" );
return true;
@@ -52,7 +53,7 @@ template< typename Problem,
typename LinearSystemSolver >
bool
tnlExplicitINSTimeStepper< Problem, LinearSystemSolver >::
-init( const MeshType& mesh )
+init( const MeshType& mesh ) //Inicializace time stepperu - vytvoreni matic podle site
{
/*cout << "Setting up the linear system...";
if( ! this->problem->setupLinearSystem( mesh, this->matrix ) )
@@ -77,7 +78,7 @@ template< typename Problem,
typename LinearSystemSolver >
void
tnlExplicitINSTimeStepper< Problem, LinearSystemSolver >::
-setProblem( ProblemType& problem )
+setProblem( ProblemType& problem ) //Nesahej
{
this -> problem = &problem;
};
@@ -114,16 +115,30 @@ solve( const RealType& time,
const RealType& stopTime,
const MeshType& mesh,
DofVectorType& dofVector,
- DofVectorType& auxiliaryDofVector )
+ DofVectorType& auxiliaryDofVector ) //Hlavni cast, kterou bude potreba menit
{
tnlAssert( this->problem != 0, );
- /*RealType t = time;
- this->linearSystemSolver->setMatrix( this->matrix );
+ RealType t = time;
+ //this->_matSolver.setMatrix(this->matrix);
while( t < stopTime )
{
+
RealType currentTau = Min( this->timeStep, stopTime - t );
+ currentTau = 0.005;
+
+ this->problem->diffuse(currentTau,mesh);
+ this->problem->project(mesh);
+ this->problem->advect(currentTau, mesh);
+ this->problem->project(mesh);
+
+ /*/if( ! this->_matSolver->template solve< DofVectorType, tnlLinearResidueGetter< typename Problem::MatrixType, DofVectorType > >( _rightHandSide, auxiliaryDofVector ) )
+ {
+ cerr << endl << "The linear system solver did not converge." << endl;
+ return false;
+ }*/
+ //this->problem->Project(auxiliaryDofVector);
- if( ! this->problem->preIterate( t,
+ /*if( ! this->problem->preIterate( t,
currentTau,
mesh,
dofVector,
@@ -143,7 +158,7 @@ solve( const RealType& time,
this->rightHandSide );
if( verbose )
cout << " Solving the linear system for time " << t << " \r" << flush;
- if( ! this->linearSystemSolver->template solve< DofVectorType, tnlLinearResidueGetter< MatrixType, DofVectorType > >( this->rightHandSide, dofVector ) )
+ if( ! this->matSolver->template solve< DofVectorType, tnlLinearResidueGetter< MatrixType, DofVectorType > >( this->rightHandSide, dofVector ) )
{
cerr << endl << "The linear system solver did not converge." << endl;
return false;
@@ -158,9 +173,9 @@ solve( const RealType& time,
{
cerr << endl << "Postiteration failed." << endl;
return false;
- }
+ }*/
t += currentTau;
- }*/
+ }
return true;
}
diff --git a/examples/incompressible-navier-stokes/tnlIncompressibleNavierStokesProblem.h b/examples/incompressible-navier-stokes/tnlIncompressibleNavierStokesProblem.h
index e45b11bbe9..41e2ae5eb3 100644
--- a/examples/incompressible-navier-stokes/tnlIncompressibleNavierStokesProblem.h
+++ b/examples/incompressible-navier-stokes/tnlIncompressibleNavierStokesProblem.h
@@ -24,6 +24,8 @@
#include "tnlINSBoundaryConditions.h"
#include "tnlExplicitINSTimeStepper.h"
+template<class T> T square(const T & val){return val*val;}
+
template< typename Mesh,
typename BoundaryCondition,
typename RightHandSide,
@@ -44,7 +46,11 @@ class tnlIncompressibleNavierStokesProblem : public tnlPDEProblem< Mesh,
using typename BaseType::MeshType;
using typename BaseType::DofVectorType;
- typedef tnlExplicitINSTimeStepper< ThisType, void > TimeStepper;
+ typedef tnlCSRMatrix< RealType, tnlHost, IndexType > MatrixType;
+ typedef tnlSORSolver<MatrixType> LinearSolver;
+ typedef tnlExplicitINSTimeStepper< ThisType, LinearSolver > TimeStepper;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+ //typedef tnlExplicitINSTimeStepper< ThisType, void > TimeStepper;
enum { Dimensions = Mesh::Dimensions };
@@ -62,7 +68,7 @@ class tnlIncompressibleNavierStokesProblem : public tnlPDEProblem< Mesh,
DofVectorType& dofs,
DofVectorType& auxDofs );
- template< typename MatrixType >
+ //template< typename MatrixType >
bool setupLinearSystem( const MeshType& mesh,
MatrixType& matrix );
@@ -83,7 +89,7 @@ class tnlIncompressibleNavierStokesProblem : public tnlPDEProblem< Mesh,
DofVectorType& _u,
DofVectorType& _fu );
- template< typename MatrixType >
+ //template< typename MatrixType >
void assemblyLinearSystem( const RealType& time,
const RealType& tau,
const MeshType& mesh,
@@ -92,12 +98,184 @@ class tnlIncompressibleNavierStokesProblem : public tnlPDEProblem< Mesh,
MatrixType& matrix,
DofVectorType& rightHandSide );
+ void set_zero_neumann(tnlSharedVector< RealType, DeviceType, IndexType > & vec)
+ {
+ /*int ex = a.width()-1, ey=a.height()-1;
+ for (int x=1; x < ex; x++) {a(x,0) = a(x,1); a(x,ey) = a(x,ey-1);}
+ for (int y=1; y < ey; y++) {a(0,y) = a(1,y); a(ex,y) = a(ex-1,y);}
+ a(0,0)=0.5*(a(0,1)+a(1,0));
+ a(ex,0)=0.5*(a(ex-1,0)+a(ex,1));
+ a(0,ey)=0.5*(a(1,ey)+a(0,ey-1));
+ a(ex,ey)=0.5*(a(ex-1,ey)+a(ex,ey-1));*/
+ }
+
+ static void SetBnd(const MeshType& mesh)
+ {
+ mesh;
+ }
+
+ void computeVelocityDivergence(IndexType cell, const tnlVector<RealType, DeviceType, IndexType> & v, const MeshType& mesh, tnlVector<RealType, DeviceType, IndexType> & rhs)
+ {
+ double diffU = v[mesh.template getFaceNextToCell<1,0>(cell)] - v[mesh.template getFaceNextToCell<-1,0>(cell)];
+ double diffV = v[mesh.template getFaceNextToCell<0,1>(cell)] - v[mesh.template getFaceNextToCell<0,-1>(cell)];
+ rhs[cell] = -0.5f*(diffU/mesh.getDimensions().x() + diffV/mesh.getDimensions().y()); // -(u_x + v_y)
+ }
+ void updateVelocityByPressureCorrection(IndexType cell, const tnlVector<RealType, DeviceType, IndexType> & v, const MeshType& mesh, tnlVector<RealType, DeviceType, IndexType> & p)
+ {
+ RealType pVal = p[cell];
+ double nx =mesh.getDimensions().x(), ny=mesh.getDimensions().y();
+ vel[mesh.template getFaceNextToCell<-1,0>(cell)] -= 0.5*nx*pVal;
+ vel[mesh.template getFaceNextToCell<+1,0>(cell)] += 0.5*nx*pVal;
+ vel[mesh.template getFaceNextToCell<0,-1>(cell)] -= 0.5*ny*pVal;
+ vel[mesh.template getFaceNextToCell<0,+1>(cell)] += 0.5*ny*pVal;
+ }
+
+ void project(const MeshType& mesh)
+ {
+ int nx = mesh.getDimensions().x(), ny=mesh.getDimensions().y();
+ for ( int i=0 ; i< nx; i++ ) for (int j=0 ; j< ny; j++ )
+ computeVelocityDivergence(mesh.getCellIndex(CoordinatesType( i, j )), vel, mesh, p_rhs);
+
+ _matSolver.setMatrix(poissonMat);
+ _matSolver.solve(p_rhs,p);
+
+ for (int x=1; x< nx-1; x++) for(int y=1; y < ny-1; y++)
+ updateVelocityByPressureCorrection(mesh.getCellIndex(CoordinatesType(x,y)),vel, mesh, p);
+ }
+ void diffuse(RealType dt, const MeshType& mesh)
+ {
+ vel0 = vel;
+ //vytvorit matici s prvky
+ int nx = mesh.getDimensions().x(), ny=mesh.getDimensions().y();
+ double a = dt*visc*nx*nx, b = dt*visc*ny*ny;
+ prepareMatrix<1,0>(diffuseMat, mesh, nx+1, ny, -a, 1+2*a+2*b);
+ prepareMatrix<0,1>(diffuseMat, mesh, nx, ny+1, -b, 1+2*a+2*b);
+ _matSolver.setMatrix(diffuseMat);
+ _matSolver.solve(vel0,vel);
+ }
+ double getCenterU(const MeshType& mesh, IndexType cell) //x,y based on cells
+ {
+ return 0.5*(vel0[mesh.template getFaceNextToCell<-1,0>(cell)] + vel0[mesh.template getFaceNextToCell<+1,0>(cell)] );
+ }
+ double getCenterV(const MeshType& mesh, IndexType cell) //x,y based on cells
+ {
+ return 0.5*(vel0[mesh.template getFaceNextToCell<0,-1>(cell)] + vel0[mesh.template getFaceNextToCell<0,+1>(cell)] );
+ }
+ double getCrossU(const MeshType& mesh, int x, int y) //x,y based (n+1)*(n+1)
+ {
+ const CoordinatesType cellCoords(x,y);
+ const CoordinatesType downCoords(x,y-1);
+ return 0.5*(vel0[mesh.template getFaceNextToCell<-1,0>(mesh.getCellIndex(cellCoords))]
+ +vel0[mesh.template getFaceNextToCell<-1,0>(mesh.getCellIndex(downCoords))]);
+ }
+ double getCrossV(const MeshType& mesh, int x, int y) //x,y based (n+1)*(n+1)
+ {
+ const CoordinatesType cellCoords(x,y);
+ const CoordinatesType leftCoords(x-1,y);
+ return 0.5*(vel0[mesh.template getFaceNextToCell<0,-1>(mesh.getCellIndex(cellCoords))]
+ +vel0[mesh.template getFaceNextToCell<0,-1>(mesh.getCellIndex(leftCoords))]);
+ }
+ void advect(RealType dt, const MeshType& mesh)
+ {
+ vel0 = vel;
+ int nx = mesh.getDimensions().x(), ny=mesh.getDimensions().y();
+ //U has dimensions (nx+1,ny)
+ RealType cx = dt*nx, cy = dt*ny; //dt/h
+ for ( int i=1 ; i< nx ; i++ ) for (int j=1 ; j< ny-1 ; j++ )
+ {
+ const CoordinatesType cellCoordinates( i, j );
+ const CoordinatesType leftCoords(i-1,j);
+ IndexType cell = mesh.getCellIndex(cellCoordinates);
+ IndexType face = mesh.template getFaceNextToCell<-1,0>(cell);
+ vel[face] = vel0[face]
+ - cx*( square(getCenterU(mesh,cell)) - square(getCenterU(mesh, mesh.getCellIndex(leftCoords))) //(u^2)_x
+ + getCrossU(mesh,i,j+1)*getCrossV(mesh,i,j+1) - getCrossU(mesh,i,j)*getCrossV(mesh,i,j)
+ );
+ }
+ for ( int i=1 ; i< nx-1 ; i++ ) for (int j=1 ; j< ny ; j++ )
+ {
+ const CoordinatesType cellCoordinates( i, j );
+ const CoordinatesType downCoords(i,j-1);
+ IndexType cell = mesh.getCellIndex(cellCoordinates);
+ IndexType face = mesh.template getFaceNextToCell<0,-1>(cell);
+ vel[face] = vel0[face] - cy*(square(getCenterV(mesh,cell)) - square(getCenterV(mesh,mesh.getCellIndex(downCoords))) //(v^2)_y
+ + getCrossU(mesh,i+1,j)*getCrossV(mesh,i+1,j) - getCrossU(mesh,i,j)*getCrossV(mesh,i,j)
+ );
+ }
+ /*for (int j=0 ; j< ny ; j++ ) for ( int i=0 ; i < nx ; i++ )
+ {
+ const CoordinatesType cellCoords(i,j);
+ IndexType cell = mesh.getCellIndex(cellCoords);
+ int ii = mesh.template getFaceNextToCell<-1,0>(cell);
+ cout << ii << " ++ "<< i << " " << j << " "<< vel[ii] << endl;
+ }
+ cout << "----------" <<endl;
+ for (int j=0 ; j< ny ; j++ ) for ( int i=0 ; i < nx ; i++ )
+ {
+ const CoordinatesType cellCoords(i,j);
+ IndexType cell = mesh.getCellIndex(cellCoords);
+ int ii = mesh.template getFaceNextToCell<0,-1>(cell);
+ cout << ii << " ++ "<< i << " " << j << " "<< vel[ii] << endl;
+ }
+ int test = 3;
+ test++;*/
+ }
+ void save(const char * filename, const MeshType& mesh)
+ {
+ //FILE * pFile = fopen (filename, "w");
+ //fprintf(pFile, "#X Y u v\n");
+ int nx = mesh.getDimensions().x(), ny=mesh.getDimensions().y(), n=nx*ny;
+ int dims[] = {nx,ny,1};
+ double *vars = new double[n*3];
+ double *vvars[] = {vars};
+
+ int varDim[] = {3};
+ int centering[] = {0};
+ const char * names[] = {"Rychlost"};
+
+ for (IndexType j=0 ; j< ny ; j++ ) for ( IndexType i=0 ; i< nx ; i++ )
+ {
+ IndexType cell = mesh.getCellIndex(typename MeshType::CoordinatesType(i,j));
+ int ii = 3*(j*nx+i);
+ vars[ii+0] = getCenterU(mesh, cell);
+ vars[ii+1] = getCenterV(mesh, cell);
+ vars[ii+2] = 0;
+ //fprintf(pFile, "%lg %lg %lg %lg\n", (RealType)i, (RealType)j, getCenterU(mesh, cell), getCenterV(mesh, cell));
+ }
+ //fclose (pFile);
+ void write_regular_mesh(const char *filename, int useBinary, int *dims,
+ int nvars, int *vardim, int *centering,
+ const char * const *varnames, double **vars);
+ write_regular_mesh(filename, 0, dims, 1, varDim, centering, names, vvars );
+ delete[] vars;
+ }
+
+ template< int nx, int ny >
+ void prepareMatrix(MatrixType & matrix, const MeshType& mesh, IndexType Nx, IndexType Ny, RealType off, RealType diag)
+ {
+ for (IndexType y = 0; y < Ny; y++) for (IndexType x = 0; x < Nx; x++)
+ {
+ IndexType i = mesh.template getFaceIndex<nx, ny>(typename MeshType::CoordinatesType(x,y));
+ if (x==0 || x ==Nx-1 || y==0 || y==Ny-1) matrix.setElement(i,i,1);
+ else
+ {
+ matrix.setElement(i,i,diag);
+ matrix.setElement(i,mesh.template getFaceIndex<nx, ny>(typename MeshType::CoordinatesType(x-1,y)),off);
+ matrix.setElement(i,mesh.template getFaceIndex<nx, ny>(typename MeshType::CoordinatesType(x+1,y)),off);
+ matrix.setElement(i,mesh.template getFaceIndex<nx, ny>(typename MeshType::CoordinatesType(x,y-1)),off);
+ matrix.setElement(i,mesh.template getFaceIndex<nx, ny>(typename MeshType::CoordinatesType(x,y+1)),off);
+ }
+ }
+ }
protected:
- tnlSharedVector< RealType, DeviceType, IndexType > p;
-
- tnlStaticArray< Dimensions, tnlSharedVector< RealType, DeviceType, IndexType > > u;
+ RealType visc;
+ MatrixType poissonMat, diffuseMat;
+ LinearSolver _matSolver;
+
+ tnlVector<RealType, DeviceType, IndexType> vel, vel0, p, p_rhs;
+ //tnlStaticArray< Dimensions, tnlSharedVector< RealType, DeviceType, IndexType > > u, v;
DifferentialOperator differentialOperator;
@@ -109,3 +287,6 @@ class tnlIncompressibleNavierStokesProblem : public tnlPDEProblem< Mesh,
#include "tnlIncompressibleNavierStokesProblem_impl.h"
#endif /* TNLINCOMPRESSIBLENAVIERSTOKESPROBLEM_H_ */
+
+
+//Refaktor, do objektu, setup na parametry, laplace podle tnlLinearDiffusion
diff --git a/examples/incompressible-navier-stokes/tnlIncompressibleNavierStokesProblem_impl.h b/examples/incompressible-navier-stokes/tnlIncompressibleNavierStokesProblem_impl.h
index 259ea842d7..646d08201e 100644
--- a/examples/incompressible-navier-stokes/tnlIncompressibleNavierStokesProblem_impl.h
+++ b/examples/incompressible-navier-stokes/tnlIncompressibleNavierStokesProblem_impl.h
@@ -35,10 +35,11 @@ bool
tnlIncompressibleNavierStokesProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
setup( const tnlParameterContainer& parameters )
{
+ visc = parameters.getParameter< RealType >( "viscosity" );
/*if( ! this->boundaryCondition.setup( parameters, "boundary-conditions-" ) ||
! this->rightHandSide.setup( parameters, "right-hand-side-" ) )
- return false;
- return true;*/
+ return false;*/
+ return true;
}
template< typename Mesh,
@@ -76,9 +77,9 @@ tnlIncompressibleNavierStokesProblem< Mesh, BoundaryCondition, RightHandSide, Di
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 );
+ //const IndexType uDofs = mesh.getNumberOfFaces() / 2;
+ /*this->u[ 0 ].bind( dofVector.getData(), uDofs );
+ this->u[ 1 ].bind( &dofVector.getData()[ uDofs ], uDofs );*/
}
/*template< typename Mesh,
@@ -107,9 +108,46 @@ setInitialCondition( const tnlParameterContainer& parameters,
{
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 );
+
+ vel.setSize(mesh.getNumberOfFaces());
+ vel0.setSize(mesh.getNumberOfFaces());
+ p.setSize(mesh.getNumberOfCells());
+ p_rhs.setSize(mesh.getNumberOfCells());
+
+ vel.setValue(0); vel0.setValue(0);
+ p.setValue(0); p_rhs.setValue(0);
+
+ RealType initVel = parameters.getParameter< RealType >( "inletVelocity" );
+ int nx = mesh.getDimensions().x(), ny=mesh.getDimensions().y();
+ typename MatrixType::RowLengthsVector rowLenghts;
+ IndexType N = mesh.getNumberOfCells();
+ IndexType Nf = mesh.getNumberOfFaces();
+ rowLenghts.setSize(N);
+ rowLenghts.setValue(5);
+ poissonMat.setDimensions( N, N ); diffuseMat.setDimensions( Nf, Nf );
+ poissonMat.setRowLengths(rowLenghts);
+ rowLenghts.setSize(Nf);
+ rowLenghts.setValue(5);
+ diffuseMat.setRowLengths(rowLenghts);
+ for ( int x=0 ; x<nx ; x++ ) for (int y=0 ; y<ny ; y++ )
+ {
+ int i = mesh.getCellIndex(typename MeshType::CoordinatesType(x,y));
+ if (y==ny-1 && x > 0 && x < nx)
+ {
+ float xi = (x-1.0)/(nx-1.0);
+ vel[mesh.template getFaceNextToCell<-1,0>(i)] = (1-square((xi-0.5)*2))*initVel; //inital conditon for vel
+ }
+ if (x==0 || y==0 || x==nx-1 || y==ny-1)
+ {
+ poissonMat.setElement(i,i,1);
+ }else{
+ poissonMat.setElement(i,i,4);
+ poissonMat.setElement(i,mesh.getCellIndex(typename MeshType::CoordinatesType(x-1,y)),-1);
+ poissonMat.setElement(i,mesh.getCellIndex(typename MeshType::CoordinatesType(x+1,y)),-1);
+ poissonMat.setElement(i,mesh.getCellIndex(typename MeshType::CoordinatesType(x,y-1)),-1);
+ poissonMat.setElement(i,mesh.getCellIndex(typename MeshType::CoordinatesType(x,y+1)),-1);
+ }
+ }
/*const tnlString& initialConditionFile = parameters.getParameter< tnlString >( "initial-condition" );
if( ! this->solution.load( initialConditionFile ) )
@@ -124,7 +162,7 @@ template< typename Mesh,
typename BoundaryCondition,
typename RightHandSide,
typename DifferentialOperator >
- template< typename MatrixType >
+ // template< typename MatrixType >
bool
tnlIncompressibleNavierStokesProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
setupLinearSystem( const MeshType& mesh,
@@ -165,7 +203,8 @@ makeSnapshot( const RealType& time,
this->bindAuxiliaryDofs( mesh, auxiliaryDofs );
//cout << "dofs = " << dofs << endl;
tnlString fileName;
- FileNameBaseNumberEnding( "u-", step, 5, ".tnl", fileName );
+ FileNameBaseNumberEnding( "u-", step, 5, ".vtk", fileName );
+ save("test.txt", mesh);
//if( ! this->solution.save( fileName ) )
// return false;
return true;
@@ -214,7 +253,7 @@ template< typename Mesh,
typename BoundaryCondition,
typename RightHandSide,
typename DifferentialOperator >
- template< typename MatrixType >
+ // template< typename MatrixType >
void
tnlIncompressibleNavierStokesProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
assemblyLinearSystem( const RealType& time,
@@ -241,4 +280,5 @@ assemblyLinearSystem( const RealType& time,
abort();*/
}
+
#endif /* TNLINCOMPRESSIBLENAVIERSTOKESPROBLEM_IMPL_H_ */
diff --git a/examples/incompressible-navier-stokes/tnlLinearDiffusion.h b/examples/incompressible-navier-stokes/tnlLinearDiffusion.h
deleted file mode 100644
index 35eebca181..0000000000
--- a/examples/incompressible-navier-stokes/tnlLinearDiffusion.h
+++ /dev/null
@@ -1,168 +0,0 @@
-#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
deleted file mode 100644
index 4838978478..0000000000
--- a/examples/incompressible-navier-stokes/tnlLinearDiffusion_impl.h
+++ /dev/null
@@ -1,272 +0,0 @@
-
-#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/install b/install
index 0c15ff9b1e..731ccf9aa5 100755
--- a/install
+++ b/install
@@ -29,15 +29,15 @@ ${CMAKE} .. -DCMAKE_BUILD_TYPE=Debug \
-DPETSC_DIR=${PETSC_DIR} \
-DWITH_TEMPLATE_EXPLICIT_INSTANTIATION=${TEMPLATE_EXPLICIT_INSTANTIATION}
make -j${CPUS} ${VERBOSE}
-make -j${CPUS} test
-make -j${CPUS} install
-
-cd ../Release
-${CMAKE} .. -DCMAKE_INSTALL_PREFIX=${HOME}/local \
- -DWITH_CUDA=${WITH_CUDA} \
- -DPETSC_DIR=${PETSC_DIR} \
- -DWITH_TEMPLATE_EXPLICIT_INSTANTIATION=${TEMPLATE_EXPLICIT_INSTANTIATION}
-make -j${CPUS} ${VERBOSE}
-make -j${CPUS} test
-make -j${CPUS} install
+#make -j${CPUS} test
+#make -j${CPUS} install
+
+#cd ../Release
+#${CMAKE} .. -DCMAKE_INSTALL_PREFIX=${HOME}/local \
+# -DWITH_CUDA=${WITH_CUDA} \
+# -DPETSC_DIR=${PETSC_DIR} \
+# -DWITH_TEMPLATE_EXPLICIT_INSTANTIATION=${TEMPLATE_EXPLICIT_INSTANTIATION}
+#make -j${CPUS} ${VERBOSE}
+#make -j${CPUS} test
+#make -j${CPUS} install
diff --git a/src/solvers/linear/tnlLinearResidueGetter_impl.h b/src/solvers/linear/tnlLinearResidueGetter_impl.h
index b41c99a027..daf106689c 100644
--- a/src/solvers/linear/tnlLinearResidueGetter_impl.h
+++ b/src/solvers/linear/tnlLinearResidueGetter_impl.h
@@ -34,7 +34,8 @@ typename tnlLinearResidueGetter< Matrix, Vector > :: RealType
RealType err = fabs( matrix. rowVectorProduct( i, x ) - b[ i ] );
res += err * err;
}
- return sqrt( res ) / bNorm;
+ if (bNorm!=0.0) return sqrt( res ) / bNorm;
+ return sqrt(res);
}
#endif /* TNLLINEARRESIDUEGETTER_IMPL_H_ */
--
GitLab