From 6de09f94cea263339e16d65f826dbe1d4e36f53b Mon Sep 17 00:00:00 2001
From: Tomas Sobotik <sobotto4@fjfi.cvut.cz>
Date: Mon, 2 Nov 2015 10:34:38 +0100
Subject: [PATCH] Added fast sweeping method
---
examples/CMakeLists.txt | 1 +
examples/fast-sweeping/CMakeLists.txt | 22 +
examples/fast-sweeping/MainBuildConfig.h | 64 +++
examples/fast-sweeping/fastSweepingConfig.h | 36 ++
examples/fast-sweeping/main.cpp | 17 +
examples/fast-sweeping/main.cu | 17 +
examples/fast-sweeping/main.h | 61 +++
examples/fast-sweeping/tnlFastSweeping.h | 81 ++++
.../fast-sweeping/tnlFastSweeping2D_impl.h | 324 ++++++++++++++
examples/hamilton-jacobi-parallel/main.h | 4 +-
.../tnlParallelEikonalSolver_impl.h | 419 ++++++++++++------
.../godunov-eikonal/parallelGodunovEikonal.h | 21 +-
.../parallelGodunovEikonal2D_impl.h | 354 +++++++++------
.../parallelGodunovEikonal3D_impl.h | 310 ++++++++++---
14 files changed, 1388 insertions(+), 343 deletions(-)
create mode 100755 examples/fast-sweeping/CMakeLists.txt
create mode 100644 examples/fast-sweeping/MainBuildConfig.h
create mode 100644 examples/fast-sweeping/fastSweepingConfig.h
create mode 100644 examples/fast-sweeping/main.cpp
create mode 100644 examples/fast-sweeping/main.cu
create mode 100644 examples/fast-sweeping/main.h
create mode 100644 examples/fast-sweeping/tnlFastSweeping.h
create mode 100644 examples/fast-sweeping/tnlFastSweeping2D_impl.h
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 4ec29b8255..0d821366f8 100755
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -4,3 +4,4 @@ add_subdirectory( heat-equation )
add_subdirectory( navier-stokes )
#add_subdirectory( hamilton-jacobi )
add_subdirectory( hamilton-jacobi-parallel )
+add_subdirectory( fast-sweeping )
diff --git a/examples/fast-sweeping/CMakeLists.txt b/examples/fast-sweeping/CMakeLists.txt
new file mode 100755
index 0000000000..7a2963cc68
--- /dev/null
+++ b/examples/fast-sweeping/CMakeLists.txt
@@ -0,0 +1,22 @@
+set( tnl_fast_sweeping_SOURCES
+# MainBuildConfig.h
+# tnlFastSweeping2D_impl.h
+# tnlFastSweeping.h
+# fastSweepingConfig.h
+ main.cpp)
+
+
+IF( BUILD_CUDA )
+ CUDA_ADD_EXECUTABLE(fast-sweeping${debugExt} main.cu)
+ELSE( BUILD_CUDA )
+ ADD_EXECUTABLE(fast-sweeping${debugExt} main.cpp)
+ENDIF( BUILD_CUDA )
+target_link_libraries (fast-sweeping${debugExt} tnl${debugExt}-${tnlVersion} )
+
+
+INSTALL( TARGETS fast-sweeping${debugExt}
+ RUNTIME DESTINATION bin
+ PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE )
+
+#INSTALL( FILES ${tnl_fast_sweeping_SOURCES}
+# DESTINATION share/tnl-${tnlVersion}/examples/fast-sweeping )
diff --git a/examples/fast-sweeping/MainBuildConfig.h b/examples/fast-sweeping/MainBuildConfig.h
new file mode 100644
index 0000000000..8ece05b9dc
--- /dev/null
+++ b/examples/fast-sweeping/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/fast-sweeping/fastSweepingConfig.h b/examples/fast-sweeping/fastSweepingConfig.h
new file mode 100644
index 0000000000..6b8b7b9a91
--- /dev/null
+++ b/examples/fast-sweeping/fastSweepingConfig.h
@@ -0,0 +1,36 @@
+/***************************************************************************
+ fastSweepingConfig.h - description
+ -------------------
+ begin : Oct 15, 2015
+ copyright : (C) 2015 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 FASTSWEEPINGCONFIG_H_
+#define FASTSWEEPINGCONFIG_H_
+
+#include <config/tnlConfigDescription.h>
+
+template< typename ConfigTag >
+class fastSweepingConfig
+{
+ public:
+ static void configSetup( tnlConfigDescription& config )
+ {
+ config.addDelimiter( "Parallel Eikonal solver settings:" );
+ config.addEntry < tnlString > ( "problem-name", "This defines particular problem.", "fast-sweeping" );
+ config.addRequiredEntry < tnlString > ( "initial-condition", "Initial condition for solver");
+ config.addEntry < tnlString > ( "mesh", "Name of mesh.", "mesh.tnl" );
+ }
+};
+
+#endif /* FASTSWEEPINGCONFIG_H_ */
diff --git a/examples/fast-sweeping/main.cpp b/examples/fast-sweeping/main.cpp
new file mode 100644
index 0000000000..8849008ff6
--- /dev/null
+++ b/examples/fast-sweeping/main.cpp
@@ -0,0 +1,17 @@
+/***************************************************************************
+ main.cpp - description
+ -------------------
+ begin : Oct 15 , 2015
+ copyright : (C) 2015 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 "main.h"
diff --git a/examples/fast-sweeping/main.cu b/examples/fast-sweeping/main.cu
new file mode 100644
index 0000000000..8849008ff6
--- /dev/null
+++ b/examples/fast-sweeping/main.cu
@@ -0,0 +1,17 @@
+/***************************************************************************
+ main.cpp - description
+ -------------------
+ begin : Oct 15 , 2015
+ copyright : (C) 2015 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 "main.h"
diff --git a/examples/fast-sweeping/main.h b/examples/fast-sweeping/main.h
new file mode 100644
index 0000000000..c21b591358
--- /dev/null
+++ b/examples/fast-sweeping/main.h
@@ -0,0 +1,61 @@
+/***************************************************************************
+ main.h - description
+ -------------------
+ begin : Oct 15 , 2015
+ copyright : (C) 2015 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 "MainBuildConfig.h"
+#include "tnlFastSweeping.h"
+#include "fastSweepingConfig.h"
+#include <solvers/tnlConfigTags.h>
+
+#include <mesh/tnlGrid.h>
+#include <core/tnlDevice.h>
+#include <time.h>
+#include <ctime>
+
+typedef MainBuildConfig BuildConfig;
+
+int main( int argc, char* argv[] )
+{
+ time_t start;
+ time_t stop;
+ time(&start);
+ std::clock_t start2= std::clock();
+ tnlParameterContainer parameters;
+ tnlConfigDescription configDescription;
+ fastSweepingConfig< BuildConfig >::configSetup( configDescription );
+
+ if( ! parseCommandLine( argc, argv, configDescription, parameters ) )
+ return false;
+
+ tnlFastSweeping<tnlGrid<2,double,tnlHost, int>, double, int> solver;
+ if(!solver.init(parameters))
+ {
+ cerr << "Solver failed to initialize." << endl;
+ return EXIT_FAILURE;
+ }
+ cout << "-------------------------------------------------------------" << endl;
+ cout << "Starting solver..." << endl;
+ solver.run();
+
+
+ time(&stop);
+ cout << "Solver stopped..." << endl;
+ cout << endl;
+ cout << "Running time was: " << difftime(stop,start) << " .... " << (std::clock() - start2) / (double)(CLOCKS_PER_SEC) << endl;
+ return EXIT_SUCCESS;
+}
+
+
diff --git a/examples/fast-sweeping/tnlFastSweeping.h b/examples/fast-sweeping/tnlFastSweeping.h
new file mode 100644
index 0000000000..1e04925ecb
--- /dev/null
+++ b/examples/fast-sweeping/tnlFastSweeping.h
@@ -0,0 +1,81 @@
+/***************************************************************************
+ tnlFastSweeping.h - description
+ -------------------
+ begin : Oct 15 , 2015
+ copyright : (C) 2015 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 TNLFASTSWEEPING_H_
+#define TNLFASTSWEEPING_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>
+#include <core/tnlDevice.h>
+#include <ctime>
+
+
+
+
+template< typename Mesh,
+ typename Real,
+ typename Index >
+class tnlFastSweeping
+{};
+
+
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+class tnlFastSweeping< 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();
+ bool init( const tnlParameterContainer& parameters );
+
+ bool initGrid();
+ bool run();
+ void updateValue(const Index i, const Index j);
+ Real fabsMin(const Real x, const Real y);
+
+
+protected:
+
+ MeshType Mesh;
+
+ DofVectorType dofVector;
+
+ RealType h;
+
+
+};
+
+
+#include "tnlFastSweeping2D_impl.h"
+
+#endif /* TNLFASTSWEEPING_H_ */
diff --git a/examples/fast-sweeping/tnlFastSweeping2D_impl.h b/examples/fast-sweeping/tnlFastSweeping2D_impl.h
new file mode 100644
index 0000000000..0abf197618
--- /dev/null
+++ b/examples/fast-sweeping/tnlFastSweeping2D_impl.h
@@ -0,0 +1,324 @@
+/***************************************************************************
+ tnlFastSweeping_impl.h - description
+ -------------------
+ begin : Oct 15 , 2015
+ copyright : (C) 2015 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 TNLFASTSWEEPING2D_IMPL_H_
+#define TNLFASTSWEEPING2D_IMPL_H_
+
+#include "tnlFastSweeping.h"
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+tnlString tnlFastSweeping< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index > :: getType()
+{
+ return tnlString( "tnlFastSweeping< " ) +
+ MeshType::getType() + ", " +
+ ::getType< Real >() + ", " +
+ ::getType< Index >() + " >";
+}
+
+
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+bool tnlFastSweeping< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index > :: init( const tnlParameterContainer& parameters )
+{
+ const tnlString& meshFile = parameters.getParameter< tnlString >( "mesh" );
+
+ if( ! Mesh.load( meshFile ) )
+ {
+ cerr << "I am not able to load the mesh from the file " << meshFile << "." << endl;
+ return false;
+ }
+
+
+ const tnlString& initialCondition = parameters.getParameter <tnlString>("initial-condition");
+ if( ! dofVector.load( initialCondition ) )
+ {
+ cerr << "I am not able to load the initial condition from the file " << meshFile << "." << endl;
+ return false;
+ }
+
+ h = Mesh.getHx();
+
+ return initGrid();
+}
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+bool tnlFastSweeping< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index > :: initGrid()
+{
+
+ Real tmp = 0.0;
+ for(Index i = 1; i < Mesh.getDimensions().x()-1; i++)
+ {
+ for(Index j = 1; j < Mesh.getDimensions().y()-1; j++)
+ {
+ tmp = Sign(dofVector[Mesh.getCellIndex(CoordinatesType(i,j))]);
+
+ if(tmp == 0.0)
+ {}
+ else if(dofVector[Mesh.getCellIndex(CoordinatesType(i+1,j))]*tmp < 0.0 ||
+ dofVector[Mesh.getCellIndex(CoordinatesType(i-1,j))]*tmp < 0.0 ||
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j+1))]*tmp < 0.0 ||
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j-1))]*tmp < 0.0 )
+ {}
+ else
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j))] = tmp*INT_MAX;
+ }
+ }
+
+
+
+ for(int i = 1; i < Mesh.getDimensions().x()-1; i++)
+ {
+ Index j = 0;
+ tmp = Sign(dofVector[Mesh.getCellIndex(CoordinatesType(i,j))]);
+
+
+ if(tmp == 0.0)
+ {}
+ else if(dofVector[Mesh.getCellIndex(CoordinatesType(i+1,j))]*tmp < 0.0 ||
+ dofVector[Mesh.getCellIndex(CoordinatesType(i-1,j))]*tmp < 0.0 ||
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j+1))]*tmp < 0.0 )
+ {}
+ else
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j))] = tmp*INT_MAX;
+ }
+
+ for(int i = 1; i < Mesh.getDimensions().x()-1; i++)
+ {
+ Index j = Mesh.getDimensions().y() - 1;
+ tmp = Sign(dofVector[Mesh.getCellIndex(CoordinatesType(i,j))]);
+
+
+ if(tmp == 0.0)
+ {}
+ else if(dofVector[Mesh.getCellIndex(CoordinatesType(i+1,j))]*tmp < 0.0 ||
+ dofVector[Mesh.getCellIndex(CoordinatesType(i-1,j))]*tmp < 0.0 ||
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j-1))]*tmp < 0.0 )
+ {}
+ else
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j))] = tmp*INT_MAX;
+ }
+
+ for(int j = 1; j < Mesh.getDimensions().y()-1; j++)
+ {
+ Index i = 0;
+ tmp = Sign(dofVector[Mesh.getCellIndex(CoordinatesType(i,j))]);
+
+
+ if(tmp == 0.0)
+ {}
+ else if(dofVector[Mesh.getCellIndex(CoordinatesType(i+1,j))]*tmp < 0.0 ||
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j+1))]*tmp < 0.0 ||
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j-1))]*tmp < 0.0 )
+ {}
+ else
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j))] = tmp*INT_MAX;
+ }
+
+ for(int j = 1; j < Mesh.getDimensions().y()-1; j++)
+ {
+ Index i = Mesh.getDimensions().x() - 1;
+ tmp = Sign(dofVector[Mesh.getCellIndex(CoordinatesType(i,j))]);
+
+
+ if(tmp == 0.0)
+ {}
+ else if(dofVector[Mesh.getCellIndex(CoordinatesType(i-1,j))]*tmp < 0.0 ||
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j+1))]*tmp < 0.0 ||
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j-1))]*tmp < 0.0 )
+ {}
+ else
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j))] = tmp*INT_MAX;
+ }
+
+
+ Index i = Mesh.getDimensions().x() - 1;
+ Index j = Mesh.getDimensions().y() - 1;
+
+ tmp = Sign(dofVector[Mesh.getCellIndex(CoordinatesType(i,j))]);
+ if(dofVector[Mesh.getCellIndex(CoordinatesType(i-1,j))]*tmp > 0.0 &&
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j-1))]*tmp > 0.0)
+
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j))] = tmp*INT_MAX;
+
+
+
+ j = 0;
+ tmp = Sign(dofVector[Mesh.getCellIndex(CoordinatesType(i,j))]);
+ if(dofVector[Mesh.getCellIndex(CoordinatesType(i-1,j))]*tmp > 0.0 &&
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j+1))]*tmp > 0.0)
+
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j))] = tmp*INT_MAX;
+
+
+
+ i = 0;
+ j = Mesh.getDimensions().y() -1;
+ tmp = Sign(dofVector[Mesh.getCellIndex(CoordinatesType(i,j))]);
+ if(dofVector[Mesh.getCellIndex(CoordinatesType(i+1,j))]*tmp > 0.0 &&
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j-1))]*tmp > 0.0)
+
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j))] = tmp*INT_MAX;
+
+
+
+ j = 0;
+ tmp = Sign(dofVector[Mesh.getCellIndex(CoordinatesType(i,j))]);
+ if(dofVector[Mesh.getCellIndex(CoordinatesType(i+1,j))]*tmp > 0.0 &&
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j+1))]*tmp > 0.0)
+
+ dofVector[Mesh.getCellIndex(CoordinatesType(i,j))] = tmp*INT_MAX;
+
+
+ dofVector.save("u-00000.tnl");
+
+ return true;
+}
+
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+bool tnlFastSweeping< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index > :: run()
+{
+
+ for(Index i = 0; i < Mesh.getDimensions().x(); i++)
+ {
+ for(Index j = 0; j < Mesh.getDimensions().y(); j++)
+ {
+ updateValue(i,j);
+ }
+ }
+
+/*---------------------------------------------------------------------------------------------------------------------------*/
+
+ for(Index i = Mesh.getDimensions().x() - 1; i > -1; i--)
+ {
+ for(Index j = 0; j < Mesh.getDimensions().y(); j++)
+ {
+ updateValue(i,j);
+ }
+ }
+
+/*---------------------------------------------------------------------------------------------------------------------------*/
+
+ for(Index i = Mesh.getDimensions().x() - 1; i > -1; i--)
+ {
+ for(Index j = Mesh.getDimensions().y() - 1; j > -1; j--)
+ {
+ updateValue(i,j);
+ }
+ }
+
+/*---------------------------------------------------------------------------------------------------------------------------*/
+ for(Index i = 0; i < Mesh.getDimensions().x(); i++)
+ {
+ for(Index j = Mesh.getDimensions().y() - 1; j > -1; j--)
+ {
+ updateValue(i,j);
+ }
+ }
+
+/*---------------------------------------------------------------------------------------------------------------------------*/
+
+
+ dofVector.save("u-00001.tnl");
+
+ return true;
+}
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+void tnlFastSweeping< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index > :: updateValue( Index i, Index j)
+{
+ Index index = Mesh.getCellIndex(CoordinatesType(i,j));
+ Real value = dofVector[index];
+ Real a,b, tmp;
+
+ if( i == 0 )
+ a = dofVector[Mesh.template getCellNextToCell<1,0>(index)];
+ else if( i == Mesh.getDimensions().x() - 1 )
+ a = dofVector[Mesh.template getCellNextToCell<-1,0>(index)];
+ else
+ {
+ a = fabsMin( dofVector[Mesh.template getCellNextToCell<-1,0>(index)],
+ dofVector[Mesh.template getCellNextToCell<1,0>(index)] );
+ }
+
+ if( j == 0 )
+ b = dofVector[Mesh.template getCellNextToCell<0,1>(index)];
+ else if( j == Mesh.getDimensions().y() - 1 )
+ b = dofVector[Mesh.template getCellNextToCell<0,-1>(index)];
+ else
+ {
+ b = fabsMin( dofVector[Mesh.template getCellNextToCell<0,-1>(index)],
+ dofVector[Mesh.template getCellNextToCell<0,1>(index)] );
+ }
+
+
+ if(fabs(a-b) >= h)
+ tmp = fabsMin(a,b) + Sign(value)*h;
+ else
+ tmp = 0.5 * (a + b + Sign(value)*sqrt(2.0 * h * h - (a - b) * (a - b) ) );
+
+
+ dofVector[index] = fabsMin(value, tmp);
+}
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+Real tnlFastSweeping< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index > :: fabsMin( Real x, Real y)
+{
+ Real fx = fabs(x);
+ Real fy = fabs(y);
+
+ Real tmpMin = Min(fx,fy);
+
+ if(tmpMin == fx)
+ return x;
+ else
+ return y;
+
+
+}
+
+
+
+
+#endif /* TNLFASTSWEEPING_IMPL_H_ */
diff --git a/examples/hamilton-jacobi-parallel/main.h b/examples/hamilton-jacobi-parallel/main.h
index 8fe0efdc6c..b0edef6915 100644
--- a/examples/hamilton-jacobi-parallel/main.h
+++ b/examples/hamilton-jacobi-parallel/main.h
@@ -22,6 +22,7 @@
#include <mesh/tnlGrid.h>
#include <core/tnlDevice.h>
#include <time.h>
+#include <ctime>
typedef MainBuildConfig BuildConfig;
@@ -30,6 +31,7 @@ int main( int argc, char* argv[] )
time_t start;
time_t stop;
time(&start);
+ std::clock_t start2= std::clock();
tnlParameterContainer parameters;
tnlConfigDescription configDescription;
parallelEikonalConfig< BuildConfig >::configSetup( configDescription );
@@ -84,7 +86,7 @@ int main( int argc, char* argv[] )
time(&stop);
cout << endl;
- cout << "Running time was: " << difftime(stop,start) << endl;
+ cout << "Running time was: " << difftime(stop,start) << " .... " << (std::clock() - start2) / (double)(CLOCKS_PER_SEC) << endl;
return EXIT_SUCCESS;
}
diff --git a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h
index 59cc75a286..a20abadec1 100644
--- a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h
+++ b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h
@@ -25,7 +25,7 @@ template< typename SchemeHost, typename SchemeDevice, typename Device>
tnlParallelEikonalSolver<SchemeHost, SchemeDevice, Device, double, int>::tnlParallelEikonalSolver()
{
cout << "a" << endl;
- this->device = tnlHostDevice;
+ this->device = tnlCudaDevice; /////////////// tnlCuda Device --- vypocet na GPU, tnlHostDevice --- vypocet na CPU
#ifdef HAVE_CUDA
if(this->device == tnlCudaDevice)
@@ -85,10 +85,10 @@ bool tnlParallelEikonalSolver<SchemeHost, SchemeDevice, Device, double, int>::in
stretchGrid();
this->stopTime /= (double)(this->gridCols);
- this->stopTime *= (1.0+1.0/((double)(this->n) - 1.0));
- cout << "Setting stopping time to " << this->stopTime << endl;
- this->stopTime = 1.5*((double)(this->n))*parameters.getParameter <double>("stop-time")*this->mesh.getHx();
+ this->stopTime *= (1.0+1.0/((double)(this->n) - 2.0));
cout << "Setting stopping time to " << this->stopTime << endl;
+ //this->stopTime = 1.5*((double)(this->n))*parameters.getParameter <double>("stop-time")*this->mesh.getHx();
+ //cout << "Setting stopping time to " << this->stopTime << endl;
cout << "Initializating scheme..." << endl;
if(!this->schemeHost.init(parameters))
@@ -251,6 +251,9 @@ void tnlParallelEikonalSolver<SchemeHost, SchemeDevice, Device, double, int>::ru
{
//cout << "subMesh: " << i << ", BC: " << getBoundaryCondition(i) << endl;
+ if(getSubgridValue(i) == currentStep+4)
+ {
+
if(getBoundaryCondition(i) & 1)
{
insertSubgrid( runSubgrid(1, getSubgrid(i),i), i);
@@ -271,27 +274,27 @@ void tnlParallelEikonalSolver<SchemeHost, SchemeDevice, Device, double, int>::ru
insertSubgrid( runSubgrid(8, getSubgrid(i),i), i);
this->calculationsCount[i]++;
}
+ }
-
- if( ((getBoundaryCondition(i) & 2) )//|| (getBoundaryCondition(i) & 1))
+ if( ((getBoundaryCondition(i) & 2) )|| (getBoundaryCondition(i) & 1)//)
/* &&(!(getBoundaryCondition(i) & 5) && !(getBoundaryCondition(i) & 10)) */)
{
//cout << "3 @ " << getBoundaryCondition(i) << endl;
insertSubgrid( runSubgrid(3, getSubgrid(i),i), i);
}
- if( ((getBoundaryCondition(i) & 4) )//|| (getBoundaryCondition(i) & 1))
+ if( ((getBoundaryCondition(i) & 4) )|| (getBoundaryCondition(i) & 1)//)
/* &&(!(getBoundaryCondition(i) & 3) && !(getBoundaryCondition(i) & 12)) */)
{
//cout << "5 @ " << getBoundaryCondition(i) << endl;
insertSubgrid( runSubgrid(5, getSubgrid(i),i), i);
}
- if( ((getBoundaryCondition(i) & 2) )//|| (getBoundaryCondition(i) & 8))
+ if( ((getBoundaryCondition(i) & 2) )|| (getBoundaryCondition(i) & 8)//)
/* &&(!(getBoundaryCondition(i) & 12) && !(getBoundaryCondition(i) & 3))*/ )
{
//cout << "10 @ " << getBoundaryCondition(i) << endl;
insertSubgrid( runSubgrid(10, getSubgrid(i),i), i);
}
- if( ((getBoundaryCondition(i) & 4) )//|| (getBoundaryCondition(i) & 8))
+ if( ((getBoundaryCondition(i) & 4) )|| (getBoundaryCondition(i) & 8)//)
/*&&(!(getBoundaryCondition(i) & 10) && !(getBoundaryCondition(i) & 5)) */)
{
//cout << "12 @ " << getBoundaryCondition(i) << endl;
@@ -977,13 +980,13 @@ tnlParallelEikonalSolver<SchemeHost, SchemeDevice, Device, double, int>::runSubg
if( time + currentTau > finalTime ) currentTau = finalTime - time;
- for( int i = 0; i < fu.getSize(); i ++ )
- {
- //cout << "Too big RHS! i = " << i << ", fu = " << fu[i] << ", u = " << u[i] << endl;
- if((u[i]+currentTau * fu[ i ])*u[i] < 0.0 && fu[i] != 0.0 && u[i] != 0.0 )
- currentTau = fabs(u[i]/(2.0*fu[i]));
-
- }
+// for( int i = 0; i < fu.getSize(); i ++ )
+// {
+// //cout << "Too big RHS! i = " << i << ", fu = " << fu[i] << ", u = " << u[i] << endl;
+// if((u[i]+currentTau * fu[ i ])*u[i] < 0.0 && fu[i] != 0.0 && u[i] != 0.0 )
+// currentTau = fabs(u[i]/(2.0*fu[i]));
+//
+// }
for( int i = 0; i < fu.getSize(); i ++ )
@@ -1019,11 +1022,11 @@ template< typename SchemeHost, typename SchemeDevice, typename Device>
__device__
void tnlParallelEikonalSolver<SchemeHost, SchemeDevice, Device, double, int>::getSubgridCUDA( const int i ,tnlParallelEikonalSolver<SchemeHost, SchemeDevice, Device, double, int >* caller, double* a)
{
- int j = threadIdx.x + threadIdx.y * blockDim.x;
- int th = (i / caller->gridCols) * caller->n*caller->n*caller->gridCols
- + (i % caller->gridCols) * caller->n
- + (j/caller->n) * caller->n*caller->gridCols
- + (j % caller->n);
+ //int j = threadIdx.x + threadIdx.y * blockDim.x;
+ int th = (blockIdx.y) * caller->n*caller->n*caller->gridCols
+ + (blockIdx.x) * caller->n
+ + threadIdx.y * caller->n*caller->gridCols
+ + threadIdx.x;
//printf("i= %d,j= %d,th= %d\n",i,j,th);
*a = caller->work_u_cuda[th];
//printf("Hi %f \n", *a);
@@ -1036,10 +1039,10 @@ __device__
void tnlParallelEikonalSolver<SchemeHost, SchemeDevice, Device, double, int>::updateSubgridCUDA( const int i ,tnlParallelEikonalSolver<SchemeHost, SchemeDevice, Device, double, int >* caller, double* a)
{
int j = threadIdx.x + threadIdx.y * blockDim.x;
- int index = (i / caller->gridCols) * caller->n*caller->n*caller->gridCols
- + (i % caller->gridCols) * caller->n
- + (j/caller->n) * caller->n*caller->gridCols
- + (j % caller->n);
+ int index = (blockIdx.y) * caller->n*caller->n*caller->gridCols
+ + (blockIdx.x) * caller->n
+ + threadIdx.y * caller->n*caller->gridCols
+ + threadIdx.x;
if( (fabs(caller->work_u_cuda[index]) > fabs(*a)) || (caller->unusedCell_cuda[index] == 1) )
{
@@ -1058,13 +1061,13 @@ void tnlParallelEikonalSolver<SchemeHost, SchemeDevice, Device, double, int>::in
{
- int j = threadIdx.x + threadIdx.y * blockDim.x;
+// int j = threadIdx.x + threadIdx.y * blockDim.x;
//printf("j = %d, u = %f\n", j,u);
- 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);
+ int index = (blockIdx.y)*this->n*this->n*this->gridCols
+ + (blockIdx.x)*this->n
+ + threadIdx.y*this->n*this->gridCols
+ + threadIdx.x;
//printf("i= %d,j= %d,index= %d\n",i,j,index);
if( (fabs(this->work_u_cuda[index]) > fabs(u)) || (this->unusedCell_cuda[index] == 1) )
@@ -1111,37 +1114,48 @@ void tnlParallelEikonalSolver<SchemeHost, SchemeDevice, Device, double, int>::ru
__syncthreads();
if(tmp !=1)
{
+// if(computeFU)
+// absVal[l]=0.0;
+// else
+// absVal[l] = fabs(u[l]);
+//
+// __syncthreads();
+//
+// if((blockDim.x == 16) && (l < 128)) absVal[l] = Max(absVal[l],absVal[l+128]);
+// __syncthreads();
+// if((blockDim.x == 16) && (l < 64)) absVal[l] = Max(absVal[l],absVal[l+64]);
+// __syncthreads();
+// if(l < 32) absVal[l] = Max(absVal[l],absVal[l+32]);
+// if(l < 16) absVal[l] = Max(absVal[l],absVal[l+16]);
+// if(l < 8) absVal[l] = Max(absVal[l],absVal[l+8]);
+// if(l < 4) absVal[l] = Max(absVal[l],absVal[l+4]);
+// if(l < 2) absVal[l] = Max(absVal[l],absVal[l+2]);
+// if(l < 1) value = Sign(u[0])*Max(absVal[l],absVal[l+1]);
+// __syncthreads();
+//
+// if(computeFU)
+// u[l] = value;
if(computeFU)
- absVal[l]=0.0;
- else
- absVal[l] = fabs(u[l]);
-
- __syncthreads();
-
- if((blockDim.x == 16) && (l < 128)) absVal[l] = Max(absVal[l],absVal[l+128]);
- __syncthreads();
- if((blockDim.x == 16) && (l < 64)) absVal[l] = Max(absVal[l],absVal[l+64]);
- __syncthreads();
- if(l < 32) absVal[l] = Max(absVal[l],absVal[l+32]);
- if(l < 16) absVal[l] = Max(absVal[l],absVal[l+16]);
- if(l < 8) absVal[l] = Max(absVal[l],absVal[l+8]);
- if(l < 4) absVal[l] = Max(absVal[l],absVal[l+4]);
- if(l < 2) absVal[l] = Max(absVal[l],absVal[l+2]);
- if(l < 1) value = Sign(u[0])*Max(absVal[l],absVal[l+1]);
- __syncthreads();
-
- if(computeFU)
- u[l] = value;
+ {
+ if(boundaryCondition == 4)
+ u[l] = u[threadIdx.y * blockDim.x] + Sign(u[0])*this->subMesh.getHx()*(threadIdx.x) ;//+ 2*Sign(u[0])*this->subMesh.getHx()*(threadIdx.x+this->n);
+ else if(boundaryCondition == 2)
+ u[l] = u[threadIdx.y * blockDim.x + blockDim.x - 1] + Sign(u[0])*this->subMesh.getHx()*(this->n - 1 - threadIdx.x);//+ 2*Sign(u[0])*this->subMesh.getHx()*(blockDim.x - threadIdx.x - 1+this->n);
+ else if(boundaryCondition == 8)
+ u[l] = u[threadIdx.x] + Sign(u[0])*this->subMesh.getHx()*(threadIdx.y) ;//+ 2*Sign(u[0])*this->subMesh.getHx()*(threadIdx.y+this->n);
+ else if(boundaryCondition == 1)
+ u[l] = u[(blockDim.y - 1)* blockDim.x + threadIdx.x] + Sign(u[0])*this->subMesh.getHx()*(this->n - 1 - threadIdx.y) ;//+ Sign(u[0])*this->subMesh.getHx()*(blockDim.y - threadIdx.y - 1 +this->n);
+ }
}
double time = 0.0;
__shared__ double currentTau;
double cfl = this->cflCondition;
double fu = 0.0;
- if(threadIdx.x * threadIdx.y == 0)
- {
- currentTau = this->tau0;
- }
+// if(threadIdx.x * threadIdx.y == 0)
+// {
+// currentTau = this->tau0;
+// }
double finalTime = this->stopTime;
__syncthreads();
if( time + currentTau > finalTime ) currentTau = finalTime - time;
@@ -1149,18 +1163,18 @@ void tnlParallelEikonalSolver<SchemeHost, SchemeDevice, Device, double, int>::ru
while( time < finalTime )
{
-
if(computeFU)
fu = schemeHost.getValueDev( this->subMesh, l, tnlStaticVector<2,int>(i,j)/*this->subMesh.getCellCoordinates(l)*/, u, time, boundaryCondition);
sharedTau[l]=cfl/fabs(fu);
if(l == 0)
+ {
if(sharedTau[0] > 1.0 * this->subMesh.getHx()) sharedTau[0] = 1.0 * this->subMesh.getHx();
-
- if(l == blockDim.x*blockDim.y - 1)
+ }
+ else if(l == blockDim.x*blockDim.y - 1)
if( time + sharedTau[l] > finalTime ) sharedTau[l] = finalTime - time;
- __syncthreads();
+
if((blockDim.x == 16) && (l < 128)) sharedTau[l] = Min(sharedTau[l],sharedTau[l+128]);
@@ -1177,7 +1191,6 @@ void tnlParallelEikonalSolver<SchemeHost, SchemeDevice, Device, double, int>::ru
u[l] += currentTau * fu;
time += currentTau;
- __syncthreads();
}
@@ -1626,98 +1639,222 @@ void /*tnlParallelEikonalSolver<SchemeHost, SchemeDevice, Device, double, int>::
extern __shared__ double u[];
int i = blockIdx.y * gridDim.x + blockIdx.x;
int l = threadIdx.y * blockDim.x + threadIdx.x;
+ int bound = caller->getBoundaryConditionCUDA(i);
- if(caller->getSubgridValueCUDA(i) != INT_MAX && caller->getSubgridValueCUDA(i) >= 0)
+ if(caller->getSubgridValueCUDA(i) != INT_MAX && bound != 0 && caller->getSubgridValueCUDA(i) > 0)
{
caller->getSubgridCUDA(i,caller, &u[l]);
- int bound = caller->getBoundaryConditionCUDA(i);
+
//if(l == 0)
//printf("i = %d, bound = %d\n",i,caller->getSubgridValueCUDA(i));
if(caller->getSubgridValueCUDA(i) == caller->currentStep+4)
{
- if(bound & 1)
- {
- caller->runSubgridCUDA(1,u,i);
- __syncthreads();
- //caller->insertSubgridCUDA(u[l],i);
- //__syncthreads();
- //caller->getSubgridCUDA(i,caller, &u[l]);
- caller->updateSubgridCUDA(i,caller, &u[l]);
- __syncthreads();
- }
- if(bound & 2 )
- {
- caller->runSubgridCUDA(2,u,i);
- __syncthreads();
- //caller->insertSubgridCUDA(u[l],i);
- //__syncthreads();
- //caller->getSubgridCUDA(i,caller, &u[l]);
- caller->updateSubgridCUDA(i,caller, &u[l]);
- __syncthreads();
- }
- if(bound & 4)
- {
- caller->runSubgridCUDA(4,u,i);
- __syncthreads();
- //caller->insertSubgridCUDA(u[l],i);
- //__syncthreads();
- //caller->getSubgridCUDA(i,caller, &u[l]);
- caller->updateSubgridCUDA(i,caller, &u[l]);
- __syncthreads();
- }
- if(bound & 8)
- {
- caller->runSubgridCUDA(8,u,i);
- __syncthreads();
- //caller->insertSubgridCUDA(u[l],i);
- //__syncthreads();
- //caller->getSubgridCUDA(i,caller, &u[l]);
- caller->updateSubgridCUDA(i,caller, &u[l]);
- __syncthreads();
- }
- }
+ if(bound & 1)
+ {
+ caller->runSubgridCUDA(1,u,i);
+ //__syncthreads();
+ //caller->insertSubgridCUDA(u[l],i);
+ //__syncthreads();
+ //caller->getSubgridCUDA(i,caller, &u[l]);
+ caller->updateSubgridCUDA(i,caller, &u[l]);
+ __syncthreads();
+ }
+ if(bound & 2 )
+ {
+ caller->runSubgridCUDA(2,u,i);
+ //__syncthreads();
+ //caller->insertSubgridCUDA(u[l],i);
+ //__syncthreads();
+ //caller->getSubgridCUDA(i,caller, &u[l]);
+ caller->updateSubgridCUDA(i,caller, &u[l]);
+ __syncthreads();
+ }
+ if(bound & 4)
+ {
+ caller->runSubgridCUDA(4,u,i);
+ //__syncthreads();
+ //caller->insertSubgridCUDA(u[l],i);
+ //__syncthreads();
+ //caller->getSubgridCUDA(i,caller, &u[l]);
+ caller->updateSubgridCUDA(i,caller, &u[l]);
+ __syncthreads();
+ }
+ if(bound & 8)
+ {
+ caller->runSubgridCUDA(8,u,i);
+ //__syncthreads();
+ //caller->insertSubgridCUDA(u[l],i);
+ //__syncthreads();
+ //caller->getSubgridCUDA(i,caller, &u[l]);
+ caller->updateSubgridCUDA(i,caller, &u[l]);
+ __syncthreads();
+ }
+
+
+
+
+
+ if( ((bound & 3 )))
+ {
+ caller->runSubgridCUDA(3,u,i);
+ //__syncthreads();
+ //caller->insertSubgridCUDA(u[l],i);
+ //__syncthreads();
+ //caller->getSubgridCUDA(i,caller, &u[l]);
+ caller->updateSubgridCUDA(i,caller, &u[l]);
+ __syncthreads();
+ }
+ if( ((bound & 5 )))
+ {
+ caller->runSubgridCUDA(5,u,i);
+ //__syncthreads();
+ //caller->insertSubgridCUDA(u[l],i);
+ //__syncthreads();
+ //caller->getSubgridCUDA(i,caller, &u[l]);
+ caller->updateSubgridCUDA(i,caller, &u[l]);
+ __syncthreads();
+ }
+ if( ((bound & 10 )))
+ {
+ caller->runSubgridCUDA(10,u,i);
+ //__syncthreads();
+ //caller->insertSubgridCUDA(u[l],i);
+ //__syncthreads();
+ //caller->getSubgridCUDA(i,caller, &u[l]);
+ caller->updateSubgridCUDA(i,caller, &u[l]);
+ __syncthreads();
+ }
+ if( (bound & 12 ))
+ {
+ caller->runSubgridCUDA(12,u,i);
+ //__syncthreads();
+ //caller->insertSubgridCUDA(u[l],i);
+ //__syncthreads();
+ //caller->getSubgridCUDA(i,caller, &u[l]);
+ caller->updateSubgridCUDA(i,caller, &u[l]);
+ __syncthreads();
+ }
+
+
- if( ((bound & 2) || (bound & 1) ))
- {
- caller->runSubgridCUDA(3,u,i);
- __syncthreads();
- //caller->insertSubgridCUDA(u[l],i);
- //__syncthreads();
- //caller->getSubgridCUDA(i,caller, &u[l]);
- caller->updateSubgridCUDA(i,caller, &u[l]);
- __syncthreads();
- }
- if( ((bound & 4) || (bound & 1) ))
- {
- caller->runSubgridCUDA(5,u,i);
- __syncthreads();
- //caller->insertSubgridCUDA(u[l],i);
- //__syncthreads();
- //caller->getSubgridCUDA(i,caller, &u[l]);
- caller->updateSubgridCUDA(i,caller, &u[l]);
- __syncthreads();
- }
- if( ((bound & 2) || (bound & 8) ))
- {
- caller->runSubgridCUDA(10,u,i);
- __syncthreads();
- //caller->insertSubgridCUDA(u[l],i);
- //__syncthreads();
- //caller->getSubgridCUDA(i,caller, &u[l]);
- caller->updateSubgridCUDA(i,caller, &u[l]);
- __syncthreads();
}
- if( (bound & 4) || (bound & 8) )
+
+
+ else
{
- caller->runSubgridCUDA(12,u,i);
- __syncthreads();
- //caller->insertSubgridCUDA(u[l],i);
- //__syncthreads();
- //caller->getSubgridCUDA(i,caller, &u[l]);
- caller->updateSubgridCUDA(i,caller, &u[l]);
- __syncthreads();
+
+
+
+
+
+
+
+
+
+ if( ((bound == 2)))
+ {
+ caller->runSubgridCUDA(2,u,i);
+ //__syncthreads();
+ //caller->insertSubgridCUDA(u[l],i);
+ //__syncthreads();
+ //caller->getSubgridCUDA(i,caller, &u[l]);
+ caller->updateSubgridCUDA(i,caller, &u[l]);
+ __syncthreads();
+ }
+ if( ((bound == 1) ))
+ {
+ caller->runSubgridCUDA(1,u,i);
+ //__syncthreads();
+ //caller->insertSubgridCUDA(u[l],i);
+ //__syncthreads();
+ //caller->getSubgridCUDA(i,caller, &u[l]);
+ caller->updateSubgridCUDA(i,caller, &u[l]);
+ __syncthreads();
+ }
+ if( ((bound == 8) ))
+ {
+ caller->runSubgridCUDA(8,u,i);
+ //__syncthreads();
+ //caller->insertSubgridCUDA(u[l],i);
+ //__syncthreads();
+ //caller->getSubgridCUDA(i,caller, &u[l]);
+ caller->updateSubgridCUDA(i,caller, &u[l]);
+ __syncthreads();
+ }
+ if( (bound == 4))
+ {
+ caller->runSubgridCUDA(4,u,i);
+ //__syncthreads();
+ //caller->insertSubgridCUDA(u[l],i);
+ //__syncthreads();
+ //caller->getSubgridCUDA(i,caller, &u[l]);
+ caller->updateSubgridCUDA(i,caller, &u[l]);
+ __syncthreads();
+ }
+
+
+
+
+
+
+
+
+
+
+ if( ((bound & 3) ))
+ {
+ caller->runSubgridCUDA(3,u,i);
+ //__syncthreads();
+ //caller->insertSubgridCUDA(u[l],i);
+ //__syncthreads();
+ //caller->getSubgridCUDA(i,caller, &u[l]);
+ caller->updateSubgridCUDA(i,caller, &u[l]);
+ __syncthreads();
+ }
+ if( ((bound & 5) ))
+ {
+ caller->runSubgridCUDA(5,u,i);
+ //__syncthreads();
+ //caller->insertSubgridCUDA(u[l],i);
+ //__syncthreads();
+ //caller->getSubgridCUDA(i,caller, &u[l]);
+ caller->updateSubgridCUDA(i,caller, &u[l]);
+ __syncthreads();
+ }
+ if( ((bound & 10) ))
+ {
+ caller->runSubgridCUDA(10,u,i);
+ //__syncthreads();
+ //caller->insertSubgridCUDA(u[l],i);
+ //__syncthreads();
+ //caller->getSubgridCUDA(i,caller, &u[l]);
+ caller->updateSubgridCUDA(i,caller, &u[l]);
+ __syncthreads();
+ }
+ if( (bound & 12) )
+ {
+ caller->runSubgridCUDA(12,u,i);
+ //__syncthreads();
+ //caller->insertSubgridCUDA(u[l],i);
+ //__syncthreads();
+ //caller->getSubgridCUDA(i,caller, &u[l]);
+ caller->updateSubgridCUDA(i,caller, &u[l]);
+ __syncthreads();
+ }
+
+
+
+
+
+
+
+
+
+
+
+
}
/*if( bound )
{
@@ -1730,15 +1867,17 @@ void /*tnlParallelEikonalSolver<SchemeHost, SchemeDevice, Device, double, int>::
__syncthreads();
}*/
-
if(l==0)
{
caller->setBoundaryConditionCUDA(i, 0);
caller->setSubgridValueCUDA(i, caller->getSubgridValueCUDA(i) - 1 );
}
+
}
+
+
}
#endif /*HAVE_CUDA*/
diff --git a/src/operators/godunov-eikonal/parallelGodunovEikonal.h b/src/operators/godunov-eikonal/parallelGodunovEikonal.h
index aa3e34e6f4..aff70555ea 100644
--- a/src/operators/godunov-eikonal/parallelGodunovEikonal.h
+++ b/src/operators/godunov-eikonal/parallelGodunovEikonal.h
@@ -223,6 +223,17 @@ public:
const Vector& u,
const RealType& time,
const IndexType boundaryCondition ) const;
+
+#ifdef HAVE_CUDA
+ __device__
+#endif
+ Real getValueDev( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const RealType* u,
+ const RealType& time,
+ const IndexType boundaryCondition) const;
+
#ifdef HAVE_CUDA
__device__ __host__
#endif
@@ -235,9 +246,9 @@ protected:
DofVectorType dofVector;
- RealType hx;
- RealType hy;
- RealType hz;
+ RealType hx, ihx;
+ RealType hy, ihy;
+ RealType hz, ihz;
RealType epsilon;
@@ -246,9 +257,9 @@ protected:
-#include <operators/godunov-eikonal/parallelGodunovEikonal1D_impl.h>
+//#include <operators/godunov-eikonal/parallelGodunovEikonal1D_impl.h>
#include <operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h>
-#include <operators/godunov-eikonal/parallelGodunovEikonal3D_impl.h>
+//#include <operators/godunov-eikonal/parallelGodunovEikonal3D_impl.h>
#endif /* GODUNOVEIKONAL_H_ */
diff --git a/src/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h b/src/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h
index 965d8d22d0..1b4906b185 100644
--- a/src/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h
+++ b/src/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h
@@ -65,7 +65,7 @@ Real parallelGodunovEikonalScheme< tnlGrid< 2,MeshReal, Device, MeshIndex >, Rea
if (x < -eps)
return (-1.0);
- if ( eps == 0.0)
+ if ( x == 0.0)
return 0.0;
return sin(/*(M_PI*x)/(2.0*eps) */(M_PI/2.0)*(x/eps));
@@ -87,7 +87,7 @@ bool parallelGodunovEikonalScheme< tnlGrid< 2,MeshReal, Device, MeshIndex >, Rea
{
const tnlString& meshFile = parameters.getParameter< tnlString >( "mesh" );
- MeshType originalMesh;
+ //MeshType originalMesh;
if( ! originalMesh.load( meshFile ) )
{
//cerr << "I am not able to load the mesh from the file " << meshFile << "." << endl;
@@ -102,9 +102,9 @@ bool parallelGodunovEikonalScheme< tnlGrid< 2,MeshReal, Device, MeshIndex >, Rea
hy = originalMesh.getHy();
ihy = 1.0/hy;
- epsilon = parameters. getParameter< double >( "epsilon" );
+ this->epsilon = parameters. getParameter< double >( "epsilon" );
- epsilon *=sqrt( hx*hx + hy*hy );
+ this->epsilon *=sqrt( hx*hx + hy*hy );
// dofVector. setSize( this->mesh.getDofs() );
@@ -123,7 +123,7 @@ template< typename MeshReal,
#endif
tnlString parallelGodunovEikonalScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index > :: getType()
{
- return tnlString( "tnlLinearDiffusion< " ) +
+ return tnlString( "parallelGodunovEikonalScheme< " ) +
MeshType::getType() + ", " +
::getType< Real >() + ", " +
::getType< Index >() + " >";
@@ -162,128 +162,218 @@ Real parallelGodunovEikonalScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Re
return 0.0;
}
- RealType acc = hx*hy*hx*hy;
- RealType nabla, xb, xf, yb, yf, signui;
+ //-----------------------------------
- signui = sign(u[cellIndex],epsilon);
+ RealType signui;
+ signui = sign(u[cellIndex],this->epsilon);
+
+
+#ifdef HAVE_CUDA
+ //printf("%d : %d ;;;; %d : %d , %f \n",threadIdx.x, mesh.getDimensions().x() , threadIdx.y,mesh.getDimensions().y(), epsilon );
+#endif
+ RealType xb = u[cellIndex];
+ RealType xf = -u[cellIndex];
+ RealType yb = u[cellIndex];
+ RealType yf = -u[cellIndex];
+ RealType a,b,c;
- if(fabs(u[cellIndex]) < acc) return 0.0;
+
+ if(coordinates.x() == mesh.getDimensions().x() - 1)
+ xf += u[mesh.template getCellNextToCell<-1,0>( cellIndex )];
+ else
+ xf += u[mesh.template getCellNextToCell<1,0>( cellIndex )];
+
+ if(coordinates.x() == 0)
+ xb -= u[mesh.template getCellNextToCell<1,0>( cellIndex )];
+ else
+ xb -= u[mesh.template getCellNextToCell<-1,0>( cellIndex )];
+
+ if(coordinates.y() == mesh.getDimensions().y() - 1)
+ yf += u[mesh.template getCellNextToCell<0,-1>( cellIndex )];
+ else
+ yf += u[mesh.template getCellNextToCell<0,1>( cellIndex )];
+
+ if(coordinates.y() == 0)
+ yb -= u[mesh.template getCellNextToCell<0,1>( cellIndex )];
+ else
+ yb -= u[mesh.template getCellNextToCell<0,-1>( cellIndex )];
+
+
+ //xb *= ihx;
+ //xf *= ihx;
+ // yb *= ihy;
+ //yf *= ihy;
if(signui > 0.0)
{
- /**/ /* if(boundaryCondition & 2)
- xf = (u[mesh.getCellXSuccessor( cellIndex )] - u[cellIndex])/hx;
- else *//*if(boundaryCondition & 4)
- xf = 0.0;
- else /**/if(coordinates.x() == mesh.getDimensions().x() - 1)
- xf = negativePart((u[mesh.template getCellNextToCell<-1,0>( cellIndex )] - u[cellIndex])*ihx);
- else
- xf = negativePart((u[mesh.template getCellNextToCell<1,0>( cellIndex )] - u[cellIndex])*ihx);
-
- /**/ /* if(boundaryCondition & 4)
- xb = (u[cellIndex] - u[mesh.getCellXPredecessor( cellIndex )])*ihx;
- else *//*if(boundaryCondition & 2)
- xb = 0.0;
- else /**/if(coordinates.x() == 0)
- xb = positivePart((u[cellIndex] - u[mesh.template getCellNextToCell<+1,0>( cellIndex )])*ihx);
- else
- xb = positivePart((u[cellIndex] - u[mesh.template getCellNextToCell<-1,0>( cellIndex )])*ihx);
-
- /**/ /* if(boundaryCondition & 1)
- yf = (u[mesh.getCellYSuccessor( cellIndex )] - u[cellIndex])*ihy;
- else *//*if(boundaryCondition & 8)
- yf = 0.0;
- else /**/if(coordinates.y() == mesh.getDimensions().y() - 1)
- yf = negativePart((u[mesh.template getCellNextToCell<0,-1>( cellIndex )] - u[cellIndex])*ihy);
- else
- yf = negativePart((u[mesh.template getCellNextToCell<0,1>( cellIndex )] - u[cellIndex])*ihy);
-
- /**/ /* if(boundaryCondition & 8)
- yb = (u[cellIndex] - u[mesh.getCellYPredecessor( cellIndex )])*ihy;
- else *//*if(boundaryCondition & 1)
- yb = 0.0;
- else /**/if(coordinates.y() == 0)
- yb = positivePart((u[cellIndex] - u[mesh.template getCellNextToCell<0,1>( cellIndex )])*ihy);
- else
- yb = positivePart((u[cellIndex] - u[mesh.template getCellNextToCell<0,-1>( cellIndex )])*ihy);
-
- 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 );
- if(fabs(1.0-nabla) < acc)
- return 0.0;
- return signui*(1.0 - nabla);
+ xf = negativePart(xf);
+
+ xb = positivePart(xb);
+
+ yf = negativePart(yf);
+
+ yb = positivePart(yb);
+
}
- else if (signui < 0.0)
+ else if(signui < 0.0)
{
- /**/ /* if(boundaryCondition & 2)
- xf = (u[mesh.getCellXSuccessor( cellIndex )] - u[cellIndex])*ihx;
- else*//* if(boundaryCondition & 4)
- xf = 0.0;
- else /**/if(coordinates.x() == mesh.getDimensions().x() - 1)
- xf = positivePart((u[mesh.template getCellNextToCell<-1,0>( cellIndex )] - u[cellIndex])*ihx);
- else
- xf = positivePart((u[mesh.template getCellNextToCell<1,0>( cellIndex )] - u[cellIndex])*ihx);
-
- /**/ /* if(boundaryCondition & 4)
- xb = (u[cellIndex] - u[mesh.getCellXPredecessor( cellIndex )])*ihx;
- else*//* if(boundaryCondition & 2)
- xb = 0.0;
- else /**/if(coordinates.x() == 0)
- xb = negativePart((u[cellIndex] - u[mesh.template getCellNextToCell<1,0>( cellIndex )])*ihx);
- else
- xb = negativePart((u[cellIndex] - u[mesh.template getCellNextToCell<-1,0>( cellIndex )])*ihx);
-
- /**/ /* if(boundaryCondition & 1)
- yf = (u[mesh.getCellYSuccessor( cellIndex )] - u[cellIndex])*ihy;
- else *//*if(boundaryCondition & 8)
- yf = 0.0;
- else /**/if(coordinates.y() == mesh.getDimensions().y() - 1)
- yf = positivePart((u[mesh.template getCellNextToCell<0,-1>( cellIndex )] - u[cellIndex])*ihy);
- else
- yf = positivePart((u[mesh.template getCellNextToCell<0,1>( cellIndex )] - u[cellIndex])*ihy);
-
- /**/ /* if(boundaryCondition & 8)
- yb = (u[cellIndex] - u[mesh.getCellYPredecessor( cellIndex )])*ihy;
- else*//* if(boundaryCondition & 1)
- yb = 0.0;
- else /**/if(coordinates.y() == 0)
- yb = negativePart((u[cellIndex] - u[mesh.template getCellNextToCell<0,1>( cellIndex )])*ihy);
- else
- yb = negativePart((u[cellIndex] - u[mesh.template getCellNextToCell<0,-1>( cellIndex )])*ihy);
-
-
- 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 );
-
- if(fabs(1.0-nabla) < acc)
- return 0.0;
- return signui*(1.0 - nabla);
+ xb = negativePart(xb);
+
+ xf = positivePart(xf);
+
+ yb = negativePart(yb);
+
+ yf = positivePart(yf);
}
+
+
+ if(xb - xf > 0.0)
+ a = xb;
else
- {
- return 0.0;
- }
+ a = xf;
+
+ if(yb - yf > 0.0)
+ b = yb;
+ else
+ b = yf;
+
+ c =(1.0 - sqrt(a*a+b*b)*ihx );
+
+ if(Sign(c) > 0.0 )
+ return Sign(u[cellIndex])*c;
+ else
+ return signui*c;
+
+ //--------------------------------------------------
+
+//
+//
+//
+// RealType acc = hx*hy*hx*hy;
+//
+// RealType nabla, xb, xf, yb, yf, signui;
+//
+// signui = sign(u[cellIndex],this->epsilon);
+//
+//
+// //if(fabs(u[cellIndex]) < acc) return 0.0;
+//
+// if(signui > 0.0)
+// {
+// /**/ /* if(boundaryCondition & 2)
+// xf = (u[mesh.getCellXSuccessor( cellIndex )] - u[cellIndex])/hx;
+// else *//*if(boundaryCondition & 4)
+// xf = 0.0;
+// else /**/if(coordinates.x() == mesh.getDimensions().x() - 1)
+// xf = negativePart((u[mesh.template getCellNextToCell<-1,0>( cellIndex )] - u[cellIndex])*ihx);
+// else
+// xf = negativePart((u[mesh.template getCellNextToCell<1,0>( cellIndex )] - u[cellIndex])*ihx);
+//
+// /**/ /* if(boundaryCondition & 4)
+// xb = (u[cellIndex] - u[mesh.getCellXPredecessor( cellIndex )])*ihx;
+// else *//*if(boundaryCondition & 2)
+// xb = 0.0;
+// else /**/if(coordinates.x() == 0)
+// xb = positivePart((u[cellIndex] - u[mesh.template getCellNextToCell<+1,0>( cellIndex )])*ihx);
+// else
+// xb = positivePart((u[cellIndex] - u[mesh.template getCellNextToCell<-1,0>( cellIndex )])*ihx);
+//
+// /**/ /* if(boundaryCondition & 1)
+// yf = (u[mesh.getCellYSuccessor( cellIndex )] - u[cellIndex])*ihy;
+// else *//*if(boundaryCondition & 8)
+// yf = 0.0;
+// else /**/if(coordinates.y() == mesh.getDimensions().y() - 1)
+// yf = negativePart((u[mesh.template getCellNextToCell<0,-1>( cellIndex )] - u[cellIndex])*ihy);
+// else
+// yf = negativePart((u[mesh.template getCellNextToCell<0,1>( cellIndex )] - u[cellIndex])*ihy);
+//
+// /**/ /* if(boundaryCondition & 8)
+// yb = (u[cellIndex] - u[mesh.getCellYPredecessor( cellIndex )])*ihy;
+// else *//*if(boundaryCondition & 1)
+// yb = 0.0;
+// else /**/if(coordinates.y() == 0)
+// yb = positivePart((u[cellIndex] - u[mesh.template getCellNextToCell<0,1>( cellIndex )])*ihy);
+// else
+// yb = positivePart((u[cellIndex] - u[mesh.template getCellNextToCell<0,-1>( cellIndex )])*ihy);
+//
+// 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 );
+// if(fabs(1.0-nabla) < acc)
+// return 0.0;
+// return signui*(1.0 - nabla);
+// }
+// else if (signui < 0.0)
+// {
+//
+// /**/ /* if(boundaryCondition & 2)
+// xf = (u[mesh.getCellXSuccessor( cellIndex )] - u[cellIndex])*ihx;
+// else*//* if(boundaryCondition & 4)
+// xf = 0.0;
+// else /**/if(coordinates.x() == mesh.getDimensions().x() - 1)
+// xf = positivePart((u[mesh.template getCellNextToCell<-1,0>( cellIndex )] - u[cellIndex])*ihx);
+// else
+// xf = positivePart((u[mesh.template getCellNextToCell<1,0>( cellIndex )] - u[cellIndex])*ihx);
+//
+// /**/ /* if(boundaryCondition & 4)
+// xb = (u[cellIndex] - u[mesh.getCellXPredecessor( cellIndex )])*ihx;
+// else*//* if(boundaryCondition & 2)
+// xb = 0.0;
+// else /**/if(coordinates.x() == 0)
+// xb = negativePart((u[cellIndex] - u[mesh.template getCellNextToCell<1,0>( cellIndex )])*ihx);
+// else
+// xb = negativePart((u[cellIndex] - u[mesh.template getCellNextToCell<-1,0>( cellIndex )])*ihx);
+//
+// /**/ /* if(boundaryCondition & 1)
+// yf = (u[mesh.getCellYSuccessor( cellIndex )] - u[cellIndex])*ihy;
+// else *//*if(boundaryCondition & 8)
+// yf = 0.0;
+// else /**/if(coordinates.y() == mesh.getDimensions().y() - 1)
+// yf = positivePart((u[mesh.template getCellNextToCell<0,-1>( cellIndex )] - u[cellIndex])*ihy);
+// else
+// yf = positivePart((u[mesh.template getCellNextToCell<0,1>( cellIndex )] - u[cellIndex])*ihy);
+//
+// /**/ /* if(boundaryCondition & 8)
+// yb = (u[cellIndex] - u[mesh.getCellYPredecessor( cellIndex )])*ihy;
+// else*//* if(boundaryCondition & 1)
+// yb = 0.0;
+// else /**/if(coordinates.y() == 0)
+// yb = negativePart((u[cellIndex] - u[mesh.template getCellNextToCell<0,1>( cellIndex )])*ihy);
+// else
+// yb = negativePart((u[cellIndex] - u[mesh.template getCellNextToCell<0,-1>( cellIndex )])*ihy);
+//
+//
+// 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 );
+//
+// if(fabs(1.0-nabla) < acc)
+// return 0.0;
+// return signui*(1.0 - nabla);
+// }
+// else
+// {
+// return 0.0;
+// }
}
@@ -312,31 +402,13 @@ Real parallelGodunovEikonalScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Re
const IndexType boundaryCondition) const
{
-/*
- if ( ((coordinates.x() == 0 && (boundaryCondition & 4)) or
- (coordinates.x() == mesh.getDimensions().x() - 1 && (boundaryCondition & 2)) or
- (coordinates.y() == 0 && (boundaryCondition & 8)) or
- (coordinates.y() == mesh.getDimensions().y() - 1 && (boundaryCondition & 1)))
- )
- {
- return 0.0;
- }
-
-*/
- //RealType acc = hx*hy*hx*hy;
-
- RealType signui;
- signui = sign(u[cellIndex],epsilon);
-
-#ifdef HAVE_CUDA
- //printf("%d : %d ;;;; %d : %d , %f \n",threadIdx.x, mesh.getDimensions().x() , threadIdx.y,mesh.getDimensions().y(), epsilon );
-#endif
+ RealType signui = sign(u[cellIndex],this->epsilon);
RealType xb = u[cellIndex];
RealType xf = -u[cellIndex];
RealType yb = u[cellIndex];
RealType yf = -u[cellIndex];
- RealType a,b;
+ RealType a,b,c;
if(coordinates.x() == mesh.getDimensions().x() - 1)
@@ -360,11 +432,6 @@ Real parallelGodunovEikonalScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Re
yb -= u[mesh.template getCellNextToCell<0,-1>( cellIndex )];
- //xb *= ihx;
- //xf *= ihx;
- // yb *= ihy;
- //yf *= ihy;
-
if(signui > 0.0)
{
xf = negativePart(xf);
@@ -399,7 +466,12 @@ Real parallelGodunovEikonalScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Re
else
b = yf;
- return signui*(1.0 - sqrt(a*a+b*b)*ihx );
+ c =(1.0 - sqrt(a*a+b*b)*ihx );
+
+ if(c > 0.0 )
+ return Sign(u[cellIndex])*c;
+ else
+ return signui*c;
}
diff --git a/src/operators/godunov-eikonal/parallelGodunovEikonal3D_impl.h b/src/operators/godunov-eikonal/parallelGodunovEikonal3D_impl.h
index 07d51e07a7..7b1346dca7 100644
--- a/src/operators/godunov-eikonal/parallelGodunovEikonal3D_impl.h
+++ b/src/operators/godunov-eikonal/parallelGodunovEikonal3D_impl.h
@@ -83,6 +83,9 @@ bool parallelGodunovEikonalScheme< tnlGrid< 3,MeshReal, Device, MeshIndex >, Rea
hx = originalMesh.getHx();
hy = originalMesh.getHy();
hz = originalMesh.getHz();
+ ihx = 1.0/hx;
+ ihy = 1.0/hy;
+ ihz = 1.0/hz;
epsilon = parameters. getParameter< double >( "epsilon" );
@@ -125,72 +128,267 @@ Real parallelGodunovEikonalScheme< tnlGrid< 3, MeshReal, Device, MeshIndex >, Re
const Real& time,
const IndexType boundaryCondition ) const
{
-/*
- RealType nabla, xb, xf, yb, yf, zb, zf, signui;
+ if ( ((coordinates.x() == 0 && (boundaryCondition & 4)) or
+ (coordinates.x() == mesh.getDimensions().x() - 1 && (boundaryCondition & 2)) or
+ (coordinates.y() == 0 && (boundaryCondition & 8)) or
+ (coordinates.y() == mesh.getDimensions().y() - 1 && (boundaryCondition & 1)) or
+ (coordinates.z() == 0 && (boundaryCondition & 32)) or
+ (coordinates.z() == mesh.getDimensions().y() - 1 && (boundaryCondition & 16)))
+
+ )
+ {
+ return 0.0;
+ }
+
+
+ //-----------------------------------
+
+ RealType signui;
+ signui = sign(u[cellIndex],this->epsilon);
+
+
+
- signui = sign(u[cellIndex],epsilon);
+ RealType xb = u[cellIndex];
+ RealType xf = -u[cellIndex];
+ RealType yb = u[cellIndex];
+ RealType yf = -u[cellIndex];
+ RealType zb = u[cellIndex];
+ RealType zf = -u[cellIndex];
+ RealType a,b,c,d;
+
+
+ if(coordinates.x() == mesh.getDimensions().x() - 1)
+ xf += u[mesh.template getCellNextToCell<-1,0,0>( cellIndex )];
+ else
+ xf += u[mesh.template getCellNextToCell<1,0,0>( cellIndex )];
+
+ if(coordinates.x() == 0)
+ xb -= u[mesh.template getCellNextToCell<1,0,0>( cellIndex )];
+ else
+ xb -= u[mesh.template getCellNextToCell<-1,0,0>( cellIndex )];
+
+ if(coordinates.y() == mesh.getDimensions().y() - 1)
+ yf += u[mesh.template getCellNextToCell<0,-1,0>( cellIndex )];
+ else
+ yf += u[mesh.template getCellNextToCell<0,1,0>( cellIndex )];
+
+ if(coordinates.y() == 0)
+ yb -= u[mesh.template getCellNextToCell<0,1,0>( cellIndex )];
+ else
+ yb -= u[mesh.template getCellNextToCell<0,-1,0>( cellIndex )];
+
+
+ if(coordinates.z() == mesh.getDimensions().z() - 1)
+ zf += u[mesh.template getCellNextToCell<0,0,-1>( cellIndex )];
+ else
+ zf += u[mesh.template getCellNextToCell<0,0,1>( cellIndex )];
+
+ if(coordinates.z() == 0)
+ zb -= u[mesh.template getCellNextToCell<0,0,1>( cellIndex )];
+ else
+ zb -= u[mesh.template getCellNextToCell<0,0,-1>( cellIndex )];
+
+
+ //xb *= ihx;
+ //xf *= ihx;
+ // yb *= ihy;
+ //yf *= ihy;
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);
+ xf = negativePart(xf);
+
+ xb = positivePart(xb);
+
+ yf = negativePart(yf);
+
+ yb = positivePart(yb);
+
+ zf = negativePart(zf);
+
+ zb = positivePart(zb);
+
}
- else if (signui < 0.0)
+ 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);
+
+ xb = negativePart(xb);
+
+ xf = positivePart(xf);
+
+ yb = negativePart(yb);
+
+ zf = positivePart(yf);
+
+ yb = negativePart(zb);
+
+ zf = positivePart(zf);
}
+
+
+ if(xb - xf > 0.0)
+ a = xb;
+ else
+ a = xf;
+
+ if(yb - yf > 0.0)
+ b = yb;
+ else
+ b = yf;
+
+ if(zb - zf > 0.0)
+ c = zb;
+ else
+ c = zf;
+
+ d =(1.0 - sqrt(a*a+b*b+c*c)*ihx );
+
+ if(Sign(d) > 0.0 )
+ return Sign(u[cellIndex])*d;
else
-*/ {
- return 0.0;
+ return signui*d;
+}
+
+
+
+template< typename MeshReal,
+ typename Device,
+ typename MeshIndex,
+ typename Real,
+ typename Index >
+
+#ifdef HAVE_CUDA
+__device__
+#endif
+Real parallelGodunovEikonalScheme< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >:: getValueDev( const MeshType& mesh,
+ const IndexType cellIndex,
+ const CoordinatesType& coordinates,
+ const Real* u,
+ const Real& time,
+ const IndexType boundaryCondition) const
+{
+
+/*
+ if ( ((coordinates.x() == 0 && (boundaryCondition & 4)) or
+ (coordinates.x() == mesh.getDimensions().x() - 1 && (boundaryCondition & 2)) or
+ (coordinates.y() == 0 && (boundaryCondition & 8)) or
+ (coordinates.y() == mesh.getDimensions().y() - 1 && (boundaryCondition & 1)))
+ )
+ {
+ return 0.0;
+ }
+
+*/
+ //RealType acc = hx*hy*hx*hy;
+
+ RealType signui;
+ signui = sign(u[cellIndex],this->epsilon);
+
+
+#ifdef HAVE_CUDA
+ //printf("%d : %d ;;;; %d : %d , %f \n",threadIdx.x, mesh.getDimensions().x() , threadIdx.y,mesh.getDimensions().y(), epsilon );
+#endif
+
+ RealType xb = u[cellIndex];
+ RealType xf = -u[cellIndex];
+ RealType yb = u[cellIndex];
+ RealType yf = -u[cellIndex];
+ RealType zb = u[cellIndex];
+ RealType zf = -u[cellIndex];
+ RealType a,b,c,d;
+
+
+ if(coordinates.x() == mesh.getDimensions().x() - 1)
+ xf += u[mesh.template getCellNextToCell<-1,0,0>( cellIndex )];
+ else
+ xf += u[mesh.template getCellNextToCell<1,0,0>( cellIndex )];
+
+ if(coordinates.x() == 0)
+ xb -= u[mesh.template getCellNextToCell<1,0,0>( cellIndex )];
+ else
+ xb -= u[mesh.template getCellNextToCell<-1,0,0>( cellIndex )];
+
+ if(coordinates.y() == mesh.getDimensions().y() - 1)
+ yf += u[mesh.template getCellNextToCell<0,-1,0>( cellIndex )];
+ else
+ yf += u[mesh.template getCellNextToCell<0,1,0>( cellIndex )];
+
+ if(coordinates.y() == 0)
+ yb -= u[mesh.template getCellNextToCell<0,1,0>( cellIndex )];
+ else
+ yb -= u[mesh.template getCellNextToCell<0,-1,0>( cellIndex )];
+
+
+ if(coordinates.z() == mesh.getDimensions().z() - 1)
+ zf += u[mesh.template getCellNextToCell<0,0,-1>( cellIndex )];
+ else
+ zf += u[mesh.template getCellNextToCell<0,0,1>( cellIndex )];
+
+ if(coordinates.z() == 0)
+ zb -= u[mesh.template getCellNextToCell<0,0,1>( cellIndex )];
+ else
+ zb -= u[mesh.template getCellNextToCell<0,0,-1>( cellIndex )];
+
+
+ //xb *= ihx;
+ //xf *= ihx;
+ // yb *= ihy;
+ //yf *= ihy;
+
+ if(signui > 0.0)
+ {
+ xf = negativePart(xf);
+
+ xb = positivePart(xb);
+
+ yf = negativePart(yf);
+
+ yb = positivePart(yb);
+
+ zf = negativePart(zf);
+
+ zb = positivePart(zb);
+
+ }
+ else if(signui < 0.0)
+ {
+
+ xb = negativePart(xb);
+
+ xf = positivePart(xf);
+
+ yb = negativePart(yb);
+
+ zf = positivePart(yf);
+
+ yb = negativePart(zb);
+
+ zf = positivePart(zf);
}
+
+ if(xb - xf > 0.0)
+ a = xb;
+ else
+ a = xf;
+
+ if(yb - yf > 0.0)
+ b = yb;
+ else
+ b = yf;
+
+ if(zb - zf > 0.0)
+ c = zb;
+ else
+ c = zf;
+
+ d =(1.0 - sqrt(a*a+b*b+c*c)*ihx );
+
+ if(Sign(d) > 0.0 )
+ return Sign(u[cellIndex])*d;
+ else
+ return signui*d;
}
--
GitLab