Index: /anuga_validation/circular_island/get_compared_graphs.py
===================================================================
--- /anuga_validation/circular_island/get_compared_graphs.py	(revision 5131)
+++ /anuga_validation/circular_island/get_compared_graphs.py	(revision 5131)
@@ -0,0 +1,31 @@
+
+
+
+# Module imports
+#from anuga.shallow_water import Domain
+#from anuga.shallow_water import Reflective_boundary
+#from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
+#from anuga.shallow_water import Dirichlet_boundary, Time_boundary
+#from anuga.shallow_water.data_manager import get_maximum_inundation_data, start_screen_catcher, copy_code_files
+from anuga.abstract_2d_finite_volumes.util import file_function, sww2csv_gauges,csv2timeseries_graphs
+#from anuga.pmesh.mesh_interface import create_mesh_from_regions
+#from anuga.utilities.polygon import read_polygon#, plot_polygons
+#from math import cos,pi,sin,tan#,sqrt
+#from Numeric import array, zeros, Float, allclose,resize,sqrt
+#from anuga.shallow_water.data_manager import csv2dict
+#from time import localtime, strftime, gmtime
+#from anuga.utilities.system_tools import get_user_name, get_host_name
+from os import sep, environ, getenv
+#
+anuga_dir = home+'anuga_validation'+sep+'circular_island_tsunami_benchmark'+sep+'anuga'+sep
+
+csv2timeseries_graphs(directories_dic={anuga_dir+sep+'boundaries'+sep:['wavetank',0, 0],
+                                       anuga_dir+'outputs'+sep+'20080229_004328_res_0.0025_nbartzis'+sep:['Fixed Wave',0,0]},
+                            output_dir=anuga_dir+'outputs'+sep+'20080229_004328_res_0.0025_nbartzis'+sep,
+                            base_name='gauge_',
+                            plot_numbers='',
+                            quantities=['stage'],
+                            extra_plot_name='',
+                            assess_all_csv_files=True,                            
+                            create_latex=False,
+                            verbose=True)
Index: /anuga_validation/circular_island/run_circular.py
===================================================================
--- /anuga_validation/circular_island/run_circular.py	(revision 5130)
+++ /anuga_validation/circular_island/run_circular.py	(revision 5131)
@@ -302,17 +302,6 @@
                    use_cache = True)
 
-#csv2timeseries_graphs(directories_dic={anuga_dir+sep+'boundaries'+sep:['wavetank',0, 0],
-#                                       anuga_dir+'outputs'+sep+'20080219_233611_res_0.05_nbartzis'+sep:['Fixed Wave',0,0]},
-#                            output_dir=anuga_dir+'outputs'+sep+'20080219_233611_res_0.05_nbartzis'+sep,
-#                            base_name='gauge_',
-#                            plot_numbers='',
-#                            quantities=['stage'],
-#                            extra_plot_name='',
-#                            assess_all_csv_files=True,                            
-#                            create_latex=False,
-#                            verbose=True)
-
-
-
-
-
+
+
+
+
Index: /anuga_validation/circular_island/sqrt_table_run_circular.py
===================================================================
--- /anuga_validation/circular_island/sqrt_table_run_circular.py	(revision 5131)
+++ /anuga_validation/circular_island/sqrt_table_run_circular.py	(revision 5131)
@@ -0,0 +1,322 @@
+"""Validation of the AnuGA implementation of the shallow water wave equation.
+
+This script sets up 2D circular island.....
+
+use 'res' parameter of 1 to make run fast.
+
+"""
+
+# Module imports
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Reflective_boundary
+from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
+from anuga.shallow_water import Dirichlet_boundary, Time_boundary
+from anuga.shallow_water.data_manager import get_maximum_inundation_data, start_screen_catcher, copy_code_files
+from anuga.abstract_2d_finite_volumes.util import file_function, sww2csv_gauges,csv2timeseries_graphs
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
+from anuga.utilities.polygon import read_polygon#, plot_polygons
+from math import cos,pi,sin,tan,sqrt
+from Numeric import array, zeros, Float, allclose,resize
+from anuga.shallow_water.data_manager import csv2dict
+from time import localtime, strftime, gmtime
+from anuga.utilities.system_tools import get_user_name, get_host_name
+from os import sep, environ, getenv
+#-------------------------
+# Create Domain from mesh
+#-------------------------
+
+user = get_user_name()
+home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent dir  
+time = strftime('%Y%m%d_%H%M%S',localtime())
+
+res=0.05
+
+dir_comment = time+'_res_'+str(res)+'_'+str(user)
+
+anuga_dir = home+'anuga_validation'+sep+'circular_island_tsunami_benchmark'+sep+'anuga'+sep
+output_dir = anuga_dir+'outputs'+sep+dir_comment+sep
+copy_code_files(output_dir,__file__)
+
+start_screen_catcher(output_dir)
+
+def get_xy(x=0,y=0,r=1,angle=1):
+    #return x and y co-ordinates 
+    x1 = x + cos(angle)*r
+    y1 = y + sin(-angle)*r
+#    print "get",x,y,r,x1,y1,"angle",angle, cos(angle),sin(angle)
+    return x1,y1,
+
+def prepare_timeboundary(filename):
+    """Converting benchmark 2 time series to NetCDF tms file.
+    This is a 'throw-away' code taylor made for files like
+    'Benchmark_2_input.txt' from the LWRU2 benchmark
+    """
+
+    textversion = filename[:-4] + '.txt'
+
+    print 'Preparing time boundary from %s' %textversion
+    from Numeric import array
+
+    fid = open(textversion)
+
+    #Skip first line
+    line = fid.readline()
+
+    #Read remaining lines
+    lines = fid.readlines()
+    fid.close()
+
+
+    N = len(lines)
+    T = zeros(N, Float)  #Time
+    Q = zeros(N, Float)  #Values
+
+    for i, line in enumerate(lines):
+        fields = line.split()
+#        print 'fields',i,line
+
+        T[i] = float(fields[0])
+        Q[i] = float(fields[1])
+
+
+    #Create tms file
+    from Scientific.IO.NetCDF import NetCDFFile
+
+    print 'Writing to', filename
+    fid = NetCDFFile(filename[:-4] + '.tms', 'w')
+
+    fid.institution = 'Geoscience Australia'
+    fid.description = 'Input wave for Benchmark 2'
+    fid.starttime = 0.0
+    fid.createDimension('number_of_timesteps', len(T))
+    fid.createVariable('time', Float, ('number_of_timesteps',))
+    fid.variables['time'][:] = T
+
+    fid.createVariable('stage', Float, ('number_of_timesteps',))
+    fid.variables['stage'][:] = Q[:]
+
+    fid.createVariable('xmomentum', Float, ('number_of_timesteps',))
+    fid.variables['xmomentum'][:] = 0.0
+
+    fid.createVariable('ymomentum', Float, ('number_of_timesteps',))
+    fid.variables['ymomentum'][:] = 0.0
+
+    fid.close()
+
+
+# angle to make circular interior polygons
+angles1 = (0.0, 22.5, 45.0, 67.5, 90.0, 112.5, 135.0, 157.5, 180.0, 202.5, 225.0,
+          247.5, 270.0, 292.5, 315.0, 337.5)
+
+center_x = 12.96
+center_y = 13.80
+
+#d=0.75
+#L=3.6 and therefore the gauges 1,2,3 and 4 are located 1.8 from the toe of the island
+# 13.8 - 3.6(radii of island) - 1.8 = 8.4
+L=8.4
+
+#create polygons (circles) to have higher resolution
+poly = []
+poly1 = []
+internal_poly=[]
+for i,angle in enumerate(angles1):
+    #convert Degs to Rads
+    angle = ((angle-180)/180.0)*pi
+
+    p = get_xy(center_x,center_y,4.0,angle)
+    poly1.append([p[0],p[1]])
+
+poly.append(poly1)
+internal_poly.append([poly1,.1*res])
+
+poly2 = []
+for i,angle in enumerate(angles1):
+    #convert Degs to Rads
+    angle = ((angle-180)/180.0)*pi
+    p = get_xy(center_x,center_y,2.6,angle)
+    poly2.append([p[0],p[1]])
+    
+poly.append(poly2)
+internal_poly.append([poly2,.01*res])
+
+poly3 = []
+for i,angle in enumerate(angles1):
+    #convert Degs to Rads
+    angle = ((angle-180)/180.0)*pi
+#    p = get_xy(center_x,center_y,1.6,angle)
+    p = get_xy(center_x,center_y,1.5,angle)
+    poly3.append([p[0],p[1]])
+    
+poly.append(poly3)
+internal_poly.append([poly3,1*res])
+
+#plot_polygons(poly,'poly.png')
+
+poly_all= [[0,L],[0,25],[30,25],[30,L]]
+create_mesh_from_regions(poly_all,
+                         boundary_tags={'wall': [0,2],'wave': [3], 'buffer':[1]},
+                         maximum_triangle_area=1*res,
+                         interior_regions=internal_poly,
+                         filename='temp.msh',
+                         use_cache=False,
+                         verbose=True)
+
+domain = Domain('temp.msh', use_cache=True, verbose=True)
+print domain.statistics()
+
+#-------------------------
+# Initial Conditions
+#-------------------------
+domain.set_quantity('friction', 0.0)
+#water_height = 0.322
+water_height = 0.0
+domain.set_quantity('stage', water_height)
+
+
+attribute_dic, title_index_dic = csv2dict('sqrt_table.csv')
+#print attribute_dic
+#print attribute_dic['number'][10]
+
+def circular_island_elevation(x,y):
+
+    list_xy_square = (center_x-x)**2+(center_y-L-y)**2
+    print 'x',min(x),max(x)
+    print 'y',min(y),max(y)
+    list_z=[]
+    for xy in list_xy_square:
+        # this finds the int that relates to the correct sqrt in the table
+        nn=int(round(xy,2)*100)
+        #the square root is to find the distance from the center to the current xy
+        distance = float(attribute_dic['sqrt'][nn])
+#        zz=-float(attribute_dic['sqrt'][nn])
+#        print nn,zz
+        
+        # slope of cone is 1/4 and has a radii of 3.6 therefore the cone is 0.9 high. 
+        # So to have the elevation positive 0.90 has been added. Furthermore the 
+        #wavetank is submerged .32m so 0.9 - 0.32=.0578
+#        z = zz*.25+0.578
+
+#        z = -distance*.25+0.9
+#        z = -distance*.25+0.578
+        z = -distance*.25+0.9-0.322
+#        z = zz*.25+0.8975
+
+        if z <0-0.322:
+            z=0-0.322
+        if z > .625-0.322:
+            z=0.625-0.322
+#        print 'hello',z
+        list_z.append(z)
+    return list_z
+
+domain.set_quantity('elevation', circular_island_elevation, verbose=True)
+
+
+#-------------------------
+# Set simulation parameters
+#-------------------------
+model_name='wavetank_model'
+domain.set_name(model_name)  # Name of output sww file
+domain.set_datadir(output_dir)
+domain.set_default_order(1)               # Apply second order scheme 
+#domain.set_all_limiters(0.9)              # Max second order scheme (old lim)
+domain.set_minimum_storable_height(0.001) # Don't store w < 0.001m
+domain.set_maximum_allowed_speed(0.1)     # Allow a little runoff (0.1 is OK)
+domain.tight_slope_limiters = 1
+domain.beta_h = 0.0
+sww_file=output_dir+sep+model_name+'.sww'
+
+
+#-------------------------
+# Boundary Conditions
+#-------------------------
+
+# Create boundary function from timeseries provided in file
+
+boundary_filename="ts2cnew1_input_20_80sec_new"
+prepare_timeboundary(boundary_filename+'.txt')
+
+function = file_function(boundary_filename+'.tms',
+                         domain, verbose=True)
+
+# Create and assign boundary objects
+Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
+
+Br = Reflective_boundary(domain)
+Bd = Dirichlet_boundary([water_height,0,0])
+#domain.set_boundary({'wave': Bts, 'wall': Br, 'buffer':Bd})
+domain.set_boundary({'wave': Bts, 'wall': Bd, 'buffer':Bd})
+#domain.set_time(20.0)
+#-------------------------
+# Evolve through time
+#-------------------------
+import time
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = 1, finaltime = 28):
+    domain.write_time()
+#    domain.write_time(track_speeds=False)
+
+for t in domain.evolve(yieldstep = 0.1, finaltime = 36
+#for t in domain.evolve(yieldstep = 0.5, finaltime = 36
+                       ,skip_initial_step = True): 
+    domain.write_time()
+   
+for t in domain.evolve(yieldstep = 1, finaltime = 45
+                       ,skip_initial_step = True): 
+    domain.write_time()
+
+print 'That took %.2f seconds' %(time.time()-t0)
+
+
+
+#define the segments to find the run-up height to compare to wave tank study
+angles = (0.0, 22.5, 45.0, 67.5, 75.0, 80.0, 85.0, 87.5, 90.0, 92.5,
+          95.0, 100.0, 105.0, 112.5, 135.0, 157.5, 180.0, 202.5, 225.0,
+          247.5, 270.0, 292.5, 315.0, 337.5)
+r1 = 1.5
+r2 = 2.6
+#degrees to check
+d=1
+
+poly_segment = []
+
+for i,angle in enumerate(angles):
+    angle = ((angle-180)/180.0)*pi
+    p1 = get_xy(center_x,center_y,r1,angle-0.01745*d)
+    p2 = get_xy(center_x,center_y,r2,angle-0.01745*d)
+    p3 = get_xy(center_x,center_y,r2,angle+0.01745*d)
+    p4 = get_xy(center_x,center_y,r1,angle+0.01745*d)
+#    print p1,p2,p3,p4,angle
+    poly_segment =[[p1[0],p1[1]],[p2[0],p2[1]],[p3[0],p3[1]],[p4[0],p4[1]]]
+#    print i, poly_segment, angle
+    #print i,poly[i]
+    
+    run_up, location = get_maximum_inundation_data(filename=sww_file,polygon=poly_segment, verbose=False)
+    
+    print 'maxd %.3f, %.6f, %.6f, %.5s, %.5s '%(((angle/pi)*180)+180, run_up*100, (run_up-water_height)*100, location[0],location[1]) 
+
+#sww_file=anuga_dir+'outputs'+sep+'20080220_010628_res_0.05_nbartzis'+sep+'wavetank_model.sww'
+#sww_file=anuga_dir+'outputs'+sep+'20080219_233611_res_0.05_nbartzis'+sep+'wavetank_model.sww'
+sww2csv_gauges(sww_file=sww_file,
+                   gauge_file='gauges.csv',
+                   quantities = ['stage','depth', 'elevation', 'xmomentum', 'ymomentum'],
+                   verbose=True,
+                   use_cache = True)
+
+#csv2timeseries_graphs(directories_dic={anuga_dir+sep+'boundaries'+sep:['wavetank',0, 0],
+#                                       anuga_dir+'outputs'+sep+'20080219_233611_res_0.05_nbartzis'+sep:['Fixed Wave',0,0]},
+#                            output_dir=anuga_dir+'outputs'+sep+'20080219_233611_res_0.05_nbartzis'+sep,
+#                            base_name='gauge_',
+#                            plot_numbers='',
+#                            quantities=['stage'],
+#                            extra_plot_name='',
+#                            assess_all_csv_files=True,                            
+#                            create_latex=False,
+#                            verbose=True)
+
+
+
+
+
