"""Anuga convergence study using a true-scale version of the Okushiri
Island tsunami wavetank experiment

This script runs a modified version of the Okushiri Island benchmark

as published at

THE THIRD INTERNATIONAL WORKSHOP ON LONG-WAVE RUNUP MODELS
June 17-18 2004
Wrigley Marine Science Center
Catalina Island, California
http://www.cee.cornell.edu/longwave/

This version up-scales the original 1:400 scale wave-tank experiment, to
"true-scale" with mesh defined using interior polygons.

The original validation data is available at
http://www.cee.cornell.edu/longwave/index.cfm?page=benchmark&problem=2
where a detailed description of the problem is also available.

Run create_okushiri_truescale.py to convert bathymetry and input boundary
condition into NetCDF format.

"""
#----------------------------------------
# Import necessary modules
#----------------------------------------

# Standard modules
from os import sep,umask
from os.path import dirname, basename
from os import mkdir, access, F_OK
from shutil import copy
import time
import sys

# Related major packages
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.abstract_2d_finite_volumes.util import file_function
from anuga.shallow_water.data_manager import copy_code_files, start_screen_catcher
from anuga.pmesh.mesh_interface import create_mesh_from_regions


import project_truescale # Definition of filenames and interior polygons

copy_code_files(project_truescale.output_run_time_dir,'run_okushiri_truescale.py',
                'project_truescale.py' )
myid = 0
numprocs = 1
start_screen_catcher(project_truescale.output_run_time_dir, myid, numprocs)

#--------------------------------------------------------------------------
# Create the triangular mesh based on overall bounding polygon with a
# tagged boundary and interior regions defined in project_truescale.py
# along with resolutions (maximal area of per triangle) for each polygon
#--------------------------------------------------------------------------

create_mesh_from_regions(project_truescale.poly_all,
                         boundary_tags={'wall': [0, 1, 3],'wave': [2]},
                         maximum_triangle_area=project_truescale.res_poly_all,
                         #interior_regions=project_truescale.interior_regions,
                         filename=project_truescale.mesh_name+'.msh',
                         verbose=True)

#------------------------------
# Create Domain from mesh
#------------------------------
domain = Domain(project_truescale.mesh_name+'.msh', use_cache=True, verbose=True)
print domain.statistics()

z = domain.get_vertex_coordinates()

# Write vertex coordinates to file
filename=project.vertex_filename
fid=open(filename,'w')
fid.write('x (m), y (m)\n')
for i in range(len(z)):
    pt=z[i]
    x=pt[0]
    y=pt[1]
    fid.write('%.6f, %.6f\n' %(x, y))
        
#-------------------------
# Initial Conditions
#-------------------------
domain.set_quantity('friction', 0.0)
domain.set_quantity('stage', 0.0)
domain.set_quantity('elevation',
                    filename=project_truescale.bathymetry_filename,
                    alpha=0.02,                    
                    verbose=True,
                    use_cache=True)


#-------------------------
# Set simulation parameters
#-------------------------
domain.set_name(project_truescale.output_filename)  # Name of output sww file
domain.set_datadir(project_truescale.output_run_time_dir) # Name of output directory
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.1)   # Don't store h < 0.1m
domain.tight_slope_limiters = 1
domain.beta_h = 0.0


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

# Create boundary function from timeseries provided in file
function = file_function(project_truescale.boundary_filename,
                         domain, verbose=True)

# Create and assign boundary objects
Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
Br = Reflective_boundary(domain)
domain.set_boundary({'wave': Bts, 'wall': Br})


# Select triangle containing ch5 for diagnostic output
# around known gauge
triangle_id = domain.get_triangle_containing_point([1808.4, 478.4])
# This should get triangle id 32833 with centroid (4.5244, 1.1972)


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

for t in domain.evolve(yieldstep=project_truescale.yieldstep, finaltime = 450):
    print domain.timestepping_statistics(track_speeds=False,
                                         triangle_id=triangle_id)
    print domain.boundary_statistics()

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