source: inundation/validation/LWRU2/lwru2.py @ 1744

"""Example of 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
"""


read_mesh = False

import sys
from os import sep
sys.path.append('..'+sep+'..'+sep)

def prepare_timeboundary(filename):
    """Converting benchmark 2 time series to NetCDF sww file.
    """


    print 'Preparing time boundary from %s' %filename
    from Numeric import array

    fid = open(filename)

    #Skip first line
    line = fid.readline()

    lines = fid.readlines()
    fid.close()

    N = len(lines)

    T = zeros(N, Float)  #Time
    Q = zeros(N, Float)  #Values

    for i, line in enumerate(lines):
        fields = line.split()

        T[i] = float(fields[0])
        Q[i] = float(fields[1])

    #Create sww file
    from Scientific.IO.NetCDF import NetCDFFile

    outfile = filename[:-4] + '.sww'
    print 'Writing to', outfile
    fid = NetCDFFile(outfile, 'w')

    fid.institution = 'Geoscience Australia'
    fid.description = 'Input wave for Benchmark 2'
    fid.starttime = 0.0
    fid.createDimension('number_of_timesteps', len(T))
    fid.createVariable('time', Float, ('number_of_timesteps',))
    fid.variables['time'][:] = T

    fid.createVariable('stage', Float, ('number_of_timesteps',))
    fid.variables['stage'][:] = Q[:]

    fid.createVariable('xmomentum', Float, ('number_of_timesteps',))
    fid.variables['xmomentum'][:] = 0.0

    fid.createVariable('ymomentum', Float, ('number_of_timesteps',))
    fid.variables['ymomentum'][:] = 0.0

    fid.close()


# Module imports
from pyvolution.shallow_water import Domain, Reflective_boundary,\
            File_boundary, Transmissive_Momentum_Set_Stage_boundary
from pyvolution.mesh_factory import rectangular
from pyvolution.pmesh2domain import pmesh_to_domain_instance
from Numeric import array, zeros, Float, allclose
import filenames
from caching import cache

#Preparing time boundary
prepare_timeboundary(filenames.boundary_filename)

#Preparing points
from pyvolution.data_manager import xya2pts
xya2pts(filenames.bathymetry_filename, verbose = True,
        z_func = lambda z: -z)


#######################
# Domain
    
if read_mesh is True:
    print 'Creating domain from', filenames.mesh_filename
    domain = pmesh_to_domain_instance(filenames.mesh_filename, Domain)    
else:   
    print 'Creating regular mesh'
    N = 100
    points, vertices, boundary = rectangular(N, N/5*3,
                                             len1=5.448, len2=3.402)
    print 'Creating domain'
    
    #domain = Domain(points, vertices, boundary)
    domain = cache(Domain, (points, vertices, boundary))



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


#domain.check_integrity()
domain.default_order = 2

print "Number of triangles = ", len(domain)
domain.store = True    #Store for visualisation purposes

import sys, os
base = os.path.basename(sys.argv[0])
domain.filename, _ = os.path.splitext(base)

#LS code to be included in set_quantity
print 'Reading points'
from pyvolution import util, least_squares
points, attributes = util.read_xya(filenames.bathymetry_filename[:-4] + '.pts')
elevation = attributes['elevation']

#Fit attributes to mesh
#vertex_attributes = least_squares.fit_to_mesh(domain.coordinates,
#                                              domain.triangles,
#                                              points,
#                                              attributes['elevation'],
#                                              alpha = 0.0001,
#                                              verbose=True)

vertex_attributes = cache(least_squares.fit_to_mesh,
                          (domain.coordinates, domain.triangles,
                           points, elevation) ,
                          {'alpha': 0.0001, 'verbose': True})

print 'Initial values'
domain.set_quantity('elevation', vertex_attributes)
#domain.set_quantity('elevation', points=points, attributes=attributes['elevation'], alpha=0.0001)
domain.set_quantity('friction', 0.0)
domain.set_quantity('stage', 0.0)



######################
# Boundary conditions
#
print 'Boundaries'
Br = Reflective_boundary(domain)

#####################
# Error with File Boundary. We need to only set
# stage and let the scheme set the momentum
from pyvolution.util import file_function
function = file_function(filenames.boundary_filename[:-4] + '.sww', domain,
                         verbose = True)
Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)

#Set boundary conditions
if read_mesh is True:
    domain.set_boundary({'wave': Bts, 'wall': Br})      
else:    
    domain.set_boundary({'left': Bts, 'right': Br, 'bottom': Br, 'top': Br})


#Evolve
import time
t0 = time.time()
#domain.visualise = True

for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5):
    domain.write_time()
    #Stage_cv = domain.quantities['stage'].centroid_values
    #print Stage_cv[0]
    #print 'Function value',function(t)

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