source: anuga_validation/conical_island/run_circular.py @ 4810

"""Validation of the AnuGA implementation of the shallow water wave equation.

This script sets up 2D circular island.....

"""

# 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
from anuga.abstract_2d_finite_volumes.util import file_function
from anuga.pmesh.mesh_interface import create_mesh_from_regions
from anuga.utilities.polygon import read_polygon, plot_polygons
#from cmath import cos,sin,cosh,pi
from math import cos,pi,sin,tan,sqrt
from Numeric import array, zeros, Float, allclose

#import project

#-------------------------
# Create Domain from mesh
#-------------------------

def get_xy(x=0,y=0,r=1,angle=1):
    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()



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)
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)

x = 12.96
y = 13.80
r1 = 1.5
r2 = 2.5
d=0.75
res=0.0005

poly = []
poly1 = []
internal_poly=[]
for i,angle in enumerate(angles1):
    #convert Degs to Rads
    angle = ((angle-180)/180.0)*pi
#    p1 = get_xy(x,y,r1,angle-0.01745*d)
#    p2 = get_xy(x,y,r2,angle-0.01745*d)
#    p3 = get_xy(x,y,r2,angle+0.01745*d)
#    p4 = get_xy(x,y,r1,angle+0.01745*d)
#    poly[i] = [[p1[0],p1[1]],[p2[0],p2[1]],[p3[0],p3[1]],[p4[0],p4[1]]]

    p = get_xy(x,y,3.6,angle)
    poly1.append([p[0],p[1]])

poly.append(poly1)
internal_poly.append([poly1,.1])

poly2 = []
for i,angle in enumerate(angles1):
    #convert Degs to Rads
    angle = ((angle-180)/180.0)*pi
    p = get_xy(x,y,2.4,angle)
    poly2.append([p[0],p[1]])
    
poly.append(poly2)
internal_poly.append([poly2,.01])

poly3 = []
for i,angle in enumerate(angles1):
    #convert Degs to Rads
    angle = ((angle-180)/180.0)*pi
    p = get_xy(x,y,1.6,angle)
    poly3.append([p[0],p[1]])
    
poly.append(poly3)
internal_poly.append([poly3,1])

poly4 = []
for i,angle in enumerate(angles1):
    #convert Degs to Rads
    angle = ((angle-180)/180.0)*pi
    p = get_xy(x,y,1.1,angle)
    poly4.append([p[0],p[1]])
    
poly.append(poly4)
internal_poly.append([poly2,.01])

plot_polygons(poly,'poly.png')

poly_all= [[0,7],[0,25],[30,25],[30,7]]
create_mesh_from_regions(poly_all,
                         boundary_tags={'wall': [0,2],'wave': [3], 'buffer':[1]},
                         maximum_triangle_area=1,
                         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
domain.set_quantity('stage', water_height)

def test():                
    fileName = "test.csv"
    file = open(fileName,"w")
    file.write("x,y,elevation \n")
    x = 12.96
    y = 13.80
    for i in range(0,30):
        for j in range(0,25):
            a = (x-i)**2
            b = (y-j)**2
            z = -sqrt((x-i)**2+(y-j)**2)
            z = z*.25+(0.8975)
            if z <0:
                z=0
            if z > .625:
                z=0.625

#            print 'x,y,f',i,j,a,b,z
            file.write("%s, %s, %s\n" %(i, j, z))
    file.close()

test()

domain.set_quantity('elevation',filename = "test.csv", alpha=0.1)


#-------------------------
# Set simulation parameters
#-------------------------
domain.set_name('test1')  # Name of output sww file
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.01) # 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

#Timings on AMD64-242 (beta_h=0)
#  tight_slope_limiters = 0:
#    3035s - 3110s
#  tight_slope_limiters = 1:
#    3000s - 3008s
#
# beta_h>0: In the order of 3200s

#-------------------------
# Boundary Conditions
#-------------------------

# Create boundary function from timeseries provided in file
#function = file_function(project.boundary_filename,
#                         domain, verbose=True)
def waveform(t): 
        return (sin(t*2*pi))

# Create and assign boundary objects
boundary_filename="ts2cnew1_input"
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)

#Bts = Transmissive_Momentum_Set_Stage_boundary(domain, waveform)
Bw = Time_boundary(domain=domain,
                   f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
Br = Reflective_boundary(domain)
Bd = Dirichlet_boundary([water_height,0,0])
domain.set_boundary({'wave': Bts, 'wall': Br, 'buffer':Bd})

#-------------------------
# Evolve through time
#-------------------------
import time
t0 = time.time()

for t in domain.evolve(yieldstep = 1, finaltime = 32):
    domain.write_time()
#    domain.write_time(track_speeds=False)

for t in domain.evolve(yieldstep = 0.1, finaltime = 36
                       ,skip_initial_step = True): 
    domain.write_time()
   
for t in domain.evolve(yieldstep = 1, finaltime = 75
                       ,skip_initial_step = True): 
    domain.write_time()

print 'That took %.2f seconds' %(time.time()-t0)



#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)
#x = 12.96
#y = 13.80
#r1 = 1.1
#r2 = 3.6
#d=2
#poly = {}
for i,angle in enumerate(angles):
#    angle = ((angle-180)/180.0)*pi
##    angle = angle1*3.14157
#    p1 = get_xy(x,y,r1,angle-0.01745*d)
#    p2 = get_xy(x,y,r2,angle-0.01745*d)
#    p3 = get_xy(x,y,r2,angle+0.01745*d)
#    p4 = get_xy(x,y,r1,angle+0.01745*d)
#    poly[i] = [[p1[0],p1[1]],[p2[0],p2[1]],[p3[0],p3[1]],[p4[0],p4[1]]]
##    print i,poly[i]
#    
    
    run_up, x_y = get_maximum_inundation_data(filename='test.sww',polygon=poly[i], verbose=False)
    
    print 'maximum_inundation_data',angle, run_up, x_y