source: anuga_work/development/cairns_2006/run_cairns.py @ 4605

"""
Author: John Jakeman
Created: 12/12/2005
Main Cairns script using new interface
"""

#-------------------------------
# Module imports
#-------------------------------
import sys, os
from anuga.shallow_water import Domain, Reflective_boundary,\
     File_boundary, Transmissive_Momentum_Set_Stage_boundary,\
     Transmissive_boundary, Dirichlet_boundary
#from anuga.pyvolution.mesh_factory import rectangular_cross
#from anuga.pmesh2domain import pmesh_to_domain_instance
from anuga.pmesh.mesh_interface import create_mesh_from_regions
from Numeric import array, zeros, Float, allclose
import cairns_project as project
#import project
from caching import cache
from anuga.utilities.polygon import read_polygon, Polygon_function

from conversion import convert_latlon_to_xycoords

#-------------------------------
# Domain
#-------------------------------
print 'Creating domain'

"""
M = 36 #float(len1)/m
N = 10 #float(len2)/n

points, elements, boundary = rectangular_cross(M, N, len1=6393693.9504678631,
                                               len2=1781111.2941870166,
                                               origin = (0.0, 0.0))

#points, elements, boundary = rectangular_cross(N, M, len1=58.0, len2=1.0, origin=(145.0,-24.0))
                                               
domain = cache(Domain, (points, elements, boundary))

"""

#cache(pmesh_to_domain_instance,
#      (project.mesh_filename, Domain),
#      dependencies = [project.mesh_filename])
domain = Domain(project.mesh_filename, use_cache = True, verbose = True)

domain.check_integrity()
print 'Number of triangles = ', len(domain)
print 'The extent is ', domain.get_extent()

#---------------------------------------------------------
#Decide which quantities are to be stored at each timestep
domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']

#-------------------------------
# Initial Conditions
#-------------------------------
print 'Initial values'

latlong_origin = 145, -24
r0 = convert_latlon_to_xycoords(-173,-20,latlong_origin)
r1 = convert_latlon_to_xycoords(-173,-15,latlong_origin)
r2 = convert_latlon_to_xycoords(-171,-15,latlong_origin)
r3 = convert_latlon_to_xycoords(-171,-20,latlong_origin)

#p0 = read_polygon('cairns_fault.xya')
#p0 = read_polygon('cairns_fault_degrees.xya')
p0 = [r0, r1, r2, r3]
domain.set_quantity('stage',Polygon_function([(p0,1.0)]))

#print domain.quantities['stage'].vertex_values

domain.set_quantity('elevation',
                    filename = project.bathymetry_filename[:-4] + '.xya',
                    alpha = 10.0,
                    verbose = True,
                    use_cache = True)

# Need to increase elevation by the same amount as the stage
# At present unsure how to achieve this
#domain.set_quantity('elevation',Polygon_function([domain.quantities['elevation'].vertex_values+1000.0]))

domain.set_quantity('friction', 0.0)

#-------------------------------
# Boundary conditions
#-------------------------------
print 'Boundaries'

#Bts = Transmissive_Momentum_Set_Stage_boundary(domain, tide_function)
#Bd = Dirichlet_boundary([0,0,0])
Br = Reflective_boundary(domain)
domain.set_boundary({'border': Br, 'border': Br, 'border': Br, 'border': Br})

#Bt = Transmissive_boundary(domain)
#domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})

#-------------------------------
# Setup domain runtime parameters
#-------------------------------

#domain.visualise = True
#domain.visualise_color_stage = True

#Set domain.visualise = False for large problem sizes
domain.visualise = False
domain.visualise_color_stage = False

base = os.path.basename(sys.argv[0])
basename, _ = os.path.splitext(base)
domain.set_name(basename)
domain.default_order = 2
domain.store = True    #Store for visualisation purposes

#-------------------------------
# Evolve
#-------------------------------
import time
t0 = time.time()
yieldstep = 60.0
finaltime = 21600.0

for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
    domain.write_time()

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