source: anuga_validation/okushiri_2005/okushiri_parallel.py @ 3640

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

This a parallel version.
"""

# 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

from anuga_parallel.parallel_api import myid, numprocs, distribute

import project



def prepare_timeboundary(filename):
    """Converting benchmark 2 time series to NetCDF tms file.
    This is a 'throw-away' code taylor made for files like
    'Benchmark_2_input.txt' from the LWRU2 benchmark
    """

    print 'Preparing time boundary from %s' %filename
    fid = open(filename)

    #Skip first line
    line = fid.readline()

    #Read remaining lines
    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 tms file
    from Scientific.IO.NetCDF import NetCDFFile

    outfile = filename[:-4] + '.tms'
    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()



#----------------------

print 'That took %.2f seconds' %(time.time()-t0)---
# Preparing time boundary
#-------------------------
if myid == 0:
    prepare_timeboundary(project.boundary_filename)


#-------------------------
# Create Domain
#-------------------------
domain = Domain(project.mesh_filename, use_cache=True, verbose=True)
domain.set_name('okushiri')
domain.set_default_order(2)
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,                    
                    verbose = True,
                    use_cache = True)

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

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

# Note this boundary can be fully specified here and
# distributed without problem.
Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)

domain.set_boundary({'wave': Bts, 'wall': Br})


#-------------------------
# Distribute domain
#-------------------------
domain = distribute(domain)


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

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

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