diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index fb837c7015557471266808ab0ec262696cef4527..68dc872c9365da1eb66dbc2778cd16f4b30dada3 100755
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1,2 +1,3 @@
add_subdirectory( make-project )
-add_subdirectory( simple-solver )
\ No newline at end of file
+add_subdirectory( simple-solver )
+add_subdirectory( navier-stokes )
\ No newline at end of file
diff --git a/examples/navier-stokes/Makefile b/examples/navier-stokes/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..aa5e5bf75c496d54089aa82ac1f7fb14a5e94d2d
--- /dev/null
+++ b/examples/navier-stokes/Makefile
@@ -0,0 +1,60 @@
+TNL_VERSION=0.1
+TNL_INSTALL_DIR=${HOME}/local/lib
+TNL_INCLUDE_DIR=${HOME}/local/include/tnl-${TNL_VERSION}
+#TNL_INCLUDE_DIR=${HOME}/workspace/tnl/src
+
+TARGET = navier-stokes
+CONFIG_FILE = $(TARGET).cfg.desc
+INSTALL_DIR = ${HOME}/local
+CXX = g++
+CUDA_CXX = nvcc
+OMP_FLAGS = -DHAVE_OPENMP -fopenmp
+#CXX_FLAGS = -std=gnu++0x -I$(TNL_INCLUDE_DIR) -O0 -g -DDEBUG $(OMP_FLAGS)
+#CXX_FLAGS = -std=gnu++0x -I$(TNL_INCLUDE_DIR) -O3 $(OMP_FLAGS)
+CXX_FLAGS = -DHAVE_NOT_CXX11 -I$(TNL_INCLUDE_DIR) -O3 $(OMP_FLAGS)
+LD_FLAGS = -L$(TNL_INSTALL_DIR) -ltnl-0.1 -lgomp
+
+SOURCES = main.cpp
+HEADERS = navierStokesSetter.h \
+ navierStokesSetter_impl.h \
+ navierStokesSolver.h \
+ navierStokesSolver_impl.h \
+ navierStokesSolverMonitor.h \
+ navierStokesSolverMonitor_impl.h
+OBJECTS = main.o
+DIST = $(SOURCES) $(HEADERS) $(CONFIG_FILE) Makefile
+
+all: $(TARGET)
+
+clean:
+ rm -f $(OBJECTS)
+ rm -f $(TARGET)-conf.h
+
+dist: $(DIST)
+ tar zcvf $(TARGET).tgz $(DIST)
+
+install: $(TARGET)
+ cp $(TARGET) $(INSTALL_DIR)/bin
+ cp $(CONFIG_FILE) $(INSTALL_DIR)/share
+ cp -rv share/* $(INSTALL_DIR)/share/$(TARGET)
+ cp make-png-from-gnuplot $(INSTALL_DIR)/bin
+ cp merge-figures $(INSTALL_DIR)/bin
+ chmod +x $(INSTALL_DIR)/bin/make-png-from-gnuplot
+ chmod +x $(INSTALL_DIR)/bin/merge-figures
+
+uninstall: $(TARGET)
+ rm -f $(INSTALL_DIR)/bin/$(TARGET)
+ rm -f $(CONFIG_FILE) $(INSTALL_DIR)/share
+ rm -rf $(INSTALL_DIR)/share/$(TARGET)
+ rm -rf $(INSTALL_DIR)/bin/make-png-from-gnuplot
+ rm -rf $(INSTALL_DIR)/bin/merge-figures
+
+$(TARGET): $(OBJECTS)
+ $(CXX) -o $(TARGET) $(OBJECTS) $(LD_FLAGS)
+
+%.o: %.cpp $(TARGET)-conf.h $(HEADERS)
+ $(CXX) -c -o $@ $(CXX_FLAGS) $<
+
+$(TARGET)-conf.h:
+ echo "#define CONFIG_FILE \"${INSTALL_DIR}/share/${CONFIG_FILE}\" " > $(TARGET)-conf.h
+
diff --git a/examples/navier-stokes/laxFridrichs.h b/examples/navier-stokes/laxFridrichs.h
new file mode 100644
index 0000000000000000000000000000000000000000..275db9813dc8170122f40c196816718c88a248c9
--- /dev/null
+++ b/examples/navier-stokes/laxFridrichs.h
@@ -0,0 +1,92 @@
+/***************************************************************************
+ laxFridrichs.h - description
+ -------------------
+ begin : Mar 1, 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 LAXFRIDRICHS_H_
+#define LAXFRIDRICHS_H_
+
+template< typename MeshType >
+class laxFridrichs
+{
+ public:
+
+ typedef typename MeshType :: RealType RealType;
+ typedef typename MeshType :: DeviceType DeviceType;
+ typedef typename MeshType :: IndexType IndexType;
+
+ template< typename ConservativeVector,
+ typename VelocityVector >
+// typename PressureVector >
+ void getExplicitRhs( const MeshType& mesh,
+ const IndexType centralVolume,
+ const ConservativeVector& rho,
+ const ConservativeVector& rho_u1,
+ const ConservativeVector& rho_u2,
+ const VelocityVector& u1,
+ const VelocityVector& u2,
+// const PressureVector& p,
+ ConservativeVector& rho_t,
+ ConservativeVector& rho_u1_t,
+ ConservativeVector& rho_u2_t,
+ const typename MeshType :: RealType viscosityCoeff = 1.0 ) const;
+};
+
+template< typename MeshType >
+ template< typename ConservativeVector,
+ typename VelocityVector >
+// typename PressureVector >
+void laxFridrichs< MeshType > :: getExplicitRhs( const MeshType& mesh,
+ const IndexType centralVolume,
+ const ConservativeVector& rho,
+ const ConservativeVector& rho_u1,
+ const ConservativeVector& rho_u2,
+ const VelocityVector& u1,
+ const VelocityVector& u2,
+ ConservativeVector& rho_t,
+ ConservativeVector& rho_u1_t,
+ ConservativeVector& rho_u2_t,
+ const typename MeshType :: RealType viscosityCoeff ) const
+{
+ const IndexType& xSize = mesh. getDimensions(). x();
+ const IndexType& ySize = mesh. getDimensions(). y();
+ const RealType hx = mesh. getSpaceStep(). x();
+ const RealType hy = mesh. getSpaceStep(). y();
+
+ const IndexType& c = centralVolume;
+ const IndexType e = mesh. getNodeNeighbour( centralVolume, 0, 1 );
+ const IndexType w = mesh. getNodeNeighbour( centralVolume, 0, -1 );
+ const IndexType n = mesh. getNodeNeighbour( centralVolume, 1, 0 );
+ const IndexType s = mesh. getNodeNeighbour( centralVolume, -1, 0 );
+
+ /****
+ * rho_t + ( rho u_1 )_x + ( rho u_2 )_y = 0
+ */
+ rho_t[ c ]= viscosityCoeff * 0.25 * ( rho[ e ] + rho[ w ] + rho[ s ] + rho[ n ] - 4.0 * rho[ c ] )
+ - ( rho[ e ] * u1[ e ] - rho[ w ] * u1[ w ] ) / ( 2.0 * hx )
+ - ( rho[ n ] * u2[ n ] - rho[ s ] * u2[ s ] ) / ( 2.0 * hy );
+
+ /****
+ * ( rho * u1 )_t + ( rho * u1 * u1 )_x + ( rho * u1 * u2 )_y = 0
+ */
+ rho_u1_t[ c ] = viscosityCoeff * 0.25 * ( rho_u1[ e ] + rho_u1[ w ] + rho_u1[ s ] + rho_u1[ n ] - 4.0 * rho_u1[ c ] )
+ - ( rho_u1[ e ] * u1[ e ] - rho_u1[ w ] * u1[ w ] ) / ( 2.0 * hx )
+ - ( rho_u1[ n ] * u2[ n ] - rho_u1[ s ] * u2[ s ] ) / ( 2.0 * hy );
+ rho_u2_t[ c ] = viscosityCoeff * 0.25 * ( rho_u2[ e ] + rho_u2[ w ] + rho_u2[ s ] + rho_u2[ n ] - 4.0 * rho_u2[ c ] )
+ - ( rho_u2[ e ] * u1[ e ] - rho_u2[ w ] * u1[ w ] ) / ( 2.0 * hx )
+ - ( rho_u2[ n ] * u2[ n ] - rho_u2[ s ] * u2[ s ] ) / ( 2.0 * hy );
+}
+
+#endif
diff --git a/examples/navier-stokes/main.cpp b/examples/navier-stokes/main.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..78d27e70fea8bd40c8a39b0256039786e0b4f073
--- /dev/null
+++ b/examples/navier-stokes/main.cpp
@@ -0,0 +1,31 @@
+/***************************************************************************
+ main.cpp - description
+ -------------------
+ begin : Jan 12, 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. *
+ * *
+ ***************************************************************************/
+
+#include <cstdlib>
+#include "navier-stokes-conf.h"
+#include "navierStokesSetter.h"
+#include <solvers/tnlSolver.h>
+
+int main( int argc, char* argv[] )
+{
+ tnlSolver< navierStokesSetter > solver;
+ if( ! solver. run( CONFIG_FILE, argc, argv ) )
+ return EXIT_FAILURE;
+ return EXIT_SUCCESS;
+}
+
+
diff --git a/examples/navier-stokes/make-png-from-gnuplot b/examples/navier-stokes/make-png-from-gnuplot
new file mode 100644
index 0000000000000000000000000000000000000000..1fb64ae9de3c4a572be6810a42e442946e93d532
--- /dev/null
+++ b/examples/navier-stokes/make-png-from-gnuplot
@@ -0,0 +1,38 @@
+#!/bin/bash
+
+function processFile()
+{
+ file=${1}
+
+ gnuplotcommand="
+ set terminal png giant size 1280,1280 crop;
+ set output '`basename $file ".gplt"`.png';
+ set pm3d map;
+ set palette defined(0.0 0.5 1.0 0, 0.02 \"light-goldenrod\", 0.04 \"yellow\", 0.08 \"red\", 0.4 \"light-blue\", 1.0 \"blue\");
+ unset key;
+ set size ratio -1;
+ set pointsize 0.4;"
+
+ if ! test x$2 = x;
+ then
+ gnuplotcommand="${gnuplotcommand} set xrange [0:$2];"
+ fi
+ if ! test x$3 = x;
+ then
+ gnuplotcommand="${gnuplotcommand} set yrange [0:$3];"
+ fi
+ if ! test x$4 = x;
+ then
+ gnuplotcommand="${gnuplotcommand} set cbrange [0:$4];"
+ fi
+
+ gnuplotcommand="${gnuplotcommand} splot '$file' w pm3d title '${1}';"
+
+ echo ${gnuplotcommand} | gnuplot
+}
+
+for file in ${1}*.gplt
+do
+ echo -ne "Creating: `basename $file ".gplt"`.png \r"
+ processFile ${file} ${2} ${3} ${4}
+done
diff --git a/examples/navier-stokes/merge-figures b/examples/navier-stokes/merge-figures
new file mode 100755
index 0000000000000000000000000000000000000000..2cd31edde2f9b21ffdef20ae5c9ef05027aec157
--- /dev/null
+++ b/examples/navier-stokes/merge-figures
@@ -0,0 +1,10 @@
+#!/bin/bash
+
+for file in ${1}*.png;
+do
+ numbers=${file//[^0-9]/}
+ echo -ne "montage $file ${2}${numbers}.png ${3}${numbers}.png -tile 3x1 -geometry +10+10 ${4}${numbers}.png \r"
+ montage $file ${2}${numbers}.png ${3}${numbers}.png -tile 3x1 -geometry +10+10 ${4}${numbers}.png
+
+
+done
diff --git a/examples/navier-stokes/navier-stokes.cfg.desc b/examples/navier-stokes/navier-stokes.cfg.desc
new file mode 100644
index 0000000000000000000000000000000000000000..3b0d8920c55dc51ccbc3eef478f5e6eaa69db8a2
--- /dev/null
+++ b/examples/navier-stokes/navier-stokes.cfg.desc
@@ -0,0 +1,45 @@
+group IO
+{
+ string input-file [Input file name.];
+ string output-file [Output file name.];
+ string log-file("twoface-flow.log") [File name for the log.];
+},[Arguments describing input and output data.];
+
+group Problem
+{
+ integer dimensions(2) [Number of the problem dimensions. Can be only 2 for now.];
+ string problem-name(!) [Problem to solve. It can be riser, cavity,...];
+ integer x-size(!) [Number of grid cells along x-axis.];
+ integer y-size(!) [Number of grid cells along y-axis.];
+ real max-inflow-velocity(1.0) [];
+ real max-outflow-velocity(0.8) [];
+ real start-up(1.0) [Start-up time of the riser.];
+ real mu(!) [Viscosity.];
+ real p0(!) [Pressure for the initial condition.];
+ real T(!) [Temperature.];
+ real R(!) [Gas constant.];
+ real gravity(9.81) [Gravity force.];
+},[Setting up the problem we solve.];
+
+group Geometry
+{
+ real width(!) [Width of the domain in meters (x-axis).];
+ real height(!) [Height of the domain in meters (y-axis).];
+ real riser-gas-inflow-width[0.2] [Width of the gas inflow on the bottom of the riser.];
+},[Parameters describing the geometry of the problem.];
+
+group Solver
+{
+ string real-type( "double" ) [Precision of the floating point arithmetics. Can be float, double, long-double.];
+ string index-type( "int" ) [Indexing type for arrays, vectors, matrices etc. Can be int or long-int.];
+ string device( "host" ) [Device to use for the computations. Can be host or cuda.];
+ string discrete-solver(!) [The main solver for the discretised problem.];
+ real merson-adaptivity(1.0e-4) [Parameter controling adaptivity of the Runge-Kutta-Merson solver.];
+ real initial-tau(1.0e-5) [Initial time step for the time dependent PDE solver. It can be changed if an adaptive solver is used.];
+ real snapshot-period(!) [Time period for writing a state of the time dependent problem to a file.];
+ real final-time(!) [Stop time of the time dependent simulation.];
+ string scheme(!) [Numerical scheme for the space discretisation.];
+ string time-discretisation("explicit") [Time discretisation of the problem. Can be only explicit for now.];
+ integer verbose(1) [Set the verbose mode. The higher number the more messages are generated.];
+},[Parameters of the solver];
+
diff --git a/examples/navier-stokes/navierStokesSetter.h b/examples/navier-stokes/navierStokesSetter.h
new file mode 100644
index 0000000000000000000000000000000000000000..6d946903a2c0760f176b3023f303ca282f4f4160
--- /dev/null
+++ b/examples/navier-stokes/navierStokesSetter.h
@@ -0,0 +1,38 @@
+/***************************************************************************
+ simpleProblemSetter.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 SIMPLEPROBLEMTYPESSETTER_H_
+#define SIMPLEPROBLEMTYPESSETTER_H_
+
+#include <config/tnlParameterContainer.h>
+#include <mesh/tnlGrid.h>
+#include "navierStokesSolver.h"
+
+template< typename SolverStarter >
+class navierStokesSetter
+{
+ public:
+ template< typename RealType,
+ typename DeviceType,
+ typename IndexType >
+ bool run( const tnlParameterContainer& parameters ) const;
+};
+
+#include "navierStokesSetter_impl.h"
+
+
+#endif /* SIMPLEPROBLEMSETTER_H_ */
diff --git a/examples/navier-stokes/navierStokesSetter_impl.h b/examples/navier-stokes/navierStokesSetter_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..bae5bb8434386384866bc5010dae83125547811b
--- /dev/null
+++ b/examples/navier-stokes/navierStokesSetter_impl.h
@@ -0,0 +1,42 @@
+/***************************************************************************
+ navierStokesSetter_impl.h - description
+ -------------------
+ begin : Mar 9, 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 NAVIERSTOKESSETTER_IMPL_H_
+#define NAVIERSTOKESSETTER_IMPL_H_
+
+template< typename SolverStarter >
+ template< typename RealType,
+ typename DeviceType,
+ typename IndexType >
+bool navierStokesSetter< SolverStarter > :: run( const tnlParameterContainer& parameters ) const
+{
+ int dimensions = parameters. GetParameter< int >( "dimensions" );
+ if( dimensions != 2 )
+ {
+ cerr << "The problem is not defined for " << dimensions << "dimensions." << endl;
+ return false;
+ }
+ SolverStarter solverStarter;
+ if( dimensions == 2 )
+ {
+ typedef tnlGrid< 2, RealType, DeviceType, IndexType > MeshType;
+ return solverStarter. run< navierStokesSolver< MeshType > >( parameters );
+ }
+}
+
+
+#endif /* NAVIERSTOKESSETTER_IMPL_H_ */
diff --git a/examples/navier-stokes/navierStokesSolver.h b/examples/navier-stokes/navierStokesSolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..b87fab0a03fe5922a880f162833635317abf8f73
--- /dev/null
+++ b/examples/navier-stokes/navierStokesSolver.h
@@ -0,0 +1,107 @@
+/***************************************************************************
+ navierStokesSolver.h - description
+ -------------------
+ begin : Jan 13, 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 NAVIERSTOKESSOLVER_H_
+#define NAVIERSTOKESSOLVER_H_
+
+#include <core/tnlLogger.h>
+#include <core/tnlHost.h>
+#include <core/tnlVector.h>
+#include <config/tnlParameterContainer.h>
+#include <matrix/tnlCSRMatrix.h>
+#include <solvers/preconditioners/tnlDummyPreconditioner.h>
+#include <solvers/tnlSolverMonitor.h>
+#include "navierStokesSolverMonitor.h"
+#include "laxFridrichs.h"
+
+template< typename Mesh >
+class navierStokesSolver
+{
+ public:
+
+ typedef typename Mesh :: RealType RealType;
+ typedef typename Mesh :: DeviceType DeviceType;
+ typedef typename Mesh :: IndexType IndexType;
+ typedef Mesh MeshType;
+ typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
+
+ typedef tnlCSRMatrix< RealType, DeviceType, IndexType > DiscreteSolverMatrixType;
+ typedef tnlDummyPreconditioner< RealType, DeviceType, IndexType > DiscreteSolverPreconditioner;
+
+ enum BoundaryConditionType { dirichlet, neumann, noSlip };
+
+ enum CfdScheme { laxFridrichsEnm, upwindEnm };
+
+ enum ProblemType { riser, cavity };
+
+ navierStokesSolver();
+
+ static tnlString getTypeStatic();
+
+ tnlString getPrologHeader() const;
+
+ void writeProlog( tnlLogger& logger,
+ const tnlParameterContainer& parameters ) const;
+
+ bool init( const tnlParameterContainer& parameters );
+
+ bool setInitialCondition( const tnlParameterContainer& parameters );
+
+ DofVectorType& getDofVector();
+
+ bool makeSnapshot( const RealType& t,
+ const IndexType step );
+
+ bool solve();
+
+ void GetExplicitRHS( const RealType& time,
+ const RealType& tau,
+ DofVectorType& _u,
+ DofVectorType& _fu );
+
+ tnlSolverMonitor< RealType, IndexType >* getSolverMonitor();
+
+ protected:
+
+ RealType regularize( const RealType& r ) const;
+
+ template< typename Vector >
+ void updatePhysicalQuantities( const Vector& rho,
+ const Vector& rho_u1,
+ const Vector& rho_u2 );
+
+ ProblemType problem;
+
+ CfdScheme scheme;
+
+ MeshType mesh;
+
+ tnlVector< RealType, DeviceType, IndexType > rho, u1, u2, p;
+
+ DofVectorType dofVector, rhsDofVector;
+
+ RealType mu, R, T, p_0, gravity,
+ maxInflowVelocity, maxOutflowVelocity, startUp;
+
+ laxFridrichs< MeshType > laxFridrichsScheme;
+
+ navierStokesSolverMonitor< RealType, IndexType > solverMonitor;
+};
+
+#include "navierStokesSolver_impl.h"
+
+#endif /* NAVIERSTOKESSOLVER_H_ */
diff --git a/examples/navier-stokes/navierStokesSolverMonitor.h b/examples/navier-stokes/navierStokesSolverMonitor.h
new file mode 100644
index 0000000000000000000000000000000000000000..23c1a129344bc18575fb8e2ea89089cf3a51fbf1
--- /dev/null
+++ b/examples/navier-stokes/navierStokesSolverMonitor.h
@@ -0,0 +1,69 @@
+/***************************************************************************
+ navierStokesSolverMonitor.h - description
+ -------------------
+ begin : Mar 13, 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 NAVIERSTOKESSOLVERMONITOR_H_
+#define NAVIERSTOKESSOLVERMONITOR_H_
+
+#include <solvers/ode/tnlODESolverMonitor.h>
+
+template< typename Real, typename Index >
+class navierStokesSolverMonitor : public tnlODESolverMonitor< Real, Index >
+{
+ public:
+
+ navierStokesSolverMonitor();
+
+ void refresh();
+
+ void setBetaMax( const Real& betaMax );
+
+ void setEpsSMax( const Real& epsSMax );
+
+ void setUSMax( const Real& u_s_max );
+
+ void setUGMax( const Real& u_g_max );
+
+ void setEps_rho_g_max( const Real& eps_rho_g_max );
+ void setEps_rho_u1_g_max( const Real& eps_rho_u1_g_max );
+ void setEps_rho_u2_g_max( const Real& eps_rho_u2_g_max );
+ void setEps_rho_g_t_max( const Real& eps_rho_g_t_max );
+ void setEps_rho_u1_g_t_max( const Real& eps_rho_u1_g_t_max );
+ void setEps_rho_u2_g_t_max( const Real& eps_rho_u2_g_t_max );
+ void setEps_rho_s_max( const Real& eps_rho_s_max );
+ void setEps_rho_u1_s_max( const Real& eps_rho_u1_s_max );
+ void setEps_rho_u2_s_max( const Real& eps_rho_u2_s_max );
+ void setEps_rho_s_t_max( const Real& eps_rho_s_t_max );
+ void setEps_rho_u1_s_t_max( const Real& eps_rho_u1_s_t_max );
+ void setEps_rho_u2_s_t_max( const Real& eps_rho_u2_s_t_max );
+
+ protected:
+
+ public:
+
+ Real beta_max, eps_s_max, u_s_max, u_g_max, G_max;
+ Real eps_rho_g_max, eps_rho_u1_g_max, eps_rho_u2_g_max,
+ eps_rho_g_t_max, eps_rho_u1_g_t_max, eps_rho_u2_g_t_max,
+ eps_rho_s_max, eps_rho_u1_s_max, eps_rho_u2_s_max,
+ eps_rho_s_t_max, eps_rho_u1_s_t_max, eps_rho_u2_s_t_max;
+
+ public:
+ Index u_s_max_i, u_s_max_j;
+};
+
+#include "navierStokesSolverMonitor_impl.h"
+
+#endif /* NAVIERSTOKESSOLVERMONITOR_H_ */
diff --git a/examples/navier-stokes/navierStokesSolverMonitor_impl.h b/examples/navier-stokes/navierStokesSolverMonitor_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..bbca7a7a4d5c7d5f8c1ae8aa2df0fe604772f554
--- /dev/null
+++ b/examples/navier-stokes/navierStokesSolverMonitor_impl.h
@@ -0,0 +1,257 @@
+/***************************************************************************
+ navierStokesSolverMonitor_impl.h - description
+ -------------------
+ begin : Mar 13, 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 TNLNAVIERSTOKESSOLVERMONITOR_IMPL_H_
+#define TNLNAVIERSTOKESSOLVERMONITOR_IMPL_H_
+
+#include <fstream>
+
+using namespace std;
+
+template< typename Real, typename Index >
+navierStokesSolverMonitor< Real, Index > :: navierStokesSolverMonitor()
+: beta_max( 0.0 ),
+ eps_s_max( 0.0 ),
+ u_s_max( 0.0 ),
+ u_g_max( 0.0 ),
+ G_max( 0.0 ),
+ eps_rho_g_max( 0.0 ),
+ eps_rho_u1_g_max( 0.0 ),
+ eps_rho_u2_g_max( 0.0 ),
+ eps_rho_g_t_max( 0.0 ),
+ eps_rho_u1_g_t_max( 0.0 ),
+ eps_rho_u2_g_t_max( 0.0 ),
+ eps_rho_s_max( 0.0 ),
+ eps_rho_u1_s_max( 0.0 ),
+ eps_rho_u2_s_max( 0.0 ),
+ eps_rho_s_t_max( 0.0 ),
+ eps_rho_u1_s_t_max( 0.0 ),
+ eps_rho_u2_s_t_max( 0.0 )
+{
+ /****
+ * Reset statistics files
+ */
+ fstream file;
+ file. open( "tau.txt", ios :: out );
+ file. close();
+ file. open( "beta-max.txt", ios :: out );
+ file. close();
+ file. open( "eps-s-max.txt", ios :: out );
+ file. close();
+ file. open( "u-s-max.txt", ios :: out );
+ file. close();
+ file. open( "u-g-max.txt", ios :: out );
+ file. close();
+ file. open( "eps-rho-g-max.txt", ios :: out );
+ file. close();
+ file. open( "eps-rho-u1-g-max.txt", ios :: out );
+ file. close();
+ file. open( "eps-rho-u2-g-max.txt", ios :: out );
+ file. close();
+ file. open( "eps-rho-g-t-max.txt", ios :: out );
+ file. close();
+ file. open( "eps-rho-u1-g-t-max.txt", ios :: out );
+ file. close();
+ file. open( "eps-rho-u2-g-t-max.txt", ios :: out );
+ file. close();
+ file. open( "eps-rho-s-max.txt", ios :: out );
+ file. close();
+ file. open( "eps-rho-u1-s-max.txt", ios :: out );
+ file. close();
+ file. open( "eps-rho-u2-s-max.txt", ios :: out );
+ file. close();
+ file. open( "eps-rho-s-t-max.txt", ios :: out );
+ file. close();
+ file. open( "eps-rho-u1-s-t-max.txt", ios :: out );
+ file. close();
+ file. open( "eps-rho-u2-s-t-max.txt", ios :: out );
+ file. close();
+}
+
+template< typename Real, typename Index >
+void navierStokesSolverMonitor< Real, Index > :: refresh()
+{
+ /****
+ * Write to statistics files
+ */
+ fstream file;
+ file. open( "tau.txt", ios :: out | ios :: app );
+ file << this -> getTime() << " " << log10( this -> getTimeStep() ) << endl;
+ file. close();
+ file. open( "beta-max.txt", ios :: out | ios :: app );
+ file << this -> getTime() << " " << log10( this -> beta_max ) << endl;
+ file. close();
+ file. open( "eps-s-max.txt", ios :: out | ios :: app );
+ file << this -> getTime() << " " << this -> eps_s_max << endl;
+ file. close();
+ file. open( "u-s-max.txt", ios :: out | ios :: app );
+ file << this -> getTime() << " " << this -> u_s_max << endl;
+ file. close();
+ file. open( "u-g-max.txt", ios :: out | ios :: app );
+ file << this -> getTime() << " " << this -> u_g_max << endl;
+ file. close();
+ file. open( "eps-rho-g-max.txt", ios :: out | ios :: app );
+ file << this -> getTime() << " " << this -> eps_rho_g_max << endl;
+ file. close();
+ file. open( "eps-rho-u1-g-max.txt", ios :: out | ios :: app );
+ file << this -> getTime() << " " << this -> eps_rho_u1_g_max << endl;
+ file. close();
+ file. open( "eps-rho-u2-g-max.txt", ios :: out | ios :: app );
+ file << this -> getTime() << " " << this -> eps_rho_u2_g_max << endl;
+ file. close();
+ file. open( "eps-rho-g-t-max.txt", ios :: out | ios :: app );
+ file << this -> getTime() << " " << this -> eps_rho_g_t_max << endl;
+ file. close();
+ file. open( "eps-rho-u1-g-t-max.txt", ios :: out | ios :: app );
+ file << this -> getTime() << " " << this -> eps_rho_u1_g_t_max << endl;
+ file. close();
+ file. open( "eps-rho-u2-g-t-max.txt", ios :: out | ios :: app );
+ file << this -> getTime() << " " << this -> eps_rho_u2_g_t_max << endl;
+ file. close();
+ file. open( "eps-rho-s-max.txt", ios :: out | ios :: app );
+ file << this -> getTime() << " " << this -> eps_rho_s_max << endl;
+ file. close();
+ file. open( "eps-rho-u1-s-max.txt", ios :: out | ios :: app );
+ file << this -> getTime() << " " << this -> eps_rho_u1_s_max << endl;
+ file. close();
+ file. open( "eps-rho-u2-s-max.txt", ios :: out | ios :: app );
+ file << this -> getTime() << " " << this -> eps_rho_u2_s_max << endl;
+ file. close();
+ file. open( "eps-rho-s-t-max.txt", ios :: out | ios :: app );
+ file << this -> getTime() << " " << this -> eps_rho_s_t_max << endl;
+ file. close();
+ file. open( "eps-rho-u1-s-t-max.txt", ios :: out | ios :: app);
+ file << this -> getTime() << " " << this -> eps_rho_u1_s_t_max << endl;
+ file. close();
+ file. open( "eps-rho-u2-s-t-max.txt", ios :: out | ios :: app );
+ file << this -> getTime() << " " << this -> eps_rho_u2_s_t_max << endl;
+ file. close();
+
+
+
+ if( this -> verbose > 0 && this -> refreshing % this -> outputPeriod == 0 )
+ {
+ /*cout << "u-max I: "<< setw( 3 ) << this -> u_s_max_i
+ << " J: "<< setw( 3 ) << this -> u_s_max_j
+ << " G_max: " << setw( 6 ) << this -> G_max;*/
+ //cout << setprecision( 5 )
+ // << "betaMax: " << setw( 8 ) << this -> beta_max
+ // << " epsSMax: " << setw( 8 ) << this -> eps_s_max;
+ }
+ tnlODESolverMonitor< Real, Index > :: refresh();
+}
+
+template< typename Real, typename Index >
+void navierStokesSolverMonitor< Real, Index > :: setBetaMax( const Real& beta_max )
+{
+ this -> beta_max = beta_max;
+}
+
+template< typename Real, typename Index >
+void navierStokesSolverMonitor< Real, Index > :: setEpsSMax( const Real& eps_s_max )
+{
+ this -> eps_s_max = eps_s_max;
+}
+
+template< typename Real, typename Index >
+void navierStokesSolverMonitor< Real, Index > :: setUSMax( const Real& u_s_max )
+{
+ this -> u_s_max = u_s_max;
+}
+
+template< typename Real, typename Index >
+void navierStokesSolverMonitor< Real, Index > :: setUGMax( const Real& u_g_max )
+{
+ this -> u_g_max = u_g_max;
+}
+
+template< typename Real, typename Index >
+void navierStokesSolverMonitor< Real, Index > :: setEps_rho_g_max( const Real& eps_rho_g_max )
+{
+ this -> eps_rho_g_max = eps_rho_g_max;
+}
+
+template< typename Real, typename Index >
+void navierStokesSolverMonitor< Real, Index > :: setEps_rho_u1_g_max( const Real& eps_rho_u1_g_max )
+{
+ this -> eps_rho_u1_g_max = eps_rho_u1_g_max;
+}
+
+template< typename Real, typename Index >
+void navierStokesSolverMonitor< Real, Index > :: setEps_rho_u2_g_max( const Real& eps_rho_u2_g_max )
+{
+ this -> eps_rho_u2_g_max = eps_rho_u2_g_max;
+}
+
+template< typename Real, typename Index >
+void navierStokesSolverMonitor< Real, Index > :: setEps_rho_g_t_max( const Real& eps_rho_g_t_max )
+{
+ this -> eps_rho_g_t_max = eps_rho_g_t_max;
+}
+
+template< typename Real, typename Index >
+void navierStokesSolverMonitor< Real, Index > :: setEps_rho_u1_g_t_max( const Real& eps_rho_u1_g_t_max )
+{
+ this -> eps_rho_u1_g_t_max = eps_rho_u1_g_t_max;
+}
+
+template< typename Real, typename Index >
+void navierStokesSolverMonitor< Real, Index > :: setEps_rho_u2_g_t_max( const Real& eps_rho_u2_g_t_max )
+{
+ this -> eps_rho_u2_g_t_max = eps_rho_u2_g_t_max;
+}
+
+template< typename Real, typename Index >
+void navierStokesSolverMonitor< Real, Index > :: setEps_rho_s_max( const Real& eps_rho_s_max )
+{
+ this -> eps_rho_s_max = eps_rho_s_max;
+}
+
+template< typename Real, typename Index >
+void navierStokesSolverMonitor< Real, Index > :: setEps_rho_u1_s_max( const Real& eps_rho_u1_s_max )
+{
+ this -> eps_rho_u1_s_max = eps_rho_u1_s_max;
+}
+
+template< typename Real, typename Index >
+void navierStokesSolverMonitor< Real, Index > :: setEps_rho_u2_s_max( const Real& eps_rho_u2_s_max )
+{
+ this -> eps_rho_u2_s_max = eps_rho_u2_s_max;
+}
+
+template< typename Real, typename Index >
+void navierStokesSolverMonitor< Real, Index > :: setEps_rho_s_t_max( const Real& eps_rho_s_t_max )
+{
+ this -> eps_rho_s_t_max = eps_rho_s_t_max;
+}
+
+template< typename Real, typename Index >
+void navierStokesSolverMonitor< Real, Index > :: setEps_rho_u1_s_t_max( const Real& eps_rho_u1_s_t_max )
+{
+ this -> eps_rho_u1_s_t_max = eps_rho_u1_s_t_max;
+}
+
+template< typename Real, typename Index >
+void navierStokesSolverMonitor< Real, Index > :: setEps_rho_u2_s_t_max( const Real& eps_rho_u2_s_t_max )
+{
+ this -> eps_rho_u2_s_t_max = eps_rho_u2_s_t_max;
+}
+
+
+
+
+#endif /* TNLNAVIERSTOKESSOLVERMONITOR_IMPL_H_ */
diff --git a/examples/navier-stokes/navierStokesSolver_impl.h b/examples/navier-stokes/navierStokesSolver_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..7655a09fe18bb06efcdb8196059c89248534de80
--- /dev/null
+++ b/examples/navier-stokes/navierStokesSolver_impl.h
@@ -0,0 +1,556 @@
+/***************************************************************************
+ navierStokesSolver.impl.h - description
+ -------------------
+ begin : Jan 13, 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 NAVIERSTOKESSOLVER_IMPL_H_
+#define NAVIERSTOKESSOLVER_IMPL_H_
+
+
+#include "navierStokesSolver.h"
+#include <stdio.h>
+#include <iostream>
+#include <core/tnlString.h>
+#include <core/mfilename.h>
+#include <core/mfuncs.h>
+#include <core/tnlSharedVector.h>
+#include <solvers/ode/tnlMersonSolver.h>
+#include <legacy/mesh/tnlGridOld.h>
+
+
+#ifdef HAVE_CUDA
+
+#include <cuda.h>
+
+template< typename Real, typename Index >
+__device__ void computeVelocityFieldCuda( const Index size,
+ const Real R,
+ const Real T,
+ const Real* rho,
+ const Real* rho_u1,
+ const Real* rho_u2,
+ Real* u1,
+ Real* u2,
+ Real* p );
+#endif
+
+template< typename Mesh >
+navierStokesSolver< Mesh > :: navierStokesSolver()
+: mu( 0.0 ),
+ R( 0.0 ),
+ T( 0.0 ),
+ p_0( 0.0 ),
+ gravity( 0.0 )
+{
+ this -> rho. setName( "rho-g" );
+ this -> u1. setName( "u1-g");
+ this -> u2. setName( "u2-g" );
+ this -> p. setName( "p" );
+ this -> mesh. setName( "navier-stokes-mesh" );
+ this -> dofVector. setName( "navier-stokes-dof-vector" );
+}
+
+template< typename Mesh >
+bool navierStokesSolver< Mesh > :: init( const tnlParameterContainer& parameters )
+{
+ cout << "Initiating solver ... " << endl;
+
+ /****
+ * Set-up problem type
+ */
+ const tnlString& problemName = parameters. GetParameter< tnlString >( "problem-name" );
+ if( problemName == "riser" )
+ problem = riser;
+ if( problemName == "cavity" )
+ problem = cavity;
+
+ /****
+ * Set-up the geometry
+ */
+ tnlTuple< 2, RealType > proportions;
+ proportions. x() = parameters. GetParameter< double >( "width" );
+ proportions. y() = parameters. GetParameter< double >( "height" );
+ if( proportions. x() <= 0 )
+ {
+ cerr << "Error: width must be positive real number! It is " << proportions. x() << " now." << endl;
+ return false;
+ }
+ if( proportions. y() <= 0 )
+ {
+ cerr << "Error: height must be positive real number! It is " << proportions. y() << " now." << endl;
+ return false;
+ }
+ this -> mesh. setLowerCorner( tnlTuple< 2, RealType >( 0, 0 ) );
+ this -> mesh. setUpperCorner( proportions );
+
+ /****
+ * Set-up the space discretization
+ */
+ tnlTuple< 2, IndexType > meshes;
+ meshes. x() = parameters. GetParameter< int >( "x-size" );
+ meshes. y() = parameters. GetParameter< int >( "y-size" );
+ if( meshes. x() <= 0 )
+ {
+ cerr << "Error: x-size must be positive integer number! It is " << meshes. x() << " now." << endl;
+ return false;
+ }
+ if( meshes. y() <= 0 )
+ {
+ cerr << "Error: y-size must be positive integer number! It is " << meshes. y() << " now." << endl;
+ return false;
+ }
+ this -> mesh. setDimensions( meshes. y(), meshes. x() );
+ RealType hx = this -> mesh. getSpaceStep(). x();
+ RealType hy = this -> mesh. getSpaceStep(). y();
+ mesh. save( tnlString( "mesh.tnl" ) );
+
+ /****
+ * Set-up model coefficients
+ */
+ this -> p_0 = parameters. GetParameter< double >( "p0" );
+ this -> mu = parameters. GetParameter< double >( "mu");
+ this -> T = parameters. GetParameter< double >( "T" );
+ this -> R = parameters. GetParameter< double >( "R" );
+ this -> gravity = parameters. GetParameter< double >( "gravity" );
+ this -> maxInflowVelocity = parameters. GetParameter< double >( "max-inflow-velocity" );
+ this -> maxOutflowVelocity = parameters. GetParameter< double >( "max-outflow-velocity" );
+ this -> startUp = parameters. GetParameter< double >( "start-up" );
+
+ /****
+ * Set-up grid functions
+ */
+ const IndexType variablesNumber = 3;
+
+ rho. setSize( mesh. getDofs() );
+ u1. setSize( mesh. getDofs() );
+ u2. setSize( mesh. getDofs() );
+ p. setSize( mesh. getDofs() );
+
+ dofVector. setSize( variablesNumber * mesh. getDofs() );
+ rhsDofVector. setLike( dofVector );
+
+ /****
+ * Set-up numerical scheme
+ */
+ const tnlString& schemeName = parameters. GetParameter< tnlString >( "scheme" );
+ if( schemeName == "lax-fridrichs" )
+ this -> scheme = laxFridrichsEnm;
+
+ return true;
+}
+
+template< typename Mesh >
+bool navierStokesSolver< Mesh > :: setInitialCondition( const tnlParameterContainer& parameters )
+{
+ tnlSharedVector< RealType, DeviceType, IndexType > rho, rho_u1, rho_u2;
+ const IndexType& dofs = mesh. getDofs();
+ rho. bind( & dofVector. getData()[ 0 ], dofs );
+ rho_u1. bind( & dofVector. getData()[ dofs ], dofs );
+ rho_u2. bind( & dofVector. getData()[ 2 * dofs ], dofs );
+
+ rho. setValue( p_0 / ( this -> R * this -> T ) );
+ rho_u1. setValue( 0.0 );
+ rho_u2. setValue( 0.0 );
+
+ const IndexType& xSize = mesh. getDimensions(). x();
+ const IndexType& ySize = mesh. getDimensions(). y();
+ const RealType hx = mesh. getSpaceStep(). x();
+ const RealType hy = mesh. getSpaceStep(). y();
+
+ for( IndexType j = 0; j < ySize; j ++ )
+ for( IndexType i = 0; i < xSize; i ++ )
+ {
+ const IndexType c = mesh. getNodeIndex( j, i );
+ const RealType x = i * hx;
+ const RealType y = j * hy;
+
+ rho. setElement( c, p_0 / ( this -> R * this -> T ) );
+ rho_u1. setElement( c, 0.0 );
+ rho_u2. setElement( c, 0.0 );
+
+ }
+ return true;
+}
+
+template< typename Mesh >
+typename navierStokesSolver< Mesh > :: DofVectorType& navierStokesSolver< Mesh > :: getDofVector()
+{
+ return this -> dofVector;
+}
+
+template< typename Mesh >
+bool navierStokesSolver< Mesh > :: makeSnapshot( const RealType& t,
+ const IndexType step )
+{
+ cout << endl << "Writing output at time " << t << " step " << step << "." << endl;
+ tnlSharedVector< RealType, DeviceType, IndexType > rho, rho_u1, rho_u2,
+ rho_t, rho_u1_t, rho_u2_t;
+ const IndexType& dofs = mesh. getDofs();
+ rho. bind( & dofVector. getData()[ 0 ], dofs );
+ rho_u1. bind( & dofVector. getData()[ dofs ], dofs );
+ rho_u2. bind( & dofVector. getData()[ 2 * dofs ], dofs );
+ rho_t. bind( & this -> rhsDofVector. getData()[ 0 ], dofs );
+ rho_u1_t. bind( & this -> rhsDofVector. getData()[ dofs ], dofs );
+ rho_u2_t. bind( & this -> rhsDofVector. getData()[ 2 * dofs ], dofs );
+
+
+ updatePhysicalQuantities( rho, rho_u1, rho_u2 );
+ tnlVector< RealType, DeviceType, IndexType > u;
+ u. setLike( u1 );
+ for( IndexType i = 0; i < this -> u1. getSize(); i ++ )
+ u[ i ] = sqrt( u1[ i ] * u1[ i ] + u2[ i ] * u2[ i ] );
+ tnlString fileName;
+ /*FileNameBaseNumberEnding( "u-1-", step, 5, ".tnl", fileName );
+ if( ! this -> u1. save( fileName ) )
+ return false;
+ FileNameBaseNumberEnding( "u-2-", step, 5, ".tnl", fileName );
+ if( ! this -> u2. save( fileName ) )
+ return false;*/
+ /*FileNameBaseNumberEnding( "p-", step, 5, ".tnl", fileName );
+ if( ! p. save( fileName ) )
+ return false;*/
+ FileNameBaseNumberEnding( "u-", step, 5, ".tnl", fileName );
+ if( ! u. save( fileName ) )
+ return false;
+ return true;
+}
+
+template< typename Mesh >
+ template< typename Vector >
+void navierStokesSolver< Mesh > :: updatePhysicalQuantities( const Vector& rho,
+ const Vector& rho_u1,
+ const Vector& rho_u2 )
+{
+ if( DeviceType :: getDevice() == tnlHostDevice )
+ {
+ const IndexType& xSize = mesh. getDimensions(). x();
+ const IndexType& ySize = mesh. getDimensions(). y();
+
+ #ifdef HAVE_OPENMP
+ #pragma omp parallel for
+ #endif
+ for( IndexType j = 0; j < ySize; j ++ )
+ for( IndexType i = 0; i < xSize; i ++ )
+ {
+ IndexType c = mesh. getNodeIndex( j, i );
+ u1[ c ] = rho_u1[ c ] / rho[ c ];
+ u2[ c ] = rho_u2[ c ] / rho[ c ];
+ p[ c ] = rho[ c ] * this -> R * this -> T;
+ }
+ }
+#ifdef HAVE_CUDA
+ if( Mesh :: DeviceType :: getDevice() == tnlCudaDevice )
+ {
+ const int cudaBlockSize = 256;
+ const int maxGridSize = 65536;
+ const int cudaBlocks = ceil( u1. getSize() / cudaBlockSize );
+ const int gridNumbers = ceil( cudaBlocks / maxGridSize );
+ dim3 blockSize, gridSize;
+ blockSize. x = cudaBlockSize;
+ for( int grid = 0; grid < gridNumbers; grid ++ )
+ {
+
+ Index size( 0 );
+ if( grid < gridNumbers )
+ {
+ gridSize. x = maxGridSize;
+ size = gridSize. x * blockSize. x;
+ }
+ else
+ {
+ gridSize. x = cudaBlocks % maxGridSize;
+ size = u1. getSize() - ( gridNumbers - 1 ) * maxGridSize * cudaBlockSize;
+ }
+ Index startIndex = grid * maxGridSize * cudaBlockSize;
+ computeVelocityFieldCuda<<< blockSize, gridSize >>>( size,
+ this -> R,
+ this -> T,
+ & rho. getData()[ startIndex ],
+ & rho_u1. getData()[ startIndex ],
+ & rho_u2. getData()[ startIndex ],
+ & u1. getData()[ startIndex ],
+ & u2. getData()[ startIndex ],
+ & p. getData()[ startIdex ] );
+
+ }
+#endif
+}
+
+template< typename Mesh >
+void navierStokesSolver< Mesh > :: GetExplicitRHS( const RealType& time,
+ const RealType& tau,
+ DofVectorType& u,
+ DofVectorType& fu )
+{
+ tnlSharedVector< RealType, DeviceType, IndexType > rho, rho_u1, rho_u2,
+ rho_t, rho_u1_t, rho_u2_t;
+
+ const IndexType& dofs = mesh. getDofs();
+ rho. bind( & u. getData()[ 0 ], dofs );
+ rho_u1. bind( & u. getData()[ dofs ], dofs );
+ rho_u2. bind( & u. getData()[ 2 * dofs ], dofs );
+
+ rho_t. bind( & fu. getData()[ 0 ], dofs );
+ rho_u1_t. bind( & fu. getData()[ dofs ], dofs );
+ rho_u2_t. bind( & fu. getData()[ 2 * dofs ], dofs );
+
+ updatePhysicalQuantities( rho, rho_u1, rho_u2 );
+
+ /****
+ * Compute RHS
+ */
+ const IndexType& xSize = mesh. getDimensions(). x();
+ const IndexType& ySize = mesh. getDimensions(). y();
+ const RealType hx = mesh. getSpaceStep(). x();
+ const RealType hy = mesh. getSpaceStep(). y();
+ RealType startUpCoefficient( 1.0 );
+ if( this -> startUp != 0.0 )
+ startUpCoefficient = Min( ( RealType ) 1.0, time / this -> startUp );
+
+ if( DeviceType :: getDevice() == tnlHostDevice )
+ {
+ /****
+ * Set the boundary conditions.
+ * Speed: DBC on inlet, NBC on outlet, 0 DBC on walls.
+ * Density: NBS on inlet, DBC on outlet, NBC on walls
+ */
+ for( IndexType i = 0; i < xSize; i ++ )
+ {
+ const IndexType c1 = mesh. getNodeIndex( 0, i );
+ const IndexType c2 = mesh. getNodeIndex( 1, i );
+ const IndexType c3 = mesh. getNodeIndex( ySize - 1, i );
+ const IndexType c4 = mesh. getNodeIndex( ySize - 2, i );
+
+ RealType x = i * hx / mesh. getUpperCorner(). x();
+
+ /****
+ * Boundary conditions at the bottom and the top
+ */
+ if( problem == cavity )
+ {
+ this -> u1[ c1 ] = 0;
+ this -> u2[ c1 ] = 0;
+ this -> u1[ c3 ] = sin( M_PI * x ) * startUpCoefficient * this -> maxOutflowVelocity;
+ this -> u2[ c3 ] = 0;
+
+ rho[ c1 ] = rho[ c2 ];
+ rho[ c3 ] = rho[ c4 ];
+ //rho[ c3 ] = this -> p_0 / ( this -> R * this -> T );
+ }
+
+ rho_u1[ c1 ] = rho[ c1 ] * this -> u1[ c1 ];
+ rho_u2[ c1 ] = rho[ c1 ] * this -> u2[ c1 ];
+ rho_u1[ c3 ] = rho[ c3 ] * this -> u1[ c3 ];
+ rho_u2[ c3 ] = rho[ c3 ] * this -> u2[ c3 ];
+ }
+ for( IndexType j = 0; j < ySize; j ++ )
+ {
+ const IndexType c1 = mesh. getNodeIndex( j, 0 );
+ const IndexType c2 = mesh. getNodeIndex( j, 1 );
+ const IndexType c3 = mesh. getNodeIndex( j, xSize - 1 );
+ const IndexType c4 = mesh. getNodeIndex( j, xSize - 2 );
+
+ RealType y = j * hy / mesh. getUpperCorner(). y();
+
+ /****
+ * Boundary conditions on the left and right
+ */
+ if( problem == cavity )
+ {
+ this -> u1[ c1 ] = 0;
+ this -> u2[ c1 ] = 0;
+ this -> u1[ c3 ] = 0;
+ this -> u2[ c3 ] = 0;
+
+ rho[ c1 ] = rho[ c2 ];
+ rho[ c3 ] = rho[ c4 ];
+ }
+ rho_u1[ c1 ] = rho[ c1 ] * this -> u1[ c1 ];
+ rho_u2[ c1 ] = rho[ c1 ] * this -> u2[ c1 ];
+ rho_u1[ c3 ] = rho[ c3 ] * this -> u1[ c3 ];
+ rho_u2[ c3 ] = rho[ c3 ] * this -> u2[ c3 ];
+
+ }
+
+ const RealType c = sqrt( this -> R * this -> T );
+ #ifdef HAVE_OPENMP
+ #pragma omp parallel for
+ #endif
+ for( IndexType j = 0; j < ySize; j ++ )
+ for( IndexType i = 0; i < xSize; i ++ )
+ {
+ IndexType c = mesh. getNodeIndex( j, i );
+ if( i == 0 || j == 0 ||
+ i == xSize - 1 || j == ySize - 1 )
+ {
+ rho_t[ c ] = rho_u1_t[ c ] = rho_u2_t[ c ] = 0.0;
+ continue;
+ }
+
+ IndexType e = mesh. getNodeIndex( j, i + 1 );
+ IndexType w = mesh. getNodeIndex( j, i - 1 );
+ IndexType n = mesh. getNodeIndex( j + 1, i );
+ IndexType s = mesh. getNodeIndex( j - 1, i );
+
+ const RealType& u = this -> u1[ c ];
+ const RealType& v = this -> u2[ c ];
+ const RealType u_sqr = u * u;
+ const RealType v_sqr = v * v;
+ switch( this -> scheme )
+ {
+ case laxFridrichsEnm:
+ this -> laxFridrichsScheme. getExplicitRhs( mesh,
+ c,
+ rho,
+ rho_u1,
+ rho_u2,
+ u1,
+ u2,
+ rho_t,
+ rho_u1_t,
+ rho_u2_t,
+ 1.0 );
+ /****
+ * Remark: if fabs( u1_diff ) and fabs( u2_diff ) are used
+ * instead of u_diff in computeInterphaseFriction, higher beta
+ * can be used (up to 100 instead of 10).
+ */
+ rho_u1_t[ c ] += -( p[ e ] - p[ w ] ) / ( 2.0 * hx );
+ rho_u2_t[ c ] += -( p[ n ] - p[ s ] ) / ( 2.0 * hy );
+ - startUpCoefficient * this -> gravity * this -> rho[ c ];
+ break;
+ }
+
+ /***
+ * Add the viscosity term
+ */
+ rho_u1_t[ c ] += this -> mu *
+ ( ( u1[ e ] - 2.0 * u1[ c ] + u1[ w ] ) / ( hx * hx ) +
+ ( u1[ n ] - 2.0 * u1[ c ] + u1[ s ] ) / ( hy * hy ) );
+ rho_u2_t[ c ] += this -> mu *
+ ( ( u2[ e ] - 2.0 * u2[ c ] + u2[ w ] ) / ( hx * hx ) +
+ ( u2[ n ] - 2.0 * u2[ c ] + u2[ s ] ) / ( hy * hy ) );
+
+ }
+ }
+ if( DeviceType :: getDevice() == tnlCudaDevice )
+ {
+#ifdef HAVE_CUDA
+ const int cudaBlockSize = 256;
+ const int maxGridSize = 65536;
+ const int cudaBlocks = ceil( u1. getSize() / cudaBlockSize );
+ const int gridNumbers = ceil( cudaBlocks / maxGridSize );
+ dim3 blockSize, gridSize;
+ blockSize. x = cudaBlockSize;
+ for( int grid = 0; grid < gridNumbers; grid ++ )
+ {
+ if( grid < gridNumbers )
+ gridSize. x = maxGridSize;
+ else
+ gridSize. x = cudaBlocks % maxGridSize;
+
+ computeRhsCuda<<< blockSize, gridSize >>>
+ ( gridIdx,
+ xSize,
+ ySize,
+ hx,
+ hy,
+ rho. getData(),
+ rho_u1. getData(),
+ rho_u2. getData(),
+ u1. getData(),
+ u2. getData(),
+ p. getData(),
+ rho_t. getData(),
+ rho_u1_t. getData(),
+ rho_u2_t. getData() );
+
+ }
+#endif
+ }
+
+ rhsDofVector = fu;
+}
+
+#ifdef HAVE_CUDA
+template< typename Real, typename Index >
+__device__ void computeVelocityFieldCuda( const Index size,
+ const Real R,
+ const Real T,
+ const Real* rho,
+ const Real* rho_u1,
+ const Real* rho_u2,
+ Real* u1,
+ Real* u2,
+ Real* p )
+{
+ int globalThreadIdx = blockIdx. x * blockDim. x + threadIdx. x;
+ if( globalThreadIdx > size )
+ return;
+
+ u1[ globalThreadIdx ] = rho_u1[ globalThreadIdx ] / rho[ globalThreadIdx ];
+ u2[ globalThreadIdx ] = rho_u2[ globalThreadIdx ] / rho[ globalThreadIdx ];
+ p[ globalThreadIdx ] = rho[ globalThreadIdx ] * R * T;
+}
+#endif
+
+template< typename Mesh >
+tnlSolverMonitor< typename navierStokesSolver< Mesh > :: RealType,
+ typename navierStokesSolver< Mesh > :: IndexType >*
+ navierStokesSolver< Mesh > :: getSolverMonitor()
+{
+ return &solverMonitor;
+}
+
+template< typename Mesh >
+tnlString navierStokesSolver< Mesh > :: getTypeStatic()
+{
+ return tnlString( "navierStokesSolver< " ) +
+ Mesh :: getTypeStatic() + " >";
+}
+
+template< typename Mesh >
+tnlString navierStokesSolver< Mesh > :: getPrologHeader() const
+{
+ return tnlString( "Navier-Stokes Problem Solver" );
+}
+
+template< typename Mesh >
+void navierStokesSolver< Mesh > :: writeProlog( tnlLogger& logger,
+ const tnlParameterContainer& parameters ) const
+{
+ logger. WriteParameter< tnlString >( "Problem name:", "problem-name", parameters );
+ const tnlString& problemName = parameters. GetParameter< tnlString >( "problem-name" );
+ if( problemName == "cavity" )
+ {
+ logger. WriteParameter< double >( "Max. inflow velocity:", "max-inflow-velocity", parameters, 1 );
+ }
+ logger. WriteParameter< double >( "Viscosity:", "mu", parameters );
+ logger. WriteParameter< double >( "Temperature:", "T", parameters );
+ logger. WriteParameter< double >( "Gass constant:", "R", parameters );
+ logger. WriteParameter< double >( "Pressure:", "p0", parameters );
+ logger. WriteParameter< double >( "Gravity:", "gravity", parameters );
+ logger. WriteSeparator();
+ logger. WriteParameter< double >( "Domain width:", mesh. getUpperCorner(). x() - mesh. getLowerCorner(). x() );
+ logger. WriteParameter< double >( "Domain height:", mesh. getUpperCorner(). y() - mesh. getLowerCorner(). y() );
+ logger. WriteSeparator();
+ logger. WriteParameter< tnlString >( "Space discretisation:", "scheme", parameters );
+ logger. WriteParameter< int >( "Meshes along x:", mesh. getDimensions(). x() );
+ logger. WriteParameter< int >( "Meshes along y:", mesh. getDimensions(). y() );
+ logger. WriteParameter< double >( "Space step along x:", mesh. getSpaceStep(). x() );
+ logger. WriteParameter< double >( "Space step along y:", mesh. getSpaceStep(). y() );
+}
+
+#endif /* NAVIERSTOKESSOLVER_IMPL_H_ */
diff --git a/examples/navier-stokes/share/examples/cavity b/examples/navier-stokes/share/examples/cavity
new file mode 100644
index 0000000000000000000000000000000000000000..3bf836a0074ce9c02acd0803f58d7a23dbad17b1
--- /dev/null
+++ b/examples/navier-stokes/share/examples/cavity
@@ -0,0 +1,64 @@
+#!/bin/bash
+
+MAX_INFLOW_VELOCITIES="1 2 3 4 5 6 7 8 9 10"
+MAX_INFLOW_VELOCITIES="1"
+
+MU="10 1 0.1 0.01 0.001 2.3263e-5"
+
+max_outflow_velocity="2"
+
+for max_inflow_velocity in $MAX_INFLOW_VELOCITIES;
+do
+
+ mkdir v-${max_inflow_velocity}
+ cd v-${max_inflow_velocity}
+
+ for mu in ${MU};
+ do
+ mkdir mu-${mu}
+ cd mu-${mu}
+
+ navier-stokes --problem-name cavity \
+ --max-inflow-velocity ${max_inflow_velocity} \
+ --max-outflow-velocity ${max_outflow_velocity} \
+ --width 1 \
+ --height 1 \
+ --mu ${mu} \
+ --R 287.4 \
+ --T 400.0 \
+ --p0 101325.0 \
+ --final-time 10.0 \
+ --snapshot-period 0.04 \
+ --scheme lax-fridrichs \
+ --x-size 51 \
+ --y-size 51 \
+ --discrete-solver merson \
+ --initial-tau 2.0e-6 \
+ --merson-adaptivity 1.0e-4 \
+ --verbose 2 \
+ --device host
+
+ tnl-view --mesh mesh.tnl --input-files u*tnl eps*tnl
+
+ make-png-from-gnuplot u-g 1 1 1.5
+ make-png-from-gnuplot u-s 1 1 1.5
+ make-png-from-gnuplot eps-s 1 1 1
+
+ rm *.gplt
+
+ mencoder "mf://u-g*.png" -mf fps=25 -o u-g.avi -ovc lavc -lavcopts vcodec=mpeg4
+ mencoder "mf://u-s*.png" -mf fps=25 -o u-s.avi -ovc lavc -lavcopts vcodec=mpeg4
+ mencoder "mf://eps-s*.png" -mf fps=25 -o eps-s.avi -ovc lavc -lavcopts vcodec=mpeg4
+
+ rm *.png
+
+ cp u-g.avi ../../../u-g-v-${max_inflow_velocity}-mu-g-${mu_g}-eps-s-${eps_s}.avi
+ cp u-s.avi ../../../u-s-v-${max_inflow_velocity}-mu-g-${mu_g}-eps-s-${eps_s}.avi
+ cp eps-s.avi ../../../eps-s-v-${max_inflow_velocity}-mu-g-${mu_g}-eps-s-${eps_s}.avi
+
+ cd ..
+
+ done
+
+ cd ..
+done
diff --git a/src/implementation/mesh/tnlGrid1D_impl.h b/src/implementation/mesh/tnlGrid1D_impl.h
index dd6bd5627badad30ef9e03f1a9f37136cc5cef8e..d5e207b50db0001f901dc9fe7e2530eabfde1439 100644
--- a/src/implementation/mesh/tnlGrid1D_impl.h
+++ b/src/implementation/mesh/tnlGrid1D_impl.h
@@ -85,9 +85,9 @@ const tnlTuple< 1, Index >& tnlGrid< 1, Real, Device, Index> :: getDimensions()
template< typename Real,
typename Device,
typename Index >
-void tnlGrid< 1, Real, Device, Index> :: setLowerCorner( const tnlTuple< 1, Real >& lowerCorner )
+void tnlGrid< 1, Real, Device, Index> :: setLowerCorner( const tnlTuple< 1, Real >& origin )
{
- this -> lowerCorner = lowerCorner;
+ this -> origin = origin;
}
template< typename Real,
@@ -95,15 +95,15 @@ template< typename Real,
typename Index >
const tnlTuple< 1, Real >& tnlGrid< 1, Real, Device, Index> :: getLowerCorner() const
{
- return this -> lowerCorner;
+ return this -> origin;
}
template< typename Real,
typename Device,
typename Index >
-void tnlGrid< 1, Real, Device, Index> :: setUpperCorner( const tnlTuple< 1, Real >& upperCorner )
+void tnlGrid< 1, Real, Device, Index> :: setUpperCorner( const tnlTuple< 1, Real >& proportions )
{
- this -> upperCorner = upperCorner;
+ this -> proportions = proportions;
}
template< typename Real,
@@ -111,7 +111,7 @@ template< typename Real,
typename Index >
const tnlTuple< 1, Real >& tnlGrid< 1, Real, Device, Index> :: getUpperCorner() const
{
- return this -> upperCorner;
+ return this -> proportions;
}
template< typename Real,
@@ -119,7 +119,7 @@ template< typename Real,
typename Index >
void tnlGrid< 1, Real, Device, Index> :: setSpaceStep( const tnlTuple< 1, Real >& spaceStep )
{
- this -> upperCorner. x() = this -> lowerCorner. x() +
+ this -> proportions. x() = this -> origin. x() +
this -> dimensions. x() *
spaceStep. x();
}
@@ -134,7 +134,7 @@ tnlTuple< 1, Real > tnlGrid< 1, Real, Device, Index> :: getSpaceStep() const
<< this -> getName() );
tnlTuple< 1, RealType > spaceStep;
spaceStep. x() =
- ( this -> upperCorner. x() - this -> lowerCorner. x() ) /
+ ( this -> proportions. x() - this -> origin. x() ) /
( Real ) ( this -> dimensions. x() - 1 );
return spaceStep;
}
@@ -166,8 +166,8 @@ bool tnlGrid< 1, Real, Device, Index> :: save( tnlFile& file ) const
{
if( ! tnlObject :: save( file ) )
return false;
- if( ! this -> lowerCorner. save( file ) ||
- ! this -> upperCorner. save( file ) ||
+ if( ! this -> origin. save( file ) ||
+ ! this -> proportions. save( file ) ||
! this -> dimensions. save( file ) )
{
cerr << "I was not able to save the domain description of the tnlGrid "
@@ -184,8 +184,8 @@ bool tnlGrid< 1, Real, Device, Index> :: load( tnlFile& file )
{
if( ! tnlObject :: load( file ) )
return false;
- if( ! this -> lowerCorner. load( file ) ||
- ! this -> upperCorner. load( file ) ||
+ if( ! this -> origin. load( file ) ||
+ ! this -> proportions. load( file ) ||
! this -> dimensions. load( file ) )
{
cerr << "I was not able to load the domain description of the tnlGrid "
diff --git a/src/implementation/mesh/tnlGrid2D_impl.h b/src/implementation/mesh/tnlGrid2D_impl.h
index 4dedfa79c22b1949b1ff40768e9c75f33658a783..b9efecc13094e51ab56c61cc30bb79f38f45b1d9 100644
--- a/src/implementation/mesh/tnlGrid2D_impl.h
+++ b/src/implementation/mesh/tnlGrid2D_impl.h
@@ -91,9 +91,9 @@ const tnlTuple< 2, Index >& tnlGrid< 2, Real, Device, Index> :: getDimensions()
template< typename Real,
typename Device,
typename Index >
-void tnlGrid< 2, Real, Device, Index> :: setLowerCorner( const tnlTuple< 2, Real >& lowerCorner )
+void tnlGrid< 2, Real, Device, Index> :: setLowerCorner( const tnlTuple< 2, Real >& origin )
{
- this -> lowerCorner = lowerCorner;
+ this -> origin = origin;
}
template< typename Real,
@@ -101,15 +101,15 @@ template< typename Real,
typename Index >
const tnlTuple< 2, Real >& tnlGrid< 2, Real, Device, Index> :: getLowerCorner() const
{
- return this -> lowerCorner;
+ return this -> origin;
}
template< typename Real,
typename Device,
typename Index >
-void tnlGrid< 2, Real, Device, Index> :: setUpperCorner( const tnlTuple< 2, Real >& upperCorner )
+void tnlGrid< 2, Real, Device, Index> :: setUpperCorner( const tnlTuple< 2, Real >& proportions )
{
- this -> upperCorner = upperCorner;
+ this -> proportions = proportions;
}
template< typename Real,
@@ -117,7 +117,7 @@ template< typename Real,
typename Index >
const tnlTuple< 2, Real >& tnlGrid< 2, Real, Device, Index> :: getUpperCorner() const
{
- return this -> upperCorner;
+ return this -> proportions;
}
template< typename Real,
@@ -125,10 +125,10 @@ template< typename Real,
typename Index >
void tnlGrid< 2, Real, Device, Index> :: setSpaceStep( const tnlTuple< 2, Real >& spaceStep )
{
- this -> upperCorner. x() = this -> lowerCorner. x() +
+ this -> proportions. x() = this -> origin. x() +
this -> dimensions. x() *
spaceStep. x();
- this -> upperCorner. y() = this -> lowerCorner. y() +
+ this -> proportions. y() = this -> origin. y() +
this -> dimensions. y() *
spaceStep. y();
@@ -148,10 +148,10 @@ tnlTuple< 2, Real > tnlGrid< 2, Real, Device, Index> :: getSpaceStep() const
tnlTuple< 2, RealType > spaceStep;
spaceStep. x() =
- ( this -> upperCorner. x() - this -> lowerCorner. x() ) /
+ ( this -> proportions. x() - this -> origin. x() ) /
( Real ) ( this -> dimensions. x() - 1 );
spaceStep. y() =
- ( this -> upperCorner. y() - this -> lowerCorner. y() ) /
+ ( this -> proportions. y() - this -> origin. y() ) /
( Real ) ( this -> dimensions. y() - 1 );
return spaceStep;
}
@@ -204,8 +204,8 @@ bool tnlGrid< 2, Real, Device, Index> :: save( tnlFile& file ) const
{
if( ! tnlObject :: save( file ) )
return false;
- if( ! this -> lowerCorner. save( file ) ||
- ! this -> upperCorner. save( file ) ||
+ if( ! this -> origin. save( file ) ||
+ ! this -> proportions. save( file ) ||
! this -> dimensions. save( file ) )
{
cerr << "I was not able to save the domain description of the tnlGrid "
@@ -222,8 +222,8 @@ bool tnlGrid< 2, Real, Device, Index> :: load( tnlFile& file )
{
if( ! tnlObject :: load( file ) )
return false;
- if( ! this -> lowerCorner. load( file ) ||
- ! this -> upperCorner. load( file ) ||
+ if( ! this -> origin. load( file ) ||
+ ! this -> proportions. load( file ) ||
! this -> dimensions. load( file ) )
{
cerr << "I was not able to load the domain description of the tnlGrid "
diff --git a/src/implementation/mesh/tnlGrid3D_impl.h b/src/implementation/mesh/tnlGrid3D_impl.h
index 5b8fdaa91b9f35f11ac20b6128e1831e326b8b93..71aeaacc6f6a456c0f8dcc116ed9f8b66066227a 100644
--- a/src/implementation/mesh/tnlGrid3D_impl.h
+++ b/src/implementation/mesh/tnlGrid3D_impl.h
@@ -93,9 +93,9 @@ const tnlTuple< 3, Index >& tnlGrid< 3, Real, Device, Index> :: getDimensions()
template< typename Real,
typename Device,
typename Index >
-void tnlGrid< 3, Real, Device, Index> :: setLowerCorner( const tnlTuple< 3, Real >& lowerCorner )
+void tnlGrid< 3, Real, Device, Index> :: setLowerCorner( const tnlTuple< 3, Real >& origin )
{
- this -> lowerCorner = lowerCorner;
+ this -> origin = origin;
}
template< typename Real,
@@ -103,15 +103,15 @@ template< typename Real,
typename Index >
const tnlTuple< 3, Real >& tnlGrid< 3, Real, Device, Index> :: getLowerCorner() const
{
- return this -> lowerCorner;
+ return this -> origin;
}
template< typename Real,
typename Device,
typename Index >
-void tnlGrid< 3, Real, Device, Index> :: setUpperCorner( const tnlTuple< 3, Real >& upperCorner )
+void tnlGrid< 3, Real, Device, Index> :: setUpperCorner( const tnlTuple< 3, Real >& proportions )
{
- this -> upperCorner = upperCorner;
+ this -> proportions = proportions;
}
template< typename Real,
@@ -119,7 +119,7 @@ template< typename Real,
typename Index >
const tnlTuple< 3, Real >& tnlGrid< 3, Real, Device, Index> :: getUpperCorner() const
{
- return this -> upperCorner;
+ return this -> proportions;
}
template< typename Real,
@@ -127,13 +127,13 @@ template< typename Real,
typename Index >
void tnlGrid< 3, Real, Device, Index> :: setSpaceStep( const tnlTuple< 3, Real >& spaceStep )
{
- this -> upperCorner. x() = this -> lowerCorner. x() +
+ this -> proportions. x() = this -> origin. x() +
this -> dimensions. x() *
spaceStep. x();
- this -> upperCorner. y() = this -> lowerCorner. y() +
+ this -> proportions. y() = this -> origin. y() +
this -> dimensions. y() *
spaceStep. y();
- this -> upperCorner. z() = this -> lowerCorner. z() +
+ this -> proportions. z() = this -> origin. z() +
this -> dimensions. z() *
spaceStep. z();
}
@@ -155,13 +155,13 @@ tnlTuple< 3, Real > tnlGrid< 3, Real, Device, Index> :: getSpaceStep() const
tnlTuple< 3, RealType > spaceStep;
spaceStep. x() =
- ( this -> upperCorner. x() - this -> lowerCorner. x() ) /
+ ( this -> proportions. x() - this -> origin. x() ) /
( Real ) ( this -> dimensions. x() - 1 );
spaceStep. y() =
- ( this -> upperCorner. y() - this -> lowerCorner. y() ) /
+ ( this -> proportions. y() - this -> origin. y() ) /
( Real ) ( this -> dimensions. y() - 1 );
spaceStep. z() =
- ( this -> upperCorner. z() - this -> lowerCorner. z() ) /
+ ( this -> proportions. z() - this -> origin. z() ) /
( Real ) ( this -> dimensions. z() - 1 );
return spaceStep;
@@ -203,8 +203,8 @@ bool tnlGrid< 3, Real, Device, Index> :: save( tnlFile& file ) const
{
if( ! tnlObject :: save( file ) )
return false;
- if( ! this -> lowerCorner. save( file ) ||
- ! this -> upperCorner. save( file ) ||
+ if( ! this -> origin. save( file ) ||
+ ! this -> proportions. save( file ) ||
! this -> dimensions. save( file ) )
{
cerr << "I was not able to save the domain description of the tnlGrid "
@@ -221,8 +221,8 @@ bool tnlGrid< 3, Real, Device, Index> :: load( tnlFile& file )
{
if( ! tnlObject :: load( file ) )
return false;
- if( ! this -> lowerCorner. load( file ) ||
- ! this -> upperCorner. load( file ) ||
+ if( ! this -> origin. load( file ) ||
+ ! this -> proportions. load( file ) ||
! this -> dimensions. load( file ) )
{
cerr << "I was not able to load the domain description of the tnlGrid "
diff --git a/src/legacy/mesh/implementation/tnlGrid1D_impl.h b/src/legacy/mesh/implementation/tnlGrid1D_impl.h
index 39a8645c89d81b56609d270e26a75ef061ec8cef..a61f274ea275b7b603df0583c6d739f554472957 100644
--- a/src/legacy/mesh/implementation/tnlGrid1D_impl.h
+++ b/src/legacy/mesh/implementation/tnlGrid1D_impl.h
@@ -86,17 +86,17 @@ bool tnlGridOld< 1, Real, Device, Index > :: setDimensions( const tnlTuple< 1, I
}
template< typename Real, typename Device, typename Index >
-bool tnlGridOld< 1, Real, Device, Index > :: setDomain( const tnlTuple< 1, Real >& lowerCorner,
- const tnlTuple< 1, Real >& upperCorner )
+bool tnlGridOld< 1, Real, Device, Index > :: setDomain( const tnlTuple< 1, Real >& origin,
+ const tnlTuple< 1, Real >& proportions )
{
- if( lowerCorner >= upperCorner )
+ if( origin >= proportions )
{
cerr << "Wrong parameters for the grid domain of " << this -> getName() << ". The low corner must by smaller than the high corner." << endl
- << "lowerCorner = " << lowerCorner << endl << "upperCorner = " << upperCorner << endl;
+ << "origin = " << origin << endl << "proportions = " << proportions << endl;
return false;
}
- domainLowerCorner = lowerCorner;
- domainUpperCorner = upperCorner;
+ domainLowerCorner = origin;
+ domainUpperCorner = proportions;
spaceSteps[ 0 ] = ( domainUpperCorner[ 0 ] - domainLowerCorner[ 0 ] ) / ( Real ) ( this -> getDimensions()[ 0 ] - 1 );
return true;
}
diff --git a/src/legacy/mesh/implementation/tnlGrid2D_impl.h b/src/legacy/mesh/implementation/tnlGrid2D_impl.h
index b9e2b7f93c94b1bc7e17ff25045b49dd84c8df83..ccd918cf20286ecc3b9b5448a84c2399444abc01 100644
--- a/src/legacy/mesh/implementation/tnlGrid2D_impl.h
+++ b/src/legacy/mesh/implementation/tnlGrid2D_impl.h
@@ -92,17 +92,17 @@ bool tnlGridOld< 2, Real, Device, Index > :: setDimensions( const tnlTuple< 2, I
}
template< typename Real, typename Device, typename Index >
-bool tnlGridOld< 2, Real, Device, Index > :: setDomain( const tnlTuple< 2, Real >& lowerCorner,
- const tnlTuple< 2, Real >& upperCorner )
+bool tnlGridOld< 2, Real, Device, Index > :: setDomain( const tnlTuple< 2, Real >& origin,
+ const tnlTuple< 2, Real >& proportions )
{
- if( lowerCorner >= upperCorner )
+ if( origin >= proportions )
{
cerr << "Wrong parameters for the grid domain of " << this -> getName() << ". The low corner must by smaller than the high corner." << endl
- << "lowerCorner = " << lowerCorner << endl << "upperCorner = " << upperCorner << endl;
+ << "origin = " << origin << endl << "proportions = " << proportions << endl;
return false;
}
- domainLowerCorner = lowerCorner;
- domainUpperCorner = upperCorner;
+ domainLowerCorner = origin;
+ domainUpperCorner = proportions;
spaceSteps[ 0 ] = ( domainUpperCorner[ 0 ] - domainLowerCorner[ 0 ] ) / ( Real ) ( this -> getDimensions()[ 0 ] - 1 );
spaceSteps[ 1 ] = ( domainUpperCorner[ 1 ] - domainLowerCorner[ 1 ] ) / ( Real ) ( this -> getDimensions()[ 1 ] - 1 );
return true;
diff --git a/src/legacy/mesh/implementation/tnlGrid3D_impl.h b/src/legacy/mesh/implementation/tnlGrid3D_impl.h
index c62452df59e997b7bd8e3ea467fcf6f3a5c5a340..f25c8430d2cd1dd330aa7158bbe3c2bd72742ed2 100644
--- a/src/legacy/mesh/implementation/tnlGrid3D_impl.h
+++ b/src/legacy/mesh/implementation/tnlGrid3D_impl.h
@@ -94,17 +94,17 @@ bool tnlGridOld< 3,Real, Device, Index > :: setDimensions( const tnlTuple< 3,Ind
}
template< typename Real, typename Device, typename Index >
-bool tnlGridOld< 3,Real, Device, Index > :: setDomain( const tnlTuple< 3,Real >& lowerCorner,
- const tnlTuple< 3,Real >& upperCorner )
+bool tnlGridOld< 3,Real, Device, Index > :: setDomain( const tnlTuple< 3,Real >& origin,
+ const tnlTuple< 3,Real >& proportions )
{
- if( lowerCorner >= upperCorner )
+ if( origin >= proportions )
{
cerr << "Wrong parameters for the grid domain of " << this -> getName() << ". The low corner must by smaller than the high corner." << endl
- << "lowerCorner = " << lowerCorner << endl << "upperCorner = " << upperCorner << endl;
+ << "origin = " << origin << endl << "proportions = " << proportions << endl;
return false;
}
- domainLowerCorner = lowerCorner;
- domainUpperCorner = upperCorner;
+ domainLowerCorner = origin;
+ domainUpperCorner = proportions;
spaceSteps[ 0 ] = ( domainUpperCorner[ 0 ] - domainLowerCorner[ 0 ] ) / ( Real ) ( this -> getDimensions()[ 0 ] - 1 );
spaceSteps[ 1 ] = ( domainUpperCorner[ 1 ] - domainLowerCorner[ 1 ] ) / ( Real ) ( this -> getDimensions()[ 1 ] - 1 );
spaceSteps[ 2 ] = ( domainUpperCorner[ 2 ] - domainLowerCorner[ 2 ] ) / ( Real ) ( this -> getDimensions()[ 2 ] - 1 );
diff --git a/src/mesh/tnlGrid.h b/src/mesh/tnlGrid.h
index cc4531018cd9405fc46159dddffe02bea061090b..cc2c67527ed699dc1bfcbac86a73d0439869538a 100644
--- a/src/mesh/tnlGrid.h
+++ b/src/mesh/tnlGrid.h
@@ -54,11 +54,11 @@ class tnlGrid< 1, Real, Device, Index> : public tnlObject
const tnlTuple< 1, Index >& getDimensions() const;
- void setLowerCorner( const tnlTuple< 1, Real >& lowerCorner );
+ void setLowerCorner( const tnlTuple< 1, Real >& origin );
const tnlTuple< 1, Real >& getLowerCorner() const;
- void setUpperCorner( const tnlTuple< 1, Real >& upperCorner );
+ void setUpperCorner( const tnlTuple< 1, Real >& proportions );
const tnlTuple< 1, Real >& getUpperCorner() const;
@@ -89,7 +89,7 @@ class tnlGrid< 1, Real, Device, Index> : public tnlObject
tnlTuple< 1, IndexType > dimensions;
- tnlTuple< 1, RealType > lowerCorner, upperCorner;
+ tnlTuple< 1, RealType > origin, proportions;
IndexType dofs;
@@ -119,11 +119,11 @@ class tnlGrid< 2, Real, Device, Index> : public tnlObject
const tnlTuple< 2, Index >& getDimensions() const;
- void setLowerCorner( const tnlTuple< 2, Real >& lowerCorner );
+ void setLowerCorner( const tnlTuple< 2, Real >& origin );
const tnlTuple< 2, Real >& getLowerCorner() const;
- void setUpperCorner( const tnlTuple< 2, Real >& upperCorner );
+ void setUpperCorner( const tnlTuple< 2, Real >& proportions );
const tnlTuple< 2, Real >& getUpperCorner() const;
@@ -158,7 +158,7 @@ class tnlGrid< 2, Real, Device, Index> : public tnlObject
tnlTuple< 2, IndexType > dimensions;
- tnlTuple< 2, RealType > lowerCorner, upperCorner;
+ tnlTuple< 2, RealType > origin, proportions;
IndexType dofs;
@@ -188,11 +188,11 @@ class tnlGrid< 3, Real, Device, Index> : public tnlObject
const tnlTuple< 3, Index >& getDimensions() const;
- void setLowerCorner( const tnlTuple< 3, Real >& lowerCorner );
+ void setLowerCorner( const tnlTuple< 3, Real >& origin );
const tnlTuple< 3, Real >& getLowerCorner() const;
- void setUpperCorner( const tnlTuple< 3, Real >& upperCorner );
+ void setUpperCorner( const tnlTuple< 3, Real >& proportions );
const tnlTuple< 3, Real >& getUpperCorner() const;
@@ -223,7 +223,7 @@ class tnlGrid< 3, Real, Device, Index> : public tnlObject
tnlTuple< 3, IndexType > dimensions;
- tnlTuple< 3, RealType > lowerCorner, upperCorner;
+ tnlTuple< 3, RealType > origin, proportions;
IndexType dofs;