"""Example of shallow water wave equation analytical solution
consists of a symmetrical converging frictionless channel.

Specific methods pertaining to the 2D shallow water equation
are imported from shallow_water
for use with the generic finite volume framework

   Copyright 2005
   Christopher Zoppou, Stephen Roberts
   ANU
   
Specific methods pertaining to the 2D shallow water equation
are imported from shallow_water
for use with the generic finite volume framework

Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the
numerical vector named conserved_quantities.
"""

#---------------
# Module imports 
import sys
from os import sep
sys.path.append('..'+sep+'pyvolution')

from shallow_water import Transmissive_boundary, Reflective_boundary, \
     Dirichlet_boundary
from shallow_water import Constant_height, Domain
from pmesh2domain import pmesh_to_domain_instance

#-------
# Domain 
filename = 'converging_channel_30846.tsh'
print 'Creating domain from', filename
domain = pmesh_to_domain_instance(filename, Domain)
print 'Number of triangles = ', len(domain)

#----------------
# Order of scheme
domain.default_order = 2
domain.smooth = True

#-------------------------------------
# Provide file name for storing output
domain.store = True     #Store for visualisation purposes
domain.format = 'sww'   #Native netcdf visualisation format
domain.filename = 'contracting_channel_second-order'

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

#------------------------------------------
# Reduction operation for get_vertex_values              
#from util import mean
#domain.reduction = mean 

#------------------------
# Set boundary Conditions
tags = {}
tags['upstream'] = Dirichlet_boundary([0.2, 1.2, 0.0])
tags['reflective'] = Reflective_boundary(domain) 
tags['transmissive'] = Transmissive_boundary(domain)
domain.set_boundary(tags)

#----------------------
# Set initial condition
domain.set_quantity('elevation', 0.0)
domain.set_quantity('stage', 0.2)
    
#----------
# Evolution
import time
t0 = time.time()
for t in domain.evolve(yieldstep = 0.1, finaltime = .2):
    domain.write_time()
    
print 'That took %.2f seconds' %(time.time()-t0)