source: trunk/anuga_validation/circular_island/run_circular.py @ 8831

"""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.shallow_water_domain 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.utilities.file_utils import copy_code_files
from anuga.shallow_water.sww_interrogate import get_maximum_inundation_data
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
from math import cos,pi,sin,tan
import numpy as num
from anuga.shallow_water.data_manager import load_csv_as_dict as 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
from anuga.config import netcdf_float
from Scientific.IO.NetCDF import NetCDFFile

#-------------------------
# 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
#output_dir = anuga_dir+'outputs'+sep+'20080222_062917_res_0.0025_nbartzis'+sep
copy_code_files(output_dir,__file__)

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

    fid = open(textversion)

    #Skip first line
    line = fid.readline()

    #Read remaining lines
    lines = fid.readlines()
    fid.close()


    N = len(lines)
    T = num.zeros(N, num.float)  #Time
    Q = num.zeros(N, num.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
    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', netcdf_float, ('number_of_timesteps',))
    fid.variables['time'][:] = T

    fid.createVariable('stage', netcdf_float, ('number_of_timesteps',))
    fid.variables['stage'][:] = Q[:]

    fid.createVariable('xmomentum', netcdf_float, ('number_of_timesteps',))
    fid.variables['xmomentum'][:] = 0.0

    fid.createVariable('ymomentum', netcdf_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

#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.5,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,3.2,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)
domain.set_quantity('friction', 0.01)#for concrete
#water_height = 0.322
water_height = 0.0

domain.set_quantity('stage', water_height)


#attribute_dic, title_index_dic = csv2dict('sqrt_table.csv')

def circular_island_elevation(x,y):
    water_depth = 0.32
    list_xy = num.sqrt((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:
    for distance in list_xy:

#    The cone has a 1/4 slope and a 7.2 diameter this make the height 0.9m
#    to get the stage at 0 elevation + cone height and - water depth
        z = -distance*.25+0.9-water_depth

        if z <0-water_depth:
            z=0-water_depth
        if z > .625-water_depth:
            z=0.625-water_depth
        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(2)               # 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

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]]]
    
    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=anuga_dir+'gauges'+sep+'gauges.csv',
#                   gauge_file='gauges.csv',
                   quantities = ['stage','depth', 'elevation', 'xmomentum', 'ymomentum'],
                   verbose=True,
                   use_cache = True)