source: trunk/anuga_validation/okushiri_2005/run_okushiri.py @ 8831

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

This script sets up Okushiri Island benchmark as published at the

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/


The validation data was downloaded and made available in this directory
for convenience but the original 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.py to process the boundary condition and build a the
mesh before running this script.

"""

# Module imports
from anuga import Domain
from anuga import Reflective_boundary
from anuga import Time_boundary
from anuga import Transmissive_momentum_set_stage_boundary
from anuga import Transmissive_n_momentum_zero_t_momentum_set_stage_boundary
from anuga import file_function

import project


#-------------------------
# Create Domain from mesh
#-------------------------
domain = Domain(project.mesh_filename, use_cache=True, verbose=True)
print domain.statistics()


#-------------------------
# Initial Conditions
#-------------------------
domain.set_quantity('friction', 0.0)
domain.set_quantity('stage', 0.0)
domain.set_quantity('elevation',
                    filename=project.bathymetry_filename,
                    alpha=0.02,                    
                    verbose=True,
                    use_cache=True)


#-------------------------
# Set simulation parameters
#-------------------------
domain.set_name(project.output_filename)  # Name of output sww file
domain.set_default_order(2)               # Apply second order scheme 
domain.set_minimum_storable_height(0.001) # Don't store h < 0.001m

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

# 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([4.521, 1.196])
# This should get triangle id 32833 with centroid (4.5244, 1.1972)
print 'triangle_id = ',triangle_id
bdry_id = domain.get_triangle_containing_point([0.000, 1.696])
# This should get triangle id 56474
print 'bdry_id = ',bdry_id


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

for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5):
    domain.write_time()
    print '   ',function(domain.get_time())[0],' bdry_d = ',bdry_id,' ',\
          domain.get_conserved_quantities(bdry_id, edge=0)[0],' ',\
          domain.get_conserved_quantities(bdry_id, edge=1)[0],' ',\
          domain.get_conserved_quantities(bdry_id, edge=2)[0]    
    

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