source: anuga_work/development/analytical_solutions/contracting_channel.py @ 5162

"""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
from mesh_factory import contracting_channel_cross

#-------
# Domain from a file
# filename = 'converging_channel_30846.tsh'
# print 'Creating domain from', filename
# domain = pmesh_to_domain_instance(filename, Domain)

######################
# Domain created within python
#
Total_length = 50
W_upstream = 5.
W_downstream = 2.5
L_1 = 5.
L_2 = 11
L_3 = Total_length - L_1 - L_2
n = 5
m = 50

points, elements, boundary = \
     contracting_channel_cross(m, n, W_upstream, W_downstream, L_1, L_2, L_3)
domain = Domain(points, elements, boundary)


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']

#------------------------------------------------
# This is for Visual Python
domain.visualise = True

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

#------------------------
# Set boundary Conditions
tags = {}
tags['left']   = Dirichlet_boundary([0.2, 1.2, 0.0])
tags['top']    = Reflective_boundary(domain)
tags['bottom'] = Reflective_boundary(domain)
tags['right']  = Transmissive_boundary(domain)
domain.set_boundary(tags)

#----------------------
# Set initial condition
domain.set_quantity('elevation', 0.0)
domain.set_quantity('stage', 0.2)

# Use the inscribed circle with safety factor of 0.9 to establish the time step
# domain.set_to_inscribed_circle(safety_factor=0.9)

#----------
# Evolution
import time
t0 = time.time()
for t in domain.evolve(yieldstep = 5.0, finaltime = 50.0):
    domain.write_time()

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

N = domain.number_of_elements

Stage = domain.quantities['stage']
stage = Stage.centroid_values
XY = domain.centroid_coordinates

# Calculate average
average_stage  = 0.0
n_points = 0
for n in range(N):
    if XY[n,0] > 35.0:
        average_stage = average_stage + stage[n]
        n_points = n_points + 1

average_stage = average_stage/n_points

#Standard Deviation
sigma = 0.0
max_stage = -999999.
min_stage = 999999
for n in range(N):
    if XY[n,0] > 35.0:
        sigma = sigma + (average_stage - stage[n])**2
        if stage[n] > max_stage:
            max_stage = stage[n]
        if stage[n] < min_stage:
            min_stage = stage[n]

import math
sigma = math.sqrt(sigma/(n_points-1))

print L_2, average_stage, sigma, max_stage, min_stage, n_points

domain.initialise_visualiser()
domain.visualiser.update_all()