Added 2D parallel iterative solver for forest-fire problem.

examples/CMakeLists.txt

@@ -4,3 +4,4 @@ add_subdirectory( navier-stokes )
 #add_subdirectory( hamilton-jacobi )
 add_subdirectory( hamilton-jacobi-parallel )
 add_subdirectory( fast-sweeping )
+add_subdirectory( hamilton-jacobi-parallel-map )

examples/hamilton-jacobi-parallel-map/CMakeLists.txt

+set( tnl_hamilton_jacobi_parallel_map_SOURCES
+#     MainBuildConfig.h
+#     tnlParallelMapSolver2D_impl.h
+#     tnlParallelMapSolver.h
+#     parallelMapConfig.h 
+     main.cpp)
+
+
+IF(  BUILD_CUDA ) 
+	CUDA_ADD_EXECUTABLE(hamilton-jacobi-parallel-map${debugExt} main.cu)
+ELSE(  BUILD_CUDA )                
+	ADD_EXECUTABLE(hamilton-jacobi-parallel-map${debugExt} main.cpp)
+ENDIF( BUILD_CUDA )
+target_link_libraries (hamilton-jacobi-parallel-map${debugExt} tnl${debugExt}-${tnlVersion} )
+
+
+INSTALL( TARGETS hamilton-jacobi-parallel-map${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_map_SOURCES}
+#         DESTINATION share/tnl-${tnlVersion}/examples/hamilton-jacobi-parallel-map )

examples/hamilton-jacobi-parallel-map/MainBuildConfig.h

+/***************************************************************************
+                          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/tnlBuildConfigTags.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_ */

examples/hamilton-jacobi-parallel-map/Makefile~

+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) -DDEBUG
+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)
+	$(CUDA_CXX) -o $(TARGET) $(OBJECTS) $(LD_FLAGS)
+
+%.o: %.cpp $(HEADERS)
+	$(CXX) -c -o $@ $(CXX_FLAGS) $<
+
+

examples/hamilton-jacobi-parallel-map/main.cpp

+/***************************************************************************
+                          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 "main.h"

examples/hamilton-jacobi-parallel-map/main.cu

+/***************************************************************************
+                          main.cu  -  description
+                             -------------------
+    begin                : Mar 30 , 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"

examples/hamilton-jacobi-parallel-map/main.h

+/***************************************************************************
+                          main.h  -  description
+                             -------------------
+    begin                : Mar 22 , 2016
+    copyright            : (C) 2016 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 "tnlParallelMapSolver.h"
+#include "parallelMapConfig.h"
+#include "MainBuildConfig.h"
+#include <solvers/tnlBuildConfigTags.h>
+#include <operators/godunov-eikonal/parallelGodunovMap.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;
+   parallelMapConfig< BuildConfig >::configSetup( configDescription );
+
+   if( ! parseCommandLine( argc, argv, configDescription, parameters ) )
+      return false;
+
+
+   tnlDeviceEnum device;
+   device = tnlHostDevice;
+
+   const int& dim = parameters.getParameter< int >( "dim" );
+
+  if(dim == 2)
+  {
+
+	   typedef parallelGodunovMapScheme< tnlGrid<2,double,tnlHost, int>, double, int > SchemeTypeHost;
+		/*#ifdef HAVE_CUDA
+		   typedef parallelGodunovMapScheme< tnlGrid<2,double,tnlCuda, int>, double, int > SchemeTypeDevice;
+		#endif
+		#ifndef HAVE_CUDA*/
+	   typedef parallelGodunovMapScheme< tnlGrid<2,double,tnlHost, int>, double, int > SchemeTypeDevice;
+		/*#endif*/
+
+	   if(device==tnlHostDevice)
+	   {
+		   typedef tnlHost Device;
+
+
+		   tnlParallelMapSolver<2,SchemeTypeHost,SchemeTypeDevice, Device> solver;
+		   if(!solver.init(parameters))
+		   {
+			   cerr << "Solver failed to initialize." << endl;
+			   return EXIT_FAILURE;
+		   }
+		   cout << "-------------------------------------------------------------" << endl;
+		   cout << "Starting solver loop..." << endl;
+		   solver.run();
+	   }
+	   else if(device==tnlCudaDevice )
+	   {
+		   typedef tnlCuda Device;
+		   //typedef parallelGodunovMapScheme< tnlGrid<2,double,Device, int>, double, int > SchemeType;
+
+		   tnlParallelMapSolver<2,SchemeTypeHost,SchemeTypeDevice, Device> solver;
+		   if(!solver.init(parameters))
+		   {
+			   cerr << "Solver failed to initialize." << endl;
+			   return EXIT_FAILURE;
+		   }
+		   cout << "-------------------------------------------------------------" << endl;
+		   cout << "Starting solver loop..." << endl;
+		   solver.run();
+	   }
+  // }
+  }
+
+
+   time(&stop);
+   cout << endl;
+   cout << "Running time was: " << difftime(stop,start) << " .... " << (std::clock() - start2) / (double)(CLOCKS_PER_SEC) << endl;
+   return EXIT_SUCCESS;
+}
+
+

examples/hamilton-jacobi-parallel-map/no-Makefile

+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) -DDEBUG
+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) $<
+
+

examples/hamilton-jacobi-parallel-map/parallelMapConfig.h

+/***************************************************************************
+                          parallelMapConfig.h  -  description
+                             -------------------
+    begin                : Mar 22 , 2016
+    copyright            : (C) 2016 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 HAMILTONJACOBIPARALLELMAPPROBLEMCONFIG_H_
+#define HAMILTONJACOBIPARALLELMAPPROBLEMCONFIG_H_
+
+#include <config/tnlConfigDescription.h>
+
+template< typename ConfigTag >
+class parallelMapConfig
+{
+   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.addRequiredEntry        < tnlString > ( "map", "Gradient map 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 );
+         config.addRequiredEntry        < int > ( "dim", "Dimension of problem.");
+      }
+};
+
+#endif /* HAMILTONJACOBIPARALLELMAPPROBLEMCONFIG_H_ */

examples/hamilton-jacobi-parallel-map/run

+#!/bin/bash
+
+#GRID_SIZES="0897"
+GRID_SIZES="0008 0015 0029 0057 0113 0225 0449"
+#GRID_SIZES="1793"
+
+dimensions=2
+
+size=2
+
+time=3
+
+for grid_size in $GRID_SIZES;
+
+do
+
+	rm -r grid-${grid_size}
+   	mkdir grid-${grid_size}
+   	cd grid-${grid_size}
+
+	tnl-grid-setup --dimensions $dimensions \
+	               --origin-x -1.0 \
+	               --origin-y -1.0 \
+	               --origin-z -1.0 \
+	               --proportions-x $size \
+	               --proportions-y $size \
+	               --proportions-z $size \
+	               --size-x ${grid_size} \
+	               --size-y ${grid_size} \
+	               --size-z ${grid_size}
+
+	tnl-init --test-function sdf-para \
+		     --offset 0.25 \
+	             --output-file init.tnl \
+		     --final-time 0.0 \
+		     --snapshot-period 0.1 \
+
+
+	tnl-init --test-function sdf-para-sdf \
+		     --offset 0.25 \
+	             --output-file sdf.tnl \
+		     --final-time 0.0 \
+		     --snapshot-period 0.1
+
+	hamilton-jacobi-parallel --initial-condition init.tnl \
+	              --cfl-condition 1.0e-1 \
+		      	  --mesh mesh.tnl \
+		     	  --initial-tau 1.0e-3 \
+		      	  --epsilon 1.0 \
+	        	  --delta 0.0 \
+	       	      --stop-time $time \
+		          --scheme godunov \
+		          --subgrid-size 8
+
+        tnl-diff --mesh mesh.tnl --mode sequence --input-files sdf.tnl u-00001.tnl --write-difference yes --output-file ../${grid_size}.diff
+	
+	cd ..
+
+done
+
+
+./tnl-err2eoc-2.py --format txt --size $size *.diff
+
+              

examples/hamilton-jacobi-parallel-map/run~

+#!/bin/bash
+
+#GRID_SIZES="0897"
+#GRID_SIZES="0449"
+GRID_SIZES="1793"
+
+dimensions=2
+
+size=2
+
+time=3
+
+for grid_size in $GRID_SIZES;
+
+do
+
+	rm -r grid-${grid_size}
+   	mkdir grid-${grid_size}
+   	cd grid-${grid_size}
+
+	tnl-grid-setup --dimensions $dimensions \
+	               --origin-x -1.0 \
+	               --origin-y -1.0 \
+	               --origin-z -1.0 \
+	               --proportions-x $size \
+	               --proportions-y $size \
+	               --proportions-z $size \
+	               --size-x ${grid_size} \
+	               --size-y ${grid_size} \
+	               --size-z ${grid_size}
+
+	tnl-init --test-function sdf-para \
+		     --offset 0.25 \
+	             --output-file init.tnl \
+		     --final-time 0.0 \
+		     --snapshot-period 0.1 \
+
+
+	tnl-init --test-function sdf-para-sdf \
+		     --offset 0.25 \
+	             --output-file sdf.tnl \
+		     --final-time 0.0 \
+		     --snapshot-period 0.1
+
+	hamilton-jacobi-parallel --initial-condition init.tnl \
+	              --cfl-condition 1.0e-1 \
+		      	  --mesh mesh.tnl \
+		     	  --initial-tau 1.0e-3 \
+		      	  --epsilon 1.0 \
+	        	  --delta 0.0 \
+	       	      --stop-time $time \
+		          --scheme godunov \
+		          --subgrid-size 8
+
+        tnl-diff --mesh mesh.tnl --mode sequence --input-files sdf.tnl u-00001.tnl --write-difference yes --output-file ../${grid_size}.diff
+	
+	cd ..
+
+done
+
+
+./tnl-err2eoc-2.py --format txt --size $size *.diff
+
+              

examples/hamilton-jacobi-parallel-map/tnl-err2eoc-2.py

+#!/usr/bin/env python
+
+import sys, string, math
+
+arguments = sys. argv[1:]
+format = "txt"
+output_file_name = "eoc-table.txt"
+input_files = []
+verbose = 1
+size = 1.0
+
+i = 0
+while i < len( arguments ):
+   if arguments[ i ] == "--format":
+      format = arguments[ i + 1 ]
+      i = i + 2
+      continue
+   if arguments[ i ] == "--output-file":
+      output_file_name = arguments[ i + 1 ]
+      i = i + 2
+      continue
+   if arguments[ i ] == "--verbose":
+       verbose = float( arguments[ i + 1 ] )
+       i = i +2
+       continue
+   if arguments[ i ] == "--size":
+       size = float( arguments[ i + 1 ] )
+       i = i +2
+       continue
+   input_files. append( arguments[ i ] )
+   i = i + 1
+
+if not verbose == 0:
+   print "Writing to " + output_file_name + " in " + format + "."
+
+h_list = []
+l1_norm_list = []
+l2_norm_list = []
+max_norm_list = []
+items = 0
+
+for file_name in input_files:
+   if not verbose == 0:
+       print "Processing file " + file_name
+   file = open( file_name, "r" )
+   
+   l1_max = 0.0
+   l_max_max = 0.0
+   file.readline();
+   file.readline();
+   for line in file. readlines():
+         data = string. split( line )
+         h_list. append( size/(float(file_name[0:len(file_name)-5] ) - 1.0) )
+         l1_norm_list. append( float( data[ 1 ] ) )
+         l2_norm_list. append( float( data[ 2 ] ) )
+         max_norm_list. append( float( data[ 3 ] ) )
+         items = items + 1
+         if not verbose == 0:
+            print line
+   file. close()
+
+h_width = 12
+err_width = 15
+file = open( output_file_name, "w" )
+if format == "latex":
+      file. write( "\\begin{tabular}{|r|l|l|l|l|l|l|}\\hline\n" )
+      file. write( "\\raisebox{-1ex}[0ex]{$h$}& \n" )
+      file. write( "\\multicolumn{2}{|c|}{\\raisebox{1ex}[3.5ex]{$\\left\| \\cdot \\right\\|_{L_1\\left(\\omega_h;\\left[0,T\\right]\\right)}^{h,\\tau}$}}& \n" )
+      file. write( "\\multicolumn{2}{|c|}{\\raisebox{1ex}[3.5ex]{$\\left\| \\cdot \\right\\|_{L_2\\left(\\omega_h;\left[0,T\\right]\\right)}^{h,\\tau}$}}& \n" )
+      file. write( "\\multicolumn{2}{|c|}{\\raisebox{1ex}[3.5ex]{$\\left\| \\cdot \\right\\|_{L_\\infty\\left(\\omega_h;\\left[0,T\\right]\\right)}^{h,\\tau}$}}\\\\ \\cline{2-7} \n" )
+      file. write( " " + string. rjust( " ", h_width ) + "&" +
+                string. rjust( "Error", err_width ) + "&" +
+                string. rjust( "{\\bf EOC}", err_width ) + "&" +
+                string. rjust( "Error", err_width ) + "&" +
+                string. rjust( "{\\bf EOC}", err_width ) + "&" +
+                string. rjust( "Error.", err_width ) + "&" +
+                string. rjust( "{\\bf EOC}", err_width ) +
+                "\\\\ \\hline \\hline \n")
+if format == "txt":
+    file. write( "+--------------+----------------+----------------+----------------+----------------+----------------+----------------+\n" )
+    file. write( "|       h      |     L1 Err.    |     L1 EOC.    |     L2 Err.    |      L2 EOC    |    MAX Err.    |     MAX EOC    |\n" )
+    file. write( "+==============+================+================+================+================+================+================+\n" )
+                  
+
+i = 0
+while i < items:
+   if i == 0:
+      if format == "latex":
+         file. write( " " + string. ljust( str( h_list[ i ] ), h_width ) + "&" +
+                      string. rjust( "%.2g" % l1_norm_list[ i ], err_width ) + "&" + 
+                      string. rjust( " ", err_width ) + "&"+ 
+                      string. rjust( "%.2g" % l2_norm_list[ i ], err_width ) + "&" +
+                      string. rjust( " ", err_width ) + "&" +
+                      string. rjust( "%.2g" % max_norm_list[ i ], err_width ) + "&" +
+                      string. rjust( " ", err_width ) + "\\\\\n" )
+      if format == "txt":
+         file. write( "| " + string. ljust( str( h_list[ i ] ), h_width ) + " |" + 
+                      string. rjust( "%.2g" % l1_norm_list[ i ], err_width ) + " |" +
+                      string. rjust( " ", err_width ) + " |" +
+                      string. rjust( "%.2g" % l2_norm_list[ i ], err_width ) + " |" +
+                      string. rjust( " ", err_width ) + " |" +
+                      string. rjust( "%.2g" % max_norm_list[ i ], err_width ) + " |" +
+                      string. rjust( " ", err_width ) + " |\n" )
+         file. write( "+--------------+----------------+----------------+----------------+----------------+----------------+----------------+\n" )
+      i = i + 1;
+      continue
+   if h_list[ i ] == h_list[ i - 1 ]:
+      print "Unable to count eoc since h[ " + \
+      str( i ) + " ] = h[ " + str( i - 1 ) + \
+      " ] = " + str( h_list[ i ] ) + ". \n"
+      file. write( " eoc error:  h[ " + \
+      str( i ) + " ] = h[ " + str( i - 1 ) + \
+      " ] = " + str( h_list[ i ] ) + ". \n" )
+   else:
+      h_ratio = math. log( h_list[ i ] / h_list[ i - 1 ] )
+      l1_ratio = math. log( l1_norm_list[ i ] / l1_norm_list[ i - 1 ] )
+      l2_ratio = math. log( l2_norm_list[ i ] / l2_norm_list[ i - 1 ] )
+      max_ratio = math. log( max_norm_list[ i ] / max_norm_list[ i - 1 ] )
+      if format == "latex":
+         file. write( " " + string. ljust( str( h_list[ i ] ), h_width ) + "&" +
+                      string. rjust( "%.2g" % l1_norm_list[ i ], err_width ) + "&" +
+                      string. rjust( "{\\bf " + "%.2g" % ( l1_ratio / h_ratio ) + "}", err_width ) + "&" +
+                      string. rjust( "%.2g" % l2_norm_list[ i ], err_width ) + "&" +
+                      string. rjust( "{\\bf " + "%.2g" % ( l2_ratio / h_ratio ) + "}", err_width ) + "&" +
+                      string. rjust( "%.2g" % max_norm_list[ i ], err_width ) + "&" +
+                      string. rjust( "{\\bf " + "%.2g" % ( max_ratio / h_ratio ) + "}", err_width ) + "\\\\\n" )
+      if format == "txt":
+         file. write( "| " + string. ljust( str( h_list[ i ] ), h_width ) + " |" +
+                      string. rjust( "%.2g" % l1_norm_list[ i ], err_width ) + " |" +
+                      string. rjust( "**" + "%.2g" % ( l1_ratio / h_ratio ) + "**", err_width ) + " |" +
+                      string. rjust( "%.2g" % l2_norm_list[ i ], err_width ) + " |" +
+                      string. rjust( "**" + "%.2g" % ( l2_ratio / h_ratio ) + "**", err_width ) + " |" +
+                      string. rjust( "%.2g" % max_norm_list[ i ], err_width ) + " |" +
+                      string. rjust( "**" + "%.2g" % ( max_ratio / h_ratio ) + "**", err_width ) + " |\n" )
+         file. write( "+--------------+----------------+----------------+----------------+----------------+----------------+----------------+\n" )
+   i = i + 1
+
+if format == "latex":
+   file. write( "\\hline \n" )
+   file. write( "\\end{tabular} \n" )
+    

examples/hamilton-jacobi-parallel-map/tnlParallelMapSolver.h

+/***************************************************************************
+                          tnlParallelMapSolver.h  -  description
+                             -------------------
+    begin                : Mar 22 , 2016
+    copyright            : (C) 2016 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 TNLPARALLELMAPSOLVER_H_
+#define TNLPARALLELMAPSOLVER_H_
+
+#include <config/tnlParameterContainer.h>
+#include <core/vectors/tnlVector.h>
+#include <core/vectors/tnlStaticVector.h>
+#include <functions/tnlMeshFunction.h>
+#include <core/tnlHost.h>
+#include <mesh/tnlGrid.h>
+#include <mesh/grids/tnlGridEntity.h>
+#include <limits.h>
+#include <core/tnlDevice.h>
+
+
+#include <ctime>
+
+#ifdef HAVE_CUDA
+#include <cuda.h>
+#include <core/tnlCuda.h>
+#endif
+
+
+template< int Dimension,
+		  typename SchemeHost,
+		  typename SchemeDevice,
+		  typename Device,
+		  typename RealType = double,
+          typename IndexType = int >
+class tnlParallelMapSolver
+{};
+
+template<typename SchemeHost, typename SchemeDevice, typename Device>
+class tnlParallelMapSolver<2, SchemeHost, SchemeDevice, Device, double, int >
+{
+public:
+
+	typedef SchemeDevice SchemeTypeDevice;
+	typedef SchemeHost SchemeTypeHost;
+	typedef Device DeviceType;
+	typedef tnlVector< double, tnlHost, int > VectorType;
+	typedef tnlVector< int, tnlHost, int > IntVectorType;
+	typedef tnlGrid< 2, double, tnlHost, int > MeshType;
+#ifdef HAVE_CUDA
+	typedef tnlVector< double, tnlHost, int > VectorTypeCUDA;
+	typedef tnlVector< int, tnlHost, int > IntVectorTypeCUDA;
+	typedef tnlGrid< 2, double, tnlHost, int > MeshTypeCUDA;
+#endif
+	tnlParallelMapSolver();
+	bool init( const tnlParameterContainer& parameters );
+	void run();
+
+	void test();
+
+/*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, int subGridID);
+
+
+	tnlMeshFunction<MeshType> u0;
+	VectorType work_u, map_stretched, map;
+	IntVectorType subgridValues, boundaryConditions, unusedCell, calculationsCount;
+	MeshType mesh, subMesh;
+
+//	tnlGridEntity< MeshType, 2, tnlGridEntityNoStencilStorage > Entity;
+
+	SchemeHost schemeHost;
+	SchemeDevice schemeDevice;
+	double delta, tau0, stopTime,cflCondition;
+	int gridRows, gridCols, gridLevels, currentStep, n;
+
+	std::clock_t start;
+	double time_diff;
+
+
+	tnlDeviceEnum device;
+
+	tnlParallelMapSolver<2, SchemeHost, SchemeDevice, Device, double, int >* getSelf()
+	{
+		return this;
+	};
+
+#ifdef HAVE_CUDA
+
+	tnlParallelMapSolver<2, SchemeHost, SchemeDevice, Device, double, int >* cudaSolver;
+
+	double* work_u_cuda;
+	double* map_stretched_cuda;
+
+	int* subgridValues_cuda;
+	int* boundaryConditions_cuda;
+	int* unusedCell_cuda;
+	int* calculationsCount_cuda;
+	double* tmpw;
+	double* tmp_map;
+	//MeshTypeCUDA mesh_cuda, subMesh_cuda;
+	//SchemeDevice scheme_cuda;
+	//double delta_cuda, tau0_cuda, stopTime_cuda,cflCondition_cuda;
+	//int gridRows_cuda, gridCols_cuda, currentStep_cuda, n_cuda;
+
+	int* runcuda;
+	int run_host;
+
+
+	__device__ void getSubgridCUDA2D( const int i, tnlParallelMapSolver<2, SchemeHost, SchemeDevice, Device, double, int >* caller, double* a);
+
+	__device__ void updateSubgridCUDA2D( const int i, tnlParallelMapSolver<2, SchemeHost, SchemeDevice, Device, double, int >* caller, double* a);
+
+	__device__ void insertSubgridCUDA2D( double u, const int i );
+
+	__device__ void runSubgridCUDA2D( int boundaryCondition, double* u, int subGridID);
+
+	/*__global__ void runCUDA();*/
+
+	//__device__ void synchronizeCUDA();
+
+	__device__ int getOwnerCUDA2D( int i) const;
+
+	__device__ int getSubgridValueCUDA2D( int i ) const;
+
+	__device__ void setSubgridValueCUDA2D( int i, int value );
+
+	__device__ int getBoundaryConditionCUDA2D( int i ) const;
+
+	__device__ void setBoundaryConditionCUDA2D( int i, int value );
+
+	//__device__ bool initCUDA( tnlParallelMapSolver<SchemeHost, SchemeDevice, Device, double, int >* cudaSolver);
+
+	/*__global__ void initRunCUDA(tnlParallelMapSolver<Scheme, double, tnlHost, int >* caller);*/
+
+#endif
+
+};
+
+
+
+
+
+
+
+	template<typename SchemeHost, typename SchemeDevice, typename Device>
+	class tnlParallelMapSolver<3, SchemeHost, SchemeDevice, Device, double, int >
+	{
+	public:
+
+		typedef SchemeDevice SchemeTypeDevice;
+		typedef SchemeHost SchemeTypeHost;
+		typedef Device DeviceType;
+		typedef tnlVector< double, tnlHost, int > VectorType;
+		typedef tnlVector< int, tnlHost, int > IntVectorType;
+		typedef tnlGrid< 3, double, tnlHost, int > MeshType;
+	#ifdef HAVE_CUDA
+		typedef tnlVector< double, tnlHost, int > VectorTypeCUDA;
+		typedef tnlVector< int, tnlHost, int > IntVectorTypeCUDA;
+		typedef tnlGrid< 3, double, tnlHost, int > MeshTypeCUDA;
+	#endif
+		tnlParallelMapSolver();
+		bool init( const tnlParameterContainer& parameters );
+		void run();
+
+		void test();
+
+	/*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, int subGridID);
+
+
+		tnlMeshFunction<MeshType> u0;
+		VectorType work_u;
+		IntVectorType subgridValues, boundaryConditions, unusedCell, calculationsCount;
+		MeshType mesh, subMesh;
+		SchemeHost schemeHost;
+		SchemeDevice schemeDevice;
+		double delta, tau0, stopTime,cflCondition;
+		int gridRows, gridCols, gridLevels, currentStep, n;
+
+		std::clock_t start;
+		double time_diff;
+
+
+		tnlDeviceEnum device;
+
+		tnlParallelMapSolver<3, SchemeHost, SchemeDevice, Device, double, int >* getSelf()
+		{
+			return this;
+		};
+
+#ifdef HAVE_CUDA
+
+	tnlParallelMapSolver<3, SchemeHost, SchemeDevice, Device, double, int >* cudaSolver;
+
+	double* work_u_cuda;
+
+	int* subgridValues_cuda;
+	int* boundaryConditions_cuda;
+	int* unusedCell_cuda;
+	int* calculationsCount_cuda;
+	double* tmpw;
+	//MeshTypeCUDA mesh_cuda, subMesh_cuda;
+	//SchemeDevice scheme_cuda;
+	//double delta_cuda, tau0_cuda, stopTime_cuda,cflCondition_cuda;
+	//int gridRows_cuda, gridCols_cuda, currentStep_cuda, n_cuda;
+
+	int* runcuda;
+	int run_host;
+
+
+	__device__ void getSubgridCUDA3D( const int i, tnlParallelMapSolver<3, SchemeHost, SchemeDevice, Device, double, int >* caller, double* a);
+
+	__device__ void updateSubgridCUDA3D( const int i, tnlParallelMapSolver<3, SchemeHost, SchemeDevice, Device, double, int >* caller, double* a);
+
+	__device__ void insertSubgridCUDA3D( double u, const int i );
+
+	__device__ void runSubgridCUDA3D( int boundaryCondition, double* u, int subGridID);
+
+	/*__global__ void runCUDA();*/
+
+	//__device__ void synchronizeCUDA();
+
+	__device__ int getOwnerCUDA3D( int i) const;
+
+	__device__ int getSubgridValueCUDA3D( int i ) const;
+
+	__device__ void setSubgridValueCUDA3D( int i, int value );
+
+	__device__ int getBoundaryConditionCUDA3D( int i ) const;
+
+	__device__ void setBoundaryConditionCUDA3D( int i, int value );
+
+	//__device__ bool initCUDA( tnlParallelMapSolver<SchemeHost, SchemeDevice, Device, double, int >* cudaSolver);
+
+	/*__global__ void initRunCUDA(tnlParallelMapSolver<Scheme, double, tnlHost, int >* caller);*/
+
+#endif
+
+};
+
+
+
+
+
+
+#ifdef HAVE_CUDA
+template <typename SchemeHost, typename SchemeDevice, typename Device>
+__global__ void runCUDA2D(tnlParallelMapSolver<2, SchemeHost, SchemeDevice, Device, double, int >* caller);
+
+template <typename SchemeHost, typename SchemeDevice, typename Device>
+__global__ void initRunCUDA2D(tnlParallelMapSolver<2, SchemeHost, SchemeDevice, Device, double, int >* caller);
+
+template <typename SchemeHost, typename SchemeDevice, typename Device>
+__global__ void initCUDA2D( tnlParallelMapSolver<2, SchemeHost, SchemeDevice, Device, double, int >* cudaSolver, double* ptr, int * ptr2, int* ptr3, double* tmp_map_ptr);
+
+template <typename SchemeHost, typename SchemeDevice, typename Device>
+__global__ void synchronizeCUDA2D(tnlParallelMapSolver<2, SchemeHost, SchemeDevice, Device, double, int >* cudaSolver);
+
+template <typename SchemeHost, typename SchemeDevice, typename Device>
+__global__ void synchronize2CUDA2D(tnlParallelMapSolver<2, SchemeHost, SchemeDevice, Device, double, int >* cudaSolver);
+
+
+
+
+
+
+
+template <typename SchemeHost, typename SchemeDevice, typename Device>
+__global__ void runCUDA3D(tnlParallelMapSolver<3, SchemeHost, SchemeDevice, Device, double, int >* caller);
+
+template <typename SchemeHost, typename SchemeDevice, typename Device>
+__global__ void initRunCUDA3D(tnlParallelMapSolver<3, SchemeHost, SchemeDevice, Device, double, int >* caller);
+
+template <typename SchemeHost, typename SchemeDevice, typename Device>
+__global__ void initCUDA3D( tnlParallelMapSolver<3, SchemeHost, SchemeDevice, Device, double, int >* cudaSolver, double* ptr, int * ptr2, int* ptr3);
+
+template <typename SchemeHost, typename SchemeDevice, typename Device>
+__global__ void synchronizeCUDA3D(tnlParallelMapSolver<3, SchemeHost, SchemeDevice, Device, double, int >* cudaSolver);
+
+template <typename SchemeHost, typename SchemeDevice, typename Device>
+__global__ void synchronize2CUDA3D(tnlParallelMapSolver<3, SchemeHost, SchemeDevice, Device, double, int >* cudaSolver);
+#endif
+
+
+__device__
+double fabsMin( double x, double y)
+{
+	double fx = abs(x);
+
+	if(Min(fx,abs(y)) == fx)
+		return x;
+	else
+		return y;
+}
+
+__device__
+double atomicFabsMin(double* address, double val)
+{
+	unsigned long long int* address_as_ull =
+						  (unsigned long long int*)address;
+	unsigned long long int old = *address_as_ull, assumed;
+	do {
+		assumed = old;
+			old = atomicCAS(address_as_ull, assumed,__double_as_longlong( fabsMin(__longlong_as_double(assumed),val) ));
+	} while (assumed != old);
+	return __longlong_as_double(old);
+}
+
+
+#include "tnlParallelMapSolver2D_impl.h"
+#endif /* TNLPARALLELMAPSOLVER_H_ */

examples/hamilton-jacobi-parallel/CMakeLists.txt

 set( tnl_hamilton_jacobi_parallel_SOURCES
 #     MainBuildConfig.h
-#     tnlParallelEikonalSolver_impl.h
+#     tnlParallelEikonalSolver2D_impl.h
+#     tnlParallelEikonalSolver3D_impl.h
 #     tnlParallelEikonalSolver.h
 #     parallelEikonalConfig.h 
      main.cpp)

examples/hamilton-jacobi-parallel/parallelEikonalConfig.h

 /***************************************************************************
-                          hamiltonJacobiProblemConfig.h  -  description
+                          parallelEikonalConfig.h  -  description
                              -------------------
     begin                : Oct 5, 2014
     copyright            : (C) 2014 by Tomas Sobotik
@@ -15,8 +15,8 @@
  *                                                                         *
  ***************************************************************************/
 
-#ifndef HAMILTONJACOBIPROBLEMCONFIG_H_
-#define HAMILTONJACOBIPROBLEMCONFIG_H_
+#ifndef HAMILTONJACOBIPARALLELEIKONALPROBLEMCONFIG_H_
+#define HAMILTONJACOBIPARALLELEIKONALPROBLEMCONFIG_H_
 
 #include <config/tnlConfigDescription.h>
 
@@ -43,4 +43,4 @@ class parallelEikonalConfig
       }
 };
 
-#endif /* HAMILTONJACOBIPROBLEMCONFIG_H_ */
+#endif /* HAMILTONJACOBIPARALLELEIKONALPROBLEMCONFIG_H_ */

src/operators/godunov-eikonal/CMakeLists.txt

 set( headers godunovEikonal.h
 	     parallelGodunovEikonal.h
+	     parallelGodunovMap.h
 	     godunovEikonal1D_impl.h
              godunovEikonal2D_impl.h
              godunovEikonal3D_impl.h
              parallelGodunovEikonal1D_impl.h
              parallelGodunovEikonal2D_impl.h
-             parallelGodunovEikonal3D_impl.h)
+             parallelGodunovEikonal3D_impl.h
+             parallelGodunovMap2D_impl.h)
 
 SET( CURRENT_DIR ${CMAKE_SOURCE_DIR}/src/implementation/operators/godunov-eikonal )
 

src/operators/godunov-eikonal/parallelGodunovEikonal.h

 /***************************************************************************
-                          godunov.h  -  description
+                          parallelGodunovEikonal.h  -  description
                              -------------------
     begin                : Dec 1 , 2014
     copyright            : (C) 2014 by Tomas Sobotik
@@ -14,8 +14,8 @@
  *                                                                         *
  ***************************************************************************/
 
-#ifndef GODUNOVEIKONAL_H_
-#define GODUNOVEIKONAL_H_
+#ifndef PARALLELGODUNOVEIKONAL_H_
+#define PARALLELGODUNOVEIKONAL_H_
 
 #include <matrices/tnlCSRMatrix.h>
 #include <solvers/preconditioners/tnlDummyPreconditioner.h>
@@ -266,4 +266,4 @@ protected:
 #include <operators/godunov-eikonal/parallelGodunovEikonal3D_impl.h>
 
 
-#endif /* GODUNOVEIKONAL_H_ */
+#endif /* PARALLELGODUNOVEIKONAL_H_ */

src/operators/godunov-eikonal/parallelGodunovEikonal1D_impl.h

@@ -14,8 +14,8 @@
  *                                                                         *
  ***************************************************************************/
 
-#ifndef GODUNOVEIKONAL1D_IMPL_H_
-#define GODUNOVEIKONAL1D_IMPL_H_
+#ifndef PARALLELGODUNOVEIKONAL1D_IMPL_H_
+#define PARALLELGODUNOVEIKONAL1D_IMPL_H_
 
 
 
@@ -167,4 +167,4 @@ Real parallelGodunovEikonalScheme< tnlGrid< 1, MeshReal, Device, MeshIndex >, Re
 }
 
 
-#endif /* GODUNOVEIKONAL1D_IMPL_H_ */
+#endif /* PARALLELGODUNOVEIKONAL1D_IMPL_H_ */

src/operators/godunov-eikonal/parallelGodunovEikonal2D_impl.h

@@ -14,8 +14,8 @@
  *                                                                         *
  ***************************************************************************/
 
-#ifndef GODUNOVEIKONAL2D_IMPL_H_
-#define GODUNOVEIKONAL2D_IMPL_H_
+#ifndef PARALLELGODUNOVEIKONAL2D_IMPL_H_
+#define PARALLELGODUNOVEIKONAL2D_IMPL_H_
 
 
 template< typename MeshReal,
@@ -481,4 +481,4 @@ Real parallelGodunovEikonalScheme< tnlGrid< 2, MeshReal, Device, MeshIndex >, Re
 }
 
 
-#endif /* GODUNOVEIKONAL2D_IMPL_H_ */
+#endif /* PARALLELGODUNOVEIKONAL2D_IMPL_H_ */