From 798f93750abb73c3f88fcdb3bf9852f010ed4d8e Mon Sep 17 00:00:00 2001
From: Tomas Sobotik <sobotto4@fjfi.cvut.cz>
Date: Wed, 11 Mar 2015 08:23:21 +0100
Subject: [PATCH] added run script and partialy fixed solver
---
examples/hamilton-jacobi-parallel/run | 64 ++++++++
.../hamilton-jacobi-parallel/tnl-err2eoc-2.py | 141 ++++++++++++++++++
.../tnlParallelEikonalSolver_impl.h | 141 ++++++++++++++++--
3 files changed, 333 insertions(+), 13 deletions(-)
create mode 100755 examples/hamilton-jacobi-parallel/run
create mode 100755 examples/hamilton-jacobi-parallel/tnl-err2eoc-2.py
diff --git a/examples/hamilton-jacobi-parallel/run b/examples/hamilton-jacobi-parallel/run
new file mode 100755
index 0000000000..3aece294a9
--- /dev/null
+++ b/examples/hamilton-jacobi-parallel/run
@@ -0,0 +1,64 @@
+#!/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
+
+
diff --git a/examples/hamilton-jacobi-parallel/tnl-err2eoc-2.py b/examples/hamilton-jacobi-parallel/tnl-err2eoc-2.py
new file mode 100755
index 0000000000..f8cde3768e
--- /dev/null
+++ b/examples/hamilton-jacobi-parallel/tnl-err2eoc-2.py
@@ -0,0 +1,141 @@
+#!/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" )
+
diff --git a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h b/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h
index 07435e3629..4194114d69 100644
--- a/examples/hamilton-jacobi-parallel/tnlParallelEikonalSolver_impl.h
+++ b/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;*/
--
GitLab