source: anuga_validation/okushiri_2005/run_okushiri.py @ 3758

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

This script sets up LWRU2 benchmark with initial condition stated

See also

http://www.cee.cornell.edu/longwave/index.cfm?page=benchmark&problem=2

Depth at western boundary is d = 13.5 cm
"""

# Module imports
from Numeric import array, zeros, Float, allclose

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

import project


#-------------------------
# Create Domain
#-------------------------

use_variable_mesh = True #Use large variable mesh generated by create_mesh.py
#use_variable_mesh = False #Use small structured mesh for speed

if use_variable_mesh is True:
    print 'Creating domain from', project.mesh_filename

    domain = Domain(project.mesh_filename, use_cache=True, verbose=True)
else:
    print 'Creating regular from regular mesh'
    N = 150
    points, vertices, boundary = rectangular_cross(N, N/5*3,
                                                   len1=5.448, len2=3.402)
    domain = Domain(points, vertices, boundary)

domain.set_name('okushiri_old_limiters_smooth=0.03_maxA=0.00001')
domain.set_default_order(2)
domain.set_minimum_storable_height(0.001)

# Set old (pre Sep 2006) defaults for limiters
domain.beta_w      = 0.9
domain.beta_w_dry  = 0.9
domain.beta_uh     = 0.9
domain.beta_uh_dry = 0.9
domain.beta_vh     = 0.9
domain.beta_vh_dry = 0.9

domain.check_integrity()
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[:-4] + '.pts',
                    #alpha=0.02, #This was the best one with setMaxArea(0.00002*base_resolution)
                    alpha=0.03,                    
                    verbose=True,
                    use_cache=True)

#-------------------------
# Boundary Conditions
#-------------------------
Br = Reflective_boundary(domain)

function = file_function(project.boundary_filename[:-4] + '.tms',
                         domain,
                         verbose=True)

Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)

if use_variable_mesh is True:
    domain.set_boundary({'wave': Bts, 'wall': Br})
else:
    domain.set_boundary({'left': Bts, 'right': Br, 'bottom': Br, 'top': Br})


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

w_max = 0
for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5):
    domain.write_time()

    w = domain.get_maximum_inundation_elevation()
    x, y = domain.get_maximum_inundation_location()
    print '  Coastline elevation = %.2f at (x,y)=(%.2f, %.2f)' %(w, x, y)
    print
    
    if w > w_max:
        w_max = w
        x_max = x
        y_max = y


print '**********************************************'
print 'Max coastline elevation = %.2f at (%.2f, %.2f)' %(w_max, x_max, y_max)
print '**********************************************'
    
   

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