added run script and partialy fixed solver

examples/hamilton-jacobi-parallel/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/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/tnlParallelEikonalSolver_impl.h

@@ -157,6 +157,30 @@ void tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::run()
 				{
 					insertSubgrid( runSubgrid(8, getSubgrid(i),i), i);
 				}
+
+
+				if( (getBoundaryCondition(i) & 1) || (getBoundaryCondition(i) & 2) )
+				{
+					//cout << "3 @ " << getBoundaryCondition(i) << endl;
+					insertSubgrid( runSubgrid(3, getSubgrid(i),i), i);
+				}
+				if( (getBoundaryCondition(i) & 1) || (getBoundaryCondition(i) & 4) )
+				{
+					//cout << "5 @ " << getBoundaryCondition(i) << endl;
+					insertSubgrid( runSubgrid(5, getSubgrid(i),i), i);
+				}
+				if( (getBoundaryCondition(i) & 8) || (getBoundaryCondition(i) & 2) )
+				{
+					//cout << "10 @ " << getBoundaryCondition(i) << endl;
+					insertSubgrid( runSubgrid(10, getSubgrid(i),i), i);
+				}
+				if( (getBoundaryCondition(i) & 8) || (getBoundaryCondition(i) & 4) )
+				{
+					//cout << "12 @ " << getBoundaryCondition(i) << endl;
+					insertSubgrid( runSubgrid(12, getSubgrid(i),i), i);
+				}
+
+
 				/*if(getBoundaryCondition(i))
 				{
 					insertSubgrid( runSubgrid(15, getSubgrid(i),i), i);
@@ -513,7 +537,7 @@ tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::runSubgrid( int boundary
 	/**/
 
 
-	bool tmp = false;
+/*	bool tmp = false;
 	for(int i = 0; i < u.getSize(); i++)
 	{
 		if(u[0]*u[i] <= 0.0)
@@ -525,44 +549,131 @@ tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::runSubgrid( int boundary
 	{}
 	else if(boundaryCondition == 4)
 	{
-		for(int i = 0; i < u.getSize(); i=subMesh.getCellYSuccessor(i))
+		int i;
+		for(i = 0; i < u.getSize() - subMesh.getDimensions().x() ; i=subMesh.getCellYSuccessor(i))
 		{
-			for(int j = i; j < subMesh.getDimensions().x(); j=subMesh.getCellXSuccessor(j))
+			int j;
+			for(j = i; j < subMesh.getDimensions().x() - 1; j=subMesh.getCellXSuccessor(j))
 			{
 				u[j] = u[i];
 			}
+			u[j] = u[i];
 		}
+		int j;
+		for(j = i; j < subMesh.getDimensions().x() - 1; j=subMesh.getCellXSuccessor(j))
+		{
+			u[j] = u[i];
+		}
+		u[j] = u[i];
 	}
 	else if(boundaryCondition == 8)
 	{
-		for(int i = 0; i < subMesh.getDimensions().x(); i=subMesh.getCellXSuccessor(i))
+		int i;
+		for(i = 0; i < subMesh.getDimensions().x() - 1; i=subMesh.getCellXSuccessor(i))
 		{
-			for(int j = i; j < u.getSize(); j=subMesh.getCellYSuccessor(j))
+			int j;
+			for(j = i; j < u.getSize() - subMesh.getDimensions().x(); j=subMesh.getCellYSuccessor(j))
 			{
 				u[j] = u[i];
 			}
+			u[j] = u[i];
+		}
+		int j;
+		for(j = i; j < u.getSize() - subMesh.getDimensions().x(); j=subMesh.getCellYSuccessor(j))
+		{
+			u[j] = u[i];
 		}
+		u[j] = u[i];
+
 	}
 	else if(boundaryCondition == 2)
 	{
-		for(int i = subMesh.getDimensions().x() - 1; i < u.getSize(); i=subMesh.getCellYSuccessor(i))
+		int i;
+		for(i = subMesh.getDimensions().x() - 1; i < u.getSize() - subMesh.getDimensions().x() ; i=subMesh.getCellYSuccessor(i))
 		{
-			for(int j = i; j > (i-1)*subMesh.getDimensions().x(); j=subMesh.getCellXPredecessor(j))
+			int j;
+			for(j = i; j > (i-1)*subMesh.getDimensions().x(); j=subMesh.getCellXPredecessor(j))
 			{
 				u[j] = u[i];
 			}
+			u[j] = u[i];
 		}
+		int j;
+		for(j = i; j > (i-1)*subMesh.getDimensions().x(); j=subMesh.getCellXPredecessor(j))
+		{
+			u[j] = u[i];
+		}
+		u[j] = u[i];
 	}
 	else if(boundaryCondition == 1)
 	{
-		for(int i = (subMesh.getDimensions().y() - 1)*subMesh.getDimensions().x(); i < u.getSize(); i=subMesh.getCellXSuccessor(i))
+		int i;
+		for(i = (subMesh.getDimensions().y() - 1)*subMesh.getDimensions().x(); i < u.getSize() - 1; i=subMesh.getCellXSuccessor(i))
 		{
-			for(int j = i; j > 0; j=subMesh.getCellYPredecessor(j))
+			int j;
+			for(j = i; j >=subMesh.getDimensions().x(); j=subMesh.getCellYPredecessor(j))
 			{
 				u[j] = u[i];
 			}
+			u[j] = u[i];
 		}
+		int j;
+		for(j = i; j >=subMesh.getDimensions().x(); j=subMesh.getCellYPredecessor(j))
+		{
+			u[j] = u[i];
+		}
+		u[j] = u[i];
+	}
+*/
+	/**/
+
+
+
+	bool tmp = false;
+	for(int i = 0; i < u.getSize(); i++)
+	{
+		if(u[0]*u[i] <= 0.0)
+			tmp=true;
 	}
+	//if(this->currentStep + 3 < getSubgridValue(subGridID))
+		//tmp = true;
+
+
+	if(tmp)
+	{}
+
+
+	//north - 1, east - 2, west - 4, south - 8
+	else if(boundaryCondition == 4)
+	{
+		for(int i = 0; i < this->n; i++)
+			for(int j = 1;j < this->n; j++)
+				if(fabs(u[i*this->n + j]) <  fabs(u[i*this->n]))
+				u[i*this->n + j] = u[i*this->n];
+	}
+	else if(boundaryCondition == 2)
+	{
+		for(int i = 0; i < this->n; i++)
+			for(int j =0 ;j < this->n -1; j++)
+				if(fabs(u[i*this->n + j]) < fabs(u[(i+1)*this->n - 1]))
+				u[i*this->n + j] = u[(i+1)*this->n - 1];
+	}
+	else if(boundaryCondition == 1)
+	{
+		for(int j = 0; j < this->n; j++)
+			for(int i = 0;i < this->n - 1; i++)
+				if(fabs(u[i*this->n + j]) < fabs(u[j + this->n*(this->n - 1)]))
+				u[i*this->n + j] = u[j + this->n*(this->n - 1)];
+	}
+	else if(boundaryCondition == 8)
+	{
+		for(int j = 0; j < this->n; j++)
+			for(int i = 1;i < this->n; i++)
+				if(fabs(u[i*this->n + j]) < fabs(u[j]))
+				u[i*this->n + j] = u[j];
+	}
+
+
 
 	/**/
 
@@ -575,8 +686,9 @@ tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::runSubgrid( int boundary
    double finalTime = this->stopTime;// + 3.0*(u.max() - u.min());
    if( time + currentTau > finalTime ) currentTau = finalTime - time;
 
-   double maxResidue( 0.0 );
-   while( time < finalTime)
+   double maxResidue( 1.0 );
+   double lastResidue( 10000.0 );
+   while( time < finalTime /*|| maxResidue > subMesh.getHx()*/)
    {
       /****
        * Compute the RHS
@@ -591,6 +703,8 @@ tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::runSubgrid( int boundary
 
       if( this -> cflCondition * maxResidue != 0.0)
     	  currentTau =  this -> cflCondition / maxResidue;
+      /*if(maxResidue > lastResidue)
+    	  currentTau *=(1.0/10.0);*/
 
 
       if( time + currentTau > finalTime ) currentTau = finalTime - time;
@@ -610,9 +724,10 @@ tnlParallelEikonalSolver<Scheme, double, tnlHost, int>::runSubgrid( int boundary
       time += currentTau;
 
       //cout << '\r' << flush;
-      //cout << maxResidue << "   " << currentTau << " @ " << time << flush;
+     //cout << maxResidue << "   " << currentTau << " @ " << time << flush;
+     lastResidue = maxResidue;
    }
-
+   //cout << "Time: " << time << ", Res: " << maxResidue <<endl;
 	/*if (u.max() > 0.0)
 		this->stopTime /=(double) this->gridCols;*/