From a0069db33399a382e2ebe714b9ab2dc526caf804 Mon Sep 17 00:00:00 2001
From: Tomas Sobotik <sobotto4@fjfi.cvut.cz>
Date: Tue, 2 Dec 2014 08:01:09 +0100
Subject: [PATCH] Alpha version of parallel solver.
---
examples/CMakeLists.txt | 1 +
.../hamilton-jacobi-parallel/CMakeLists.txt | 17 +
.../MainBuildConfig.h | 64 +++
examples/hamilton-jacobi-parallel/Makefile | 41 ++
examples/hamilton-jacobi-parallel/main.cpp | 53 ++
.../parallelEikonalConfig.h | 45 ++
.../tnlParallelEikonalSolver.h | 84 ++++
.../tnlParallelEikonalSolver_impl.h | 469 ++++++++++++++++++
.../operators/godunov-eikonal/CMakeLists.txt | 5 +-
.../parallelGodunovEikonal1D_impl.h | 170 +++++++
.../parallelGodunovEikonal2D_impl.h | 224 +++++++++
.../parallelGodunovEikonal3D_impl.h | 197 ++++++++
.../solvers/ode/tnlEulerSolver_impl.h | 4 +-
src/operators/godunov-eikonal/CMakeLists.txt | 1 +
.../godunov-eikonal/parallelGodunovEikonal.h | 212 ++++++++
15 files changed, 1584 insertions(+), 3 deletions(-)
create mode 100755 examples/hamilton-jacobi-parallel/CMakeLists.txt
create mode 100644 examples/hamilton-jacobi-parallel/MainBuildConfig.h
create mode 100644 examples/hamilton-jacobi-parallel/Makefile
create mode 100644 examples/hamilton-jacobi-parallel/main.cpp
create mode 100644 examples/hamilton-jacobi-parallel/parallelEikonalConfig.h
create mode 100644 examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver.h
create mode 100644 examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h
create mode 100644 src/implementation/operators/godunov-eikonal/parallelGodunovEikonal1D_impl.h
create mode 100644 src/implementation/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h
create mode 100644 src/implementation/operators/godunov-eikonal/parallelGodunovEikonal3D_impl.h
create mode 100644 src/operators/godunov-eikonal/parallelGodunovEikonal.h
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 09c7d2cf34..b5fa429c12 100755
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -3,3 +3,4 @@ add_subdirectory( simple-solver )
add_subdirectory( heat-equation )
add_subdirectory( navier-stokes )
add_subdirectory( hamilton-jacobi )
+add_subdirectory( hamilton-jacobi-parallel )
diff --git a/examples/hamilton-jacobi-parallel/CMakeLists.txt b/examples/hamilton-jacobi-parallel/CMakeLists.txt
new file mode 100755
index 0000000000..31bc1cfca6
--- /dev/null
+++ b/examples/hamilton-jacobi-parallel/CMakeLists.txt
@@ -0,0 +1,17 @@
+set( tnl_hamilton_jacobi_parallel_SOURCES
+ MainBuildConfig.h
+ tnlParallelEikonalSolver_impl.h
+ tnlParallelEikonalSolver.h
+ main.cpp
+ parallelEikonalConfig.h )
+
+ADD_EXECUTABLE(hamilton-jacobi-parallel${debugExt} main.cpp)
+target_link_libraries (hamilton-jacobi-parallel${debugExt} tnl${debugExt}-${tnlVersion} )
+
+INSTALL( TARGETS hamilton-jacobi-parallel${debugExt}
+ RUNTIME DESTINATION bin
+ PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE )
+
+INSTALL( FILES ${tnl_hamilton_jacobi_parallel_SOURCES}
+ Makefile
+ DESTINATION share/tnl-${tnlVersion}/examples/hamilton-jacobi-parallel )
diff --git a/examples/hamilton-jacobi-parallel/MainBuildConfig.h b/examples/hamilton-jacobi-parallel/MainBuildConfig.h
new file mode 100644
index 0000000000..8ece05b9dc
--- /dev/null
+++ b/examples/hamilton-jacobi-parallel/MainBuildConfig.h
@@ -0,0 +1,64 @@
+/***************************************************************************
+ MainBuildConfig.h - description
+ -------------------
+ begin : Jul 7, 2014
+ copyright : (C) 2014 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 MAINBUILDCONFIG_H_
+#define MAINBUILDCONFIG_H_
+
+#include <solvers/tnlConfigTags.h>
+
+class MainBuildConfig
+{
+ public:
+
+ static void print() { cerr << "MainBuildConfig" << endl; }
+};
+
+/****
+ * Turn off support for float and long double.
+ */
+template<> struct tnlConfigTagReal< MainBuildConfig, float > { enum { enabled = false }; };
+template<> struct tnlConfigTagReal< MainBuildConfig, long double > { enum { enabled = false }; };
+
+/****
+ * Turn off support for short int and long int indexing.
+ */
+template<> struct tnlConfigTagIndex< MainBuildConfig, short int >{ enum { enabled = false }; };
+template<> struct tnlConfigTagIndex< MainBuildConfig, long int >{ enum { enabled = false }; };
+
+/****
+ * Use of tnlGrid is enabled for allowed dimensions and Real, Device and Index types.
+ */
+template< int Dimensions, typename Real, typename Device, typename Index >
+ struct tnlConfigTagMesh< MainBuildConfig, tnlGrid< Dimensions, Real, Device, Index > >
+ { enum { enabled = tnlConfigTagDimensions< MainBuildConfig, Dimensions >::enabled &&
+ tnlConfigTagReal< MainBuildConfig, Real >::enabled &&
+ tnlConfigTagDevice< MainBuildConfig, Device >::enabled &&
+ tnlConfigTagIndex< MainBuildConfig, Index >::enabled }; };
+
+/****
+ * Please, chose your preferred time discretisation here.
+ */
+template<> struct tnlConfigTagTimeDiscretisation< MainBuildConfig, tnlExplicitTimeDiscretisationTag >{ enum { enabled = true }; };
+template<> struct tnlConfigTagTimeDiscretisation< MainBuildConfig, tnlSemiImplicitTimeDiscretisationTag >{ enum { enabled = false}; };
+template<> struct tnlConfigTagTimeDiscretisation< MainBuildConfig, tnlImplicitTimeDiscretisationTag >{ enum { enabled = false }; };
+
+/****
+ * Only the Runge-Kutta-Merson solver is enabled by default.
+ */
+template<> struct tnlConfigTagExplicitSolver< MainBuildConfig, tnlExplicitEulerSolverTag >{ enum { enabled = false }; };
+
+#endif /* MAINBUILDCONFIG_H_ */
diff --git a/examples/hamilton-jacobi-parallel/Makefile b/examples/hamilton-jacobi-parallel/Makefile
new file mode 100644
index 0000000000..6b0e04e313
--- /dev/null
+++ b/examples/hamilton-jacobi-parallel/Makefile
@@ -0,0 +1,41 @@
+TNL_VERSION=0.1
+TNL_INSTALL_DIR=${HOME}/local/lib
+TNL_INCLUDE_DIR=${HOME}/local/include/tnl-${TNL_VERSION}
+
+TARGET = hamiltonJacobiParallelSolver
+#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) -O3 $(OMP_FLAGS)
+LD_FLAGS = -L$(TNL_INSTALL_DIR) -ltnl-0.1 -lgomp
+
+SOURCES = main.cpp
+HEADERS =
+OBJECTS = main.o
+DIST = $(SOURCES) 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
+
+uninstall: $(TARGET)
+ rm -f $(INSTALL_DIR)/bin/$(TARGET)
+ rm -f $(CONFIG_FILE) $(INSTALL_DIR)/share
+
+$(TARGET): $(OBJECTS)
+ $(CXX) -o $(TARGET) $(OBJECTS) $(LD_FLAGS)
+
+%.o: %.cpp $(HEADERS)
+ $(CXX) -c -o $@ $(CXX_FLAGS) $<
+
+
diff --git a/examples/hamilton-jacobi-parallel/main.cpp b/examples/hamilton-jacobi-parallel/main.cpp
new file mode 100644
index 0000000000..1561662b30
--- /dev/null
+++ b/examples/hamilton-jacobi-parallel/main.cpp
@@ -0,0 +1,53 @@
+/***************************************************************************
+ main.cpp - description
+ -------------------
+ begin : Jul 8 , 2014
+ copyright : (C) 2014 by Tomas Sobotik
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * 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 "tnlParallelEikonalSolver.h"
+#include "parallelEikonalConfig.h"
+#include "MainBuildConfig.h"
+#include <solvers/tnlConfigTags.h>
+#include <operators/godunov-eikonal/parallelGodunovEikonal.h>
+#include <mesh/tnlGrid.h>
+
+typedef MainBuildConfig BuildConfig;
+
+int main( int argc, char* argv[] )
+{
+ tnlParameterContainer parameters;
+ tnlConfigDescription configDescription;
+ parallelEikonalConfig< BuildConfig >::configSetup( configDescription );
+
+ if( ! ParseCommandLine( argc, argv, configDescription, parameters ) )
+ return false;
+
+ //if (parameters.GetParameter <tnlString>("scheme") == "godunov")
+ //{
+ typedef parallelGodunovEikonalScheme< tnlGrid<2,double,tnlHost, int>, double, int > SchemeType;
+ tnlParallelEikonalSolver<SchemeType> solver;
+ if(!solver.init(parameters))
+ {
+ cerr << "Solver failed to initialize." << endl;
+ return EXIT_FAILURE;
+ }
+ cout << "-------------------------------------------------------------" << endl;
+ cout << "Starting solver loop..." << endl;
+ solver.run();
+ cout << "a" <<endl;
+ // }
+
+ return EXIT_SUCCESS;
+}
+
+
diff --git a/examples/hamilton-jacobi-parallel/parallelEikonalConfig.h b/examples/hamilton-jacobi-parallel/parallelEikonalConfig.h
new file mode 100644
index 0000000000..76784f7a38
--- /dev/null
+++ b/examples/hamilton-jacobi-parallel/parallelEikonalConfig.h
@@ -0,0 +1,45 @@
+/***************************************************************************
+ hamiltonJacobiProblemConfig.h - description
+ -------------------
+ begin : Oct 5, 2014
+ copyright : (C) 2014 by Tomas Sobotik
+ email :
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * 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 HAMILTONJACOBIPROBLEMCONFIG_H_
+#define HAMILTONJACOBIPROBLEMCONFIG_H_
+
+#include <config/tnlConfigDescription.h>
+
+template< typename ConfigTag >
+class parallelEikonalConfig
+{
+ public:
+ static void configSetup( tnlConfigDescription& config )
+ {
+ config.addDelimiter( "Parallel Eikonal solver settings:" );
+ config.addEntry < tnlString > ( "problem-name", "This defines particular problem.", "hamilton-jacobi-parallel" );
+ config.addEntry < tnlString > ( "scheme", "This defines scheme used for discretization.", "godunov" );
+ config.addEntryEnum( "godunov" );
+ config.addEntryEnum( "upwind" );
+ config.addRequiredEntry < tnlString > ( "initial-condition", "Initial condition for solver");
+ config.addEntry < tnlString > ( "mesh", "Name of mesh.", "mesh.tnl" );
+ config.addEntry < double > ( "epsilon", "This defines epsilon for smoothening of sign().", 0.0 );
+ config.addEntry < double > ( "delta", " Allowed difference on subgrid boundaries", 0.0 );
+ config.addRequiredEntry < double > ( "stop-time", " Final time for solver");
+ config.addRequiredEntry < double > ( "initial-tau", " initial tau for solver" );
+ config.addEntry < double > ( "cfl-condition", " CFL condition", 0.0 );
+ config.addEntry < int > ( "subgrid-size", "Subgrid size.", 16 );
+ }
+};
+
+#endif /* HAMILTONJACOBIPROBLEMCONFIG_H_ */
diff --git a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver.h b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver.h
new file mode 100644
index 0000000000..5278e8e800
--- /dev/null
+++ b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver.h
@@ -0,0 +1,84 @@
+/***************************************************************************
+ tnlParallelEikonalSolver.h - description
+ -------------------
+ begin : Nov 28 , 2014
+ copyright : (C) 2014 by Tomas Sobotik
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * 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 TNLPARALLELEIKONALSOLVER_H_
+#define TNLPARALLELEIKONALSOLVER_H_
+
+#include <config/tnlParameterContainer.h>
+#include <core/vectors/tnlVector.h>
+#include <core/vectors/tnlStaticVector.h>
+#include <core/tnlHost.h>
+#include <mesh/tnlGrid.h>
+#include <limits.h>
+
+
+template< typename Scheme,
+ typename RealType = double,
+ typename Device = tnlHost,
+ typename IndexType = int >
+class tnlParallelEikonalSolver
+{};
+
+template< typename Scheme>
+class tnlParallelEikonalSolver<Scheme, double, tnlHost, int >
+{
+public:
+
+ typedef tnlVector< double, tnlHost, int > VectorType;
+ typedef tnlVector< int, tnlHost, int > IntVectorType;
+ typedef tnlGrid< 2, double, tnlHost, int > MeshType;
+
+ tnlParallelEikonalSolver();
+ bool init( const tnlParameterContainer& parameters );
+ void run();
+
+private:
+
+ void synchronize();
+
+ int getOwner( int i) const;
+
+ int getSubgridValue( int i ) const;
+
+ void setSubgridValue( int i, int value );
+
+ int getBoundaryCondition( int i ) const;
+
+ void setBoundaryCondition( int i, int value );
+
+ void stretchGrid();
+
+ void contractGrid();
+
+ VectorType getSubgrid( const int i ) const;
+
+ void insertSubgrid( VectorType u, const int i );
+
+ VectorType runSubgrid( int boundaryCondition, VectorType u);
+
+
+ VectorType u0, work_u;
+ IntVectorType subgridValues, boundaryConditions;
+ MeshType mesh, subMesh;
+ Scheme scheme;
+ double delta, tau0, stopTime,cflCondition;
+ int gridRows, gridCols, currentStep, n;
+
+};
+
+#include "tnlParallelEikonalSolver_impl.h"
+
+#endif /* TNLPARALLELEIKONALSOLVER_H_ */
diff --git a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h
new file mode 100644
index 0000000000..f6e0b4a1a7
--- /dev/null
+++ b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h
@@ -0,0 +1,469 @@
+/***************************************************************************
+ tnlParallelEikonalSolver_impl.h - description
+ -------------------
+ begin : Nov 28 , 2014
+ copyright : (C) 2014 by Tomas Sobotik
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * 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 TNLPARALLELEIKONALSOLVER_IMPL_H_
+#define TNLPARALLELEIKONALSOLVER_IMPL_H_
+
+
+#include "tnlParallelEikonalSolver.h"
+
+template< typename Scheme>
+tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::tnlParallelEikonalSolver()
+{
+}
+
+template< typename Scheme>
+bool tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::init( const tnlParameterContainer& parameters )
+{
+ cout << "Initializating solver..." << endl;
+ const tnlString& meshLocation = parameters.GetParameter <tnlString>("mesh");
+ this->mesh.load( meshLocation );
+
+ this->n = parameters.GetParameter <int>("subgrid-size");
+
+ this->subMesh.setDimensions( this->n, this->n );
+ this->subMesh.setDomain( tnlStaticVector<2,double>(0.0, 0.0),
+ tnlStaticVector<2,double>(this->mesh.getHx()*(double)(this->n), this->mesh.getHy()*(double)(this->n)) );
+
+
+ const tnlString& initialCondition = parameters.GetParameter <tnlString>("initial-condition");
+ this->u0.load( initialCondition );
+
+ this->delta = parameters.GetParameter <double>("delta");
+ this->delta *= sqrt(this->mesh.getHx()*this->mesh.getHx() + this->mesh.getHy()*this->mesh.getHy());
+
+ this->tau0 = parameters.GetParameter <double>("initial-tau");
+ this->stopTime = parameters.GetParameter <double>("stop-time");
+
+ this->cflCondition = parameters.GetParameter <double>("cfl-condition");
+
+ stretchGrid();
+ this->stopTime /= (double)(this->gridCols);
+
+ cout << "Initializating scheme..." << endl;
+ if(!this->scheme.init(parameters))
+ {
+ cerr << "Scheme failed to initialize." << endl;
+ return false;
+ }
+ cout << "Scheme initialized." << endl;
+
+
+ VectorType* tmp = new VectorType[subgridValues.getSize()];
+ bool containsCurve = false;
+
+ for(int i = 0; i < this->subgridValues.getSize(); i++)
+ {
+
+ if(! tmp[i].setSize(this->n * this->n))
+ cout << "Could not allocate tmp["<< i <<"] array." << endl;
+ tmp[i] = getSubgrid(i);
+ containsCurve = false;
+
+ for(int j = 0; j < tmp[i].getSize(); j++)
+ {
+ if(tmp[i][0]*tmp[i][j] <= 0.0)
+ {
+ containsCurve = true;
+ //cout << tmp[i][0] << ":" << tmp[i][j] << endl;
+ j=tmp[i].getSize();
+ }
+
+ }
+ if(containsCurve)
+ {
+ //cout << "Computing initial SDF on subgrid " << i << "." << endl;
+ insertSubgrid(runSubgrid(0, tmp[i]), i);
+ setSubgridValue(i, 4);
+ //cout << "Computed initial SDF on subgrid " << i << "." << endl;
+ }
+ containsCurve = false;
+
+ }
+
+ this->currentStep = 1;
+ synchronize();
+ cout << "Solver initialized." << endl;
+ return true;
+}
+
+template< typename Scheme >
+void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::run()
+{
+ while (this->subgridValues.max() > 0)
+ {
+
+ for(int i = 0; i < this->subgridValues.getSize(); i++)
+ {
+ if(getSubgridValue(i) != INT_MAX && getSubgridValue(i) > 0)
+ {
+ if(getBoundaryCondition(i) & 1)
+ {
+ insertSubgrid( runSubgrid(1, getSubgrid(i)), i);
+ setBoundaryCondition(i, getBoundaryCondition(i) - 1);
+ }
+ if(getBoundaryCondition(i) & 2)
+ {
+ insertSubgrid( runSubgrid(2, getSubgrid(i)), i);
+ setBoundaryCondition(i, getBoundaryCondition(i) - 2);
+ }
+ if(getBoundaryCondition(i) & 4)
+ {
+ insertSubgrid( runSubgrid(4, getSubgrid(i)), i);
+ setBoundaryCondition(i, getBoundaryCondition(i) - 4);
+ }
+ if(getBoundaryCondition(i) & 8)
+ {
+ insertSubgrid( runSubgrid(8, getSubgrid(i)), i);
+ setBoundaryCondition(i, getBoundaryCondition(i) - 8);
+ }
+
+ setSubgridValue(i, getSubgridValue(i)-1);
+ }
+ }
+
+ synchronize();
+ }
+
+ contractGrid();
+ this->u0.save("u-00001.tnl");
+ cout << "Solver finished" << endl;
+
+}
+//north - 1, east - 2, west - 4, south - 8
+template< typename Scheme >
+void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::synchronize()
+{
+ cout << "Synchronizig..." << endl;
+ int tmp1, tmp2;
+ int grid1, grid2;
+
+ for(int j = 0; j < this->gridRows - 1; j++)
+ {
+ for (int i = 0; i < this->gridCols*this->n; i++)
+ {
+ tmp1 = this->gridCols*this->n*((this->n-1)+j*this->n) + i;
+ tmp2 = this->gridCols*this->n*((this->n)+j*this->n) + i;
+ grid1 = getSubgridValue(getOwner(tmp1));
+ grid2 = getSubgridValue(getOwner(tmp2));
+ if ((this->work_u[tmp1] < this->work_u[tmp2] - this->delta || grid2 == INT_MAX|| grid2 == -INT_MAX) && (grid1 != INT_MAX && grid1 != -INT_MAX))
+ {
+ this->work_u[tmp2] = this->work_u[tmp1];
+ if(getSubgridValue(getOwner(tmp2)) == INT_MAX)
+ {
+ setSubgridValue(getOwner(tmp2), -INT_MAX);
+ }
+ setBoundaryCondition(getOwner(tmp2), getBoundaryCondition(getOwner(tmp2))+1);
+ }
+ else if ((this->work_u[tmp1] > this->work_u[tmp2] + this->delta || grid1 == INT_MAX || grid1 == -INT_MAX) && (grid2 != INT_MAX && grid2 != -INT_MAX))
+ {
+ this->work_u[tmp1] = this->work_u[tmp2];
+ if(getSubgridValue(getOwner(tmp1)) == INT_MAX)
+ {
+ setSubgridValue(getOwner(tmp1), -INT_MAX);
+ }
+ setBoundaryCondition(getOwner(tmp1), getBoundaryCondition(getOwner(tmp1))+8);
+ }
+ /*else if ( (grid1 != INT_MAX && grid1 != -INT_MAX) || (grid2 != INT_MAX && grid2 != -INT_MAX))
+ {
+ if(getSubgridValue(getOwner(tmp2)) == INT_MAX)
+ {
+ setSubgridValue(getOwner(tmp2), -INT_MAX);
+ }
+ else if(getSubgridValue(getOwner(tmp1)) == INT_MAX)
+ {
+ setSubgridValue(getOwner(tmp1), -INT_MAX);
+ }
+ }*/
+ }
+
+ }
+
+
+ for(int i = 1; i < this->gridCols; i++)
+ {
+ for (int j = 0; j < this->gridRows*this->n; j++)
+ {
+ tmp1 = this->gridCols*this->n*j + i*this->n;
+ tmp2 = this->gridCols*this->n*j + i*this->n + 1;
+ if ((this->work_u[tmp1] < this->work_u[tmp2] - this->delta || grid2 == INT_MAX|| grid2 == -INT_MAX) && (grid1 != INT_MAX && grid1 != -INT_MAX))
+ {
+ this->work_u[tmp2] = this->work_u[tmp1];
+ if(getSubgridValue(getOwner(tmp2)) == INT_MAX)
+ {
+ setSubgridValue(getOwner(tmp2), -INT_MAX);
+ }
+ setBoundaryCondition(getOwner(tmp2), getBoundaryCondition(getOwner(tmp2))+2);
+ }
+ else if ((this->work_u[tmp1] > this->work_u[tmp2] + this->delta || grid1 == INT_MAX || grid1 == -INT_MAX) && (grid2 != INT_MAX && grid2 != -INT_MAX))
+ {
+ this->work_u[tmp1] = this->work_u[tmp2];
+ if(getSubgridValue(getOwner(tmp1)) == INT_MAX)
+ {
+ setSubgridValue(getOwner(tmp1), -INT_MAX);
+ }
+ setBoundaryCondition(getOwner(tmp1), getBoundaryCondition(getOwner(tmp1))+4);
+ }
+ /*else if ( (grid1 != INT_MAX && grid1 != -INT_MAX) || (grid2 != INT_MAX && grid2 != -INT_MAX))
+ {
+ if(getSubgridValue(getOwner(tmp2)) == INT_MAX)
+ {
+ setSubgridValue(getOwner(tmp2), -INT_MAX);
+ }
+ else if(getSubgridValue(getOwner(tmp1)) == INT_MAX)
+ {
+ setSubgridValue(getOwner(tmp1), -INT_MAX);
+ }
+ }*/
+ }
+
+ }
+
+ int stepValue = this->currentStep + 3;
+ this->currentStep++;
+ for (int i = 0; i < this->subgridValues.getSize(); i++)
+ {
+ if( getSubgridValue(i) == -INT_MAX )
+ setSubgridValue(i, stepValue);
+ }
+
+ cout << "Grid synchronized at step " << (this->currentStep - 1 ) << endl;
+
+}
+
+
+template< typename Scheme >
+int tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::getOwner(int i) const
+{
+
+ return (i / (this->gridCols*this->n*this->n))*this->gridCols + (i % (this->gridCols*this->n))/this->n;
+}
+
+template< typename Scheme >
+int tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::getSubgridValue( int i ) const
+{
+ return this->subgridValues[i];
+}
+
+template< typename Scheme >
+void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::setSubgridValue(int i, int value)
+{
+ this->subgridValues[i] = value;
+}
+
+template< typename Scheme >
+int tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::getBoundaryCondition( int i ) const
+{
+ return this->boundaryConditions[i];
+}
+
+template< typename Scheme >
+void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::setBoundaryCondition(int i, int value)
+{
+ this->boundaryConditions[i] = value;
+}
+
+template< typename Scheme >
+void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::stretchGrid()
+{
+ cout << "Stretching grid..." << endl;
+
+
+ this->gridCols = ceil((double)(this->mesh.getDimensions().x())/(double)(this->n));
+ this->gridRows = ceil((double)(this->mesh.getDimensions().y())/(double)(this->n));
+
+ this->subgridValues.setSize(this->gridCols*this->gridRows);
+ this->boundaryConditions.setSize(this->gridCols*this->gridRows);
+
+ for(int i = 0; i < this->subgridValues.getSize(); i++ )
+ {
+ this->subgridValues[i] = INT_MAX;
+ this->boundaryConditions[i] = 0;
+ }
+
+ int stretchedSize = this->n*this->n*this->gridCols*this->gridRows;
+
+ if(!this->work_u.setSize(stretchedSize))
+ cerr << "Could not allocate memory for stretched grid." << endl;
+
+
+ for(int i = 0; i < stretchedSize; i++)
+ {
+ int k = i/this->n - i/(this->n*this->gridCols) + this->n*this->gridCols*(i/(this->n*this->n*this->gridCols));
+ int j=(i % (this->n*this->gridCols)) - ( (this->mesh.getDimensions().x() - this->n)/(this->n - 1) + this->mesh.getDimensions().x());
+
+ if(j > 0)
+ k += j;
+
+ this->work_u[i] = this->u0[i-k];
+ }
+
+ cout << "Grid stretched." << endl;
+}
+
+template< typename Scheme >
+void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::contractGrid()
+{
+ cout << "Contracting grid..." << endl;
+ int stretchedSize = this->n*this->n*this->gridCols*this->gridRows;
+
+ for(int i = 0; i < stretchedSize; i++)
+ {
+ int k = i/this->n - i/(this->n*this->gridCols) + this->n*this->gridCols*(i/(this->n*this->n*this->gridCols));
+ int j=(i % (this->n*this->gridCols)) - ((this->mesh.getDimensions().x() - this->n)/(this->n - 1) + this->mesh.getDimensions().x());
+
+ if(!(j > 0))
+ this->u0[i-k] = this->work_u[i];
+ }
+
+ //MISSING FILL FOR BOTTOM
+ cout << "Grid contracted" << endl;
+}
+
+template< typename Scheme >
+typename tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::VectorType
+tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::getSubgrid( const int i ) const
+{
+ VectorType u;
+ u.setSize(this->n*this->n);
+
+ for( int j = 0; j < this->n*this->n; j++)
+ {
+ u[j] = this->work_u[(i / this->gridCols)*this->n*this->n*this->gridCols + (i % this->gridCols)*this->n + (j/this->n)*this->n*this->gridCols + (j % this->n)];
+ }
+ return u;
+}
+
+template< typename Scheme >
+void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::insertSubgrid( VectorType u, const int i )
+{
+ for( int j = 0; j < this->n*this->n; j++)
+ {
+ int index = (i / this->gridCols)*this->n*this->n*this->gridCols + (i % this->gridCols)*this->n + (j/this->n)*this->n*this->gridCols + (j % this->n);
+ //OMP LOCK index
+ if(this->work_u[index] > u[j])
+ this->work_u[index] = u[j];
+ //OMP UNLOCK index
+ }
+}
+
+template< typename Scheme >
+typename tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::VectorType
+tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::runSubgrid( int boundaryCondition, VectorType u)
+{
+ VectorType fu;
+
+ fu.setLike(u);
+ fu.setValue( 0.0 );
+
+/*
+ * Insert Euler-Solver Here
+ */
+
+ double time = 0.0;
+ double currentTau = this->tau0;
+ if( time + currentTau > this -> stopTime ) currentTau = this -> stopTime - time;
+ /*this->resetIterations();
+ this->setResidue( this->getConvergenceResidue() + 1.0 );
+
+ this -> refreshSolverMonitor();*/
+
+ /****
+ * Start the main loop
+ */
+ while( time < this->stopTime )
+ {
+ /****
+ * Compute the RHS
+ */
+ //this->problem->getExplicitRHS( time, currentTau, u, fu );
+
+ for( int i = 0; i < fu.getSize(); i ++ )
+ {
+ fu[ i ] = scheme.getValue( this->subMesh, i, this->subMesh.getCellCoordinates(i), u, time, boundaryCondition );
+ }
+
+
+ //RealType lastResidue = this->getResidue();
+ double maxResidue( 0.0 );
+ //cout << fu. absMax() << endl;
+ if( this -> cflCondition != 0.0 )
+ {
+ maxResidue = fu. absMax();
+ if( currentTau * maxResidue > this -> cflCondition )
+ {
+ currentTau *= 0.9;
+ continue;
+ }
+ }
+ //RealType newResidue( 0.0 );
+ //computeNewTimeLevel( u, currentTau, newResidue );
+ for( int i = 0; i < fu.getSize(); i ++ )
+ {
+ double add = currentTau * fu[ i ];
+ u[ i ] += add;
+ //localResidue += fabs( add );
+ }
+ //this->setResidue( newResidue );
+
+ /****
+ * When time is close to stopTime the new residue
+ * may be inaccurate significantly.
+ */
+ //if( currentTau + time == this -> stopTime ) this->setResidue( lastResidue );
+ time += currentTau;
+
+ /*if( ! this->nextIteration() )
+ return false;*/
+
+ /****
+ * Compute the new time step.
+ */
+ if( this -> cflCondition != 0.0 )
+ currentTau /= 0.95;
+
+ if( time + currentTau > this -> stopTime )
+ currentTau = this -> stopTime - time; //we don't want to keep such tau
+ //else this -> tau = currentTau;
+
+ //this -> refreshSolverMonitor();
+
+ /****
+ * Check stop conditions.
+ */
+ /*if( time >= this->getStopTime() ||
+ ( this -> getConvergenceResidue() != 0.0 && this->getResidue() < this -> getConvergenceResidue() ) )
+ {
+ //this -> refreshSolverMonitor();
+ return true;
+ }*/
+
+
+ }
+
+
+ VectorType solution;
+ solution.setLike(u);
+ for( int i = 0; i < u.getSize(); i ++ )
+ {
+ solution[i]=u[i];
+ }
+
+ return solution;
+}
+
+
+#endif /* TNLPARALLELEIKONALSOLVER_IMPL_H_ */
diff --git a/src/implementation/operators/godunov-eikonal/CMakeLists.txt b/src/implementation/operators/godunov-eikonal/CMakeLists.txt
index 6dfff4c79a..6e0c100777 100755
--- a/src/implementation/operators/godunov-eikonal/CMakeLists.txt
+++ b/src/implementation/operators/godunov-eikonal/CMakeLists.txt
@@ -1,6 +1,9 @@
set( headers godunovEikonal1D_impl.h
godunovEikonal2D_impl.h
- godunovEikonal3D_impl.h)
+ godunovEikonal3D_impl.h
+ parallelGodunovEikonal1D_impl.h
+ parallelGodunovEikonal2D_impl.h
+ parallelGodunovEikonal3D_impl.h)
SET( CURRENT_DIR ${CMAKE_SOURCE_DIR}/src/implementation/operators/godunov-eikonal )
diff --git a/src/implementation/operators/godunov-eikonal/parallelGodunovEikonal1D_impl.h b/src/implementation/operators/godunov-eikonal/parallelGodunovEikonal1D_impl.h
new file mode 100644
index 0000000000..899549178b
--- /dev/null
+++ b/src/implementation/operators/godunov-eikonal/parallelGodunovEikonal1D_impl.h
@@ -0,0 +1,170 @@
+/***************************************************************************
+ parallelGodunovEikonal1D_impl.h - description
+ -------------------
+ begin : Dec 1 , 2014
+ copyright : (C) 2014 by Tomas Sobotik
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * 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 GODUNOVEIKONAL1D_IMPL_H_
+#define GODUNOVEIKONAL1D_IMPL_H_
+
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+Real parallelGodunovEikonalScheme< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index > :: positivePart(const Real arg) const
+{
+ if(arg > 0.0)
+ return arg;
+ return 0.0;
+}
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+Real parallelGodunovEikonalScheme< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index > :: negativePart(const Real arg) const
+{
+ if(arg < 0.0)
+ return arg;
+ return 0.0;
+}
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+Real parallelGodunovEikonalScheme< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index > :: sign(const Real x, const Real eps) const
+{
+ if(x > eps)
+ return 1.0;
+ else if (x < -eps)
+ return (-1.0);
+
+ if(eps == 0.0)
+ return 0.0;
+
+ return sin( ( M_PI * x ) / (2 * eps) );
+}
+
+
+
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+bool parallelGodunovEikonalScheme< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index > :: init( const tnlParameterContainer& parameters )
+{
+
+
+ const tnlString& meshFile = parameters.GetParameter< tnlString >( "mesh" );
+ if( ! this->originalMesh.load( meshFile ) )
+ {
+ cerr << "I am not able to load the mesh from the file " << meshFile << "." << endl;
+ return false;
+ }
+
+
+ h = originalMesh.getHx();
+ cout << "h = " << h << endl;
+
+ epsilon = parameters. GetParameter< double >( "epsilon" );
+
+ if(epsilon != 0.0)
+ epsilon *=h;
+
+ /*dofVector. setSize( this->mesh.getDofs() );*/
+
+ return true;
+}
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString parallelGodunovEikonalScheme< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index > :: getType()
+{
+ return tnlString( "parallelGodunovEikonalScheme< " ) +
+ 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 parallelGodunovEikonalScheme< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >:: getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType boundaryCondition ) const
+{
+
+ RealType nabla, xb, xf, signui;
+
+ signui = sign(u[cellIndex],epsilon);
+
+ if(signui > 0.0)
+ {
+ xf = negativePart((u[mesh.getCellXSuccessor( cellIndex )] - u[cellIndex])/h);
+ xb = positivePart((u[cellIndex] - u[mesh.getCellXPredecessor( cellIndex )])/h);
+
+ if(xb + xf > 0.0)
+ xf = 0.0;
+ else
+ xb = 0.0;
+
+ nabla = sqrt (xf*xf + xb*xb );
+
+ return signui*(1.0 - nabla);
+ }
+ else if (signui < 0.0)
+ {
+ xf = positivePart((u[mesh.getCellXSuccessor( cellIndex )] - u[cellIndex])/h);
+ xb = negativePart((u[cellIndex] - u[mesh.getCellXPredecessor( cellIndex )])/h);
+
+ if(xb + xf > 0.0)
+ xb = 0.0;
+ else
+ xf = 0.0;
+
+ nabla = sqrt (xf*xf + xb*xb );
+
+ return signui*(1.0 - nabla);
+ }
+ else
+ {
+ return 0.0;
+ }
+
+}
+
+
+#endif /* GODUNOVEIKONAL1D_IMPL_H_ */
diff --git a/src/implementation/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h b/src/implementation/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h
new file mode 100644
index 0000000000..ba64ea36fb
--- /dev/null
+++ b/src/implementation/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h
@@ -0,0 +1,224 @@
+/***************************************************************************
+ parallelGodunovEikonal2D_impl.h - description
+ -------------------
+ begin : Dec 1 , 2014
+ copyright : (C) 2014 by Tomas Sobotik
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * 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 GODUNOVEIKONAL2D_IMPL_H_
+#define GODUNOVEIKONAL2D_IMPL_H_
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+Real parallelGodunovEikonalScheme< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index >:: positivePart(const Real arg) const
+{
+ if(arg > 0.0)
+ return arg;
+ return 0.0;
+}
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+Real parallelGodunovEikonalScheme< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index > :: negativePart(const Real arg) const
+{
+ if(arg < 0.0)
+ return arg;
+ return 0.0;
+}
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+Real parallelGodunovEikonalScheme< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index > :: sign(const Real x, const Real eps) const
+{
+ if(x > eps)
+ return 1.0;
+ else if (x < -eps)
+ return (-1.0);
+
+ if ( eps == 0.0)
+ return 0.0;
+
+ return sin((M_PI*x)/(2.0*eps));
+}
+
+
+
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+bool parallelGodunovEikonalScheme< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index > :: init( const tnlParameterContainer& parameters )
+{
+
+ const tnlString& meshFile = parameters.GetParameter< tnlString >( "mesh" );
+ MeshType originalMesh;
+ if( ! originalMesh.load( meshFile ) )
+ {
+ cerr << "I am not able to load the mesh from the file " << meshFile << "." << endl;
+ return false;
+ }
+
+
+
+
+ hx = originalMesh.getHx();
+ hy = originalMesh.getHy();
+
+ epsilon = parameters. GetParameter< double >( "epsilon" );
+
+ if(epsilon != 0.0)
+ epsilon *=sqrt( hx*hx + hy*hy );
+
+// dofVector. setSize( this->mesh.getDofs() );
+
+ return true;
+
+}
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString parallelGodunovEikonalScheme< 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 >
+template< typename Vector >
+#ifdef HAVE_CUDA
+__device__ __host__
+#endif
+Real parallelGodunovEikonalScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >:: getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType boundaryCondition ) const
+{
+
+ if ( (coordinates.x() == 0 && boundaryCondition == 4) ||
+ (coordinates.x() == mesh.getDimensions().x() - 1 && boundaryCondition == 2) ||
+ (coordinates.y() == 0 && boundaryCondition == 1) ||
+ (coordinates.y() == mesh.getDimensions().y() - 1 && boundaryCondition == 8) )
+ return 0.0;
+
+ RealType nabla, xb, xf, yb, yf, signui;
+
+ signui = sign(u[cellIndex],epsilon);
+
+ if(signui > 0.0)
+ {
+ if(coordinates.x() == mesh.getDimensions().x() - 1)
+ xf = negativePart((u[mesh.getCellXPredecessor( cellIndex )] - u[cellIndex])/hx);
+ else
+ xf = negativePart((u[mesh.getCellXSuccessor( cellIndex )] - u[cellIndex])/hx);
+
+ if(coordinates.x() == 0)
+ xb = positivePart((u[cellIndex] - u[mesh.getCellXSuccessor( cellIndex )])/hx);
+ else
+ xb = positivePart((u[cellIndex] - u[mesh.getCellXPredecessor( cellIndex )])/hx);
+
+ if(coordinates.y() == mesh.getDimensions().y() - 1)
+ yf = negativePart((u[mesh.getCellYPredecessor( cellIndex )] - u[cellIndex])/hy);
+ else
+ yf = negativePart((u[mesh.getCellYSuccessor( cellIndex )] - u[cellIndex])/hy);
+
+ if(coordinates.y() == 0)
+ yb = positivePart((u[cellIndex] - u[mesh.getCellYSuccessor( cellIndex )])/hy);
+ else
+ yb = positivePart((u[cellIndex] - u[mesh.getCellYPredecessor( cellIndex )])/hy);
+
+ if(xb + xf > 0.0)
+ xf = 0.0;
+ else
+ xb = 0.0;
+
+ if(yb + yf > 0.0)
+ yf = 0.0;
+ else
+ yb = 0.0;
+
+ nabla = sqrt (xf*xf + xb*xb + yf*yf + yb*yb );
+
+ return signui*(1.0 - nabla);
+ }
+ else if (signui < 0.0)
+ {
+ if(coordinates.x() == mesh.getDimensions().x() - 1)
+ xf = positivePart((u[mesh.getCellXPredecessor( cellIndex )] - u[cellIndex])/hx);
+ else
+ xf = positivePart((u[mesh.getCellXSuccessor( cellIndex )] - u[cellIndex])/hx);
+
+ if(coordinates.x() == 0)
+ xb = negativePart((u[cellIndex] - u[mesh.getCellXSuccessor( cellIndex )])/hx);
+ else
+ xb = negativePart((u[cellIndex] - u[mesh.getCellXPredecessor( cellIndex )])/hx);
+
+ if(coordinates.y() == mesh.getDimensions().y() - 1)
+ yf = positivePart((u[mesh.getCellYPredecessor( cellIndex )] - u[cellIndex])/hy);
+ else
+ yf = positivePart((u[mesh.getCellYSuccessor( cellIndex )] - u[cellIndex])/hy);
+
+ if(coordinates.y() == 0)
+ yb = negativePart((u[cellIndex] - u[mesh.getCellYSuccessor( cellIndex )])/hy);
+ else
+ yb = negativePart((u[cellIndex] - u[mesh.getCellYPredecessor( cellIndex )])/hy);
+
+ if(xb + xf > 0.0)
+ xb = 0.0;
+ else
+ xf = 0.0;
+
+ if(yb + yf > 0.0)
+ yb = 0.0;
+ else
+ yf = 0.0;
+
+ nabla = sqrt (xf*xf + xb*xb + yf*yf + yb*yb );
+
+ return signui*(1.0 - nabla);
+ }
+ else
+ {
+ return 0.0;
+ }
+
+}
+
+
+#endif /* GODUNOVEIKONAL2D_IMPL_H_ */
diff --git a/src/implementation/operators/godunov-eikonal/parallelGodunovEikonal3D_impl.h b/src/implementation/operators/godunov-eikonal/parallelGodunovEikonal3D_impl.h
new file mode 100644
index 0000000000..f160223ad1
--- /dev/null
+++ b/src/implementation/operators/godunov-eikonal/parallelGodunovEikonal3D_impl.h
@@ -0,0 +1,197 @@
+/***************************************************************************
+ parallelGodunovEikonal3D_impl.h - description
+ -------------------
+ begin : Dec 1 , 2014
+ copyright : (C) 2014 by Tomas Sobotik
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * 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 GODUNOVEIKONAL3D_IMPL_H_
+#define GODUNOVEIKONAL3D_IMPL_H_
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+Real parallelGodunovEikonalScheme< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, Index > :: positivePart(const Real arg) const
+{
+ if(arg > 0.0)
+ return arg;
+ return 0.0;
+}
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+Real parallelGodunovEikonalScheme< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, Index > :: negativePart(const Real arg) const
+{
+ if(arg < 0.0)
+ return arg;
+ return 0.0;
+}
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+Real parallelGodunovEikonalScheme< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, Index > :: sign(const Real x, const Real eps) const
+{
+ if(x > eps)
+ return 1.0;
+ else if (x < -eps)
+ return (-1.0);
+
+ if ( eps == 0.0)
+ return 0.0;
+
+ return sin((M_PI*x)/(2.0*eps));
+}
+
+
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+bool parallelGodunovEikonalScheme< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, Index > :: init( const tnlParameterContainer& parameters )
+{
+ const tnlString& meshFile = parameters.GetParameter< tnlString >( "mesh" );
+ if( ! this->originalMesh.load( meshFile ) )
+ {
+ cerr << "I am not able to load the mesh from the file " << meshFile << "." << endl;
+ return false;
+ }
+
+
+ hx = originalMesh.getHx();
+ hy = originalMesh.getHy();
+ hz = originalMesh.getHz();
+
+ epsilon = parameters. GetParameter< double >( "epsilon" );
+
+ if(epsilon != 0.0)
+ epsilon *=sqrt( hx*hx + hy*hy +hz*hz );
+
+ // dofVector. setSize( this->mesh.getDofs() );
+
+ return true;
+
+}
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString parallelGodunovEikonalScheme< 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 parallelGodunovEikonalScheme< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: getValue( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Vector& u,
+ const Real& time,
+ const IndexType boundaryCondition ) const
+{
+
+ RealType nabla, xb, xf, yb, yf, zb, zf, signui;
+
+ signui = sign(u[cellIndex],epsilon);
+
+ if(signui > 0.0)
+ {
+ xf = negativePart((u[mesh.getCellXSuccessor( cellIndex )] - u[cellIndex])/hx);
+ xb = positivePart((u[cellIndex] - u[mesh.getCellXPredecessor( cellIndex )])/hx);
+ yf = negativePart((u[mesh.getCellYSuccessor( cellIndex )] - u[cellIndex])/hy);
+ yb = positivePart((u[cellIndex] - u[mesh.getCellYPredecessor( cellIndex )])/hy);
+ zf = negativePart((u[mesh.getCellZSuccessor( cellIndex )] - u[cellIndex])/hz);
+ zb = positivePart((u[cellIndex] - u[mesh.getCellZPredecessor( cellIndex )])/hz);
+
+ if(xb + xf > 0.0)
+ xf = 0.0;
+ else
+ xb = 0.0;
+
+ if(yb + yf > 0.0)
+ yf = 0.0;
+ else
+ yb = 0.0;
+
+ if(zb + zf > 0.0)
+ zf = 0.0;
+ else
+ zb = 0.0;
+
+ nabla = sqrt (xf*xf + xb*xb + yf*yf + yb*yb + zf*zf + zb*zb );
+
+ return signui*(1.0 - nabla);
+ }
+ else if (signui < 0.0)
+ {
+ xf = positivePart((u[mesh.getCellXSuccessor( cellIndex )] - u[cellIndex])/hx);
+ xb = negativePart((u[cellIndex] - u[mesh.getCellXPredecessor( cellIndex )])/hx);
+ yf = positivePart((u[mesh.getCellYSuccessor( cellIndex )] - u[cellIndex])/hy);
+ yb = negativePart((u[cellIndex] - u[mesh.getCellYPredecessor( cellIndex )])/hy);
+ zf = positivePart((u[mesh.getCellZSuccessor( cellIndex )] - u[cellIndex])/hz);
+ zb = negativePart((u[cellIndex] - u[mesh.getCellZPredecessor( cellIndex )])/hz);
+
+ if(xb + xf > 0.0)
+ xb = 0.0;
+ else
+ xf = 0.0;
+
+ if(yb + yf > 0.0)
+ yb = 0.0;
+ else
+ yf = 0.0;
+
+ if(zb + zf > 0.0)
+ zb = 0.0;
+ else
+ zf = 0.0;
+
+ nabla = sqrt (xf*xf + xb*xb + yf*yf + yb*yb + zf*zf + zb*zb );
+
+ return signui*(1.0 - nabla);
+ }
+ else
+ {
+ return 0.0;
+ }
+
+}
+
+
+#endif /* GODUNOVEIKONAL3D_IMPL_H_ */
diff --git a/src/implementation/solvers/ode/tnlEulerSolver_impl.h b/src/implementation/solvers/ode/tnlEulerSolver_impl.h
index 3fb5f16170..eb2fdb28b1 100644
--- a/src/implementation/solvers/ode/tnlEulerSolver_impl.h
+++ b/src/implementation/solvers/ode/tnlEulerSolver_impl.h
@@ -81,8 +81,8 @@ bool tnlEulerSolver< Problem > :: solve( DofVectorType& u )
/****
* Set necessary parameters
*/
- RealType& time = this->time;
- RealType currentTau = this->tau;
+ double time = 0.0;
+ double currentTau = this->tau0;
if( time + currentTau > this -> getStopTime() ) currentTau = this -> getStopTime() - time;
if( currentTau == 0.0 ) return true;
this->resetIterations();
diff --git a/src/operators/godunov-eikonal/CMakeLists.txt b/src/operators/godunov-eikonal/CMakeLists.txt
index 79145e6fe7..554890c93e 100755
--- a/src/operators/godunov-eikonal/CMakeLists.txt
+++ b/src/operators/godunov-eikonal/CMakeLists.txt
@@ -1,4 +1,5 @@
SET( headers godunovEikonal.h
+ parallelGodunovEikonal.h
)
INSTALL( FILES ${headers} DESTINATION include/tnl-${tnlVersion}/operators/godunov-eikonal )
diff --git a/src/operators/godunov-eikonal/parallelGodunovEikonal.h b/src/operators/godunov-eikonal/parallelGodunovEikonal.h
new file mode 100644
index 0000000000..0834e22080
--- /dev/null
+++ b/src/operators/godunov-eikonal/parallelGodunovEikonal.h
@@ -0,0 +1,212 @@
+/***************************************************************************
+ godunov.h - description
+ -------------------
+ begin : Dec 1 , 2014
+ copyright : (C) 2014 by Tomas Sobotik
+ ***************************************************************************/
+
+/***************************************************************************
+ * *
+ * 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 GODUNOVEIKONAL_H_
+#define GODUNOVEIKONAL_H_
+
+#include <matrices/tnlCSRMatrix.h>
+#include <solvers/preconditioners/tnlDummyPreconditioner.h>
+#include <solvers/tnlSolverMonitor.h>
+#include <core/tnlLogger.h>
+#include <core/vectors/tnlVector.h>
+#include <core/vectors/tnlSharedVector.h>
+#include <core/mfilename.h>
+#include <mesh/tnlGrid.h>
+
+
+template< typename Mesh,
+ typename Real,
+ typename Index >
+class parallelGodunovEikonalScheme
+{
+};
+
+
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class parallelGodunovEikonalScheme< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index >
+{
+
+public:
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef tnlGrid< 1, Real, Device, Index > MeshType;
+ typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+
+
+
+ static tnlString getType();
+
+ RealType positivePart(const RealType arg) const;
+
+ RealType negativePart(const RealType arg) const;
+
+ RealType sign(const RealType x, const RealType eps) const;
+
+ 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 IndexType boundaryCondition) const;
+
+ bool init( const tnlParameterContainer& parameters );
+
+
+protected:
+
+ MeshType originalMesh;
+
+ DofVectorType dofVector;
+
+ RealType h;
+
+ RealType epsilon;
+
+
+};
+
+
+
+
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class parallelGodunovEikonalScheme< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index >
+{
+
+public:
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef tnlGrid< 2, Real, Device, Index > MeshType;
+ typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+
+ static tnlString getType();
+
+ RealType positivePart(const RealType arg) const;
+
+ RealType negativePart(const RealType arg) const;
+
+ RealType sign(const RealType x, const Real eps) const;
+
+ 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 IndexType boundaryCondition ) const;
+
+ bool init( const tnlParameterContainer& parameters );
+
+
+protected:
+
+ MeshType originalMesh;
+
+ DofVectorType dofVector;
+
+ RealType hx;
+ RealType hy;
+
+ RealType epsilon;
+
+
+};
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class parallelGodunovEikonalScheme< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, Index >
+{
+
+public:
+ typedef Real RealType;
+ typedef Device DeviceType;
+ typedef Index IndexType;
+ typedef tnlGrid< 3, Real, Device, Index > MeshType;
+ typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;
+ typedef typename MeshType::CoordinatesType CoordinatesType;
+
+
+ static tnlString getType();
+
+ RealType positivePart(const RealType arg) const;
+
+ RealType negativePart(const RealType arg) const;
+
+ RealType sign(const RealType x, const Real eps) const;
+
+ 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 IndexType boundaryCondition ) const;
+
+ bool init( const tnlParameterContainer& parameters );
+
+
+protected:
+
+ MeshType originalMesh;
+
+ DofVectorType dofVector;
+
+ RealType hx;
+ RealType hy;
+ RealType hz;
+
+ RealType epsilon;
+
+
+};
+
+
+
+#include <implementation/operators/godunov-eikonal/parallelGodunovEikonal1D_impl.h>
+#include <implementation/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h>
+#include <implementation/operators/godunov-eikonal/parallelGodunovEikonal3D_impl.h>
+
+
+#endif /* GODUNOVEIKONAL_H_ */
--
GitLab