source: trunk/anuga_work/development/2010-projects/anuga_1d/sww/1dprofile_dispersion.py @ 8252

import os
from math import *
import random
import numpy as np
import time
import anuga
import csv
from Scientific.IO.NetCDF import NetCDFFile
from anuga_1d.sww.sww_domain import *
from anuga_1d.config import g, epsilon
from anuga_1d.base.generic_mesh import uniform_mesh


#===========================================================================


# Define functions for initial quantities
## initial area for 1D SWE
def initial_area(x):
    return (stage(x)-sea_bed(x))*width(x)

#--------------------------------------------
def width(x):
    return 1.0
#---------------------------------------------
# sea bed or bathymetry data

def elevation():
    x=[]
    elevation=[]

    f=open('gebco_1m.txt')
    lines=f.readlines()
    f.close
    for line in lines[2:]:
        val=[float(v) for v in line.split()]
        y=val[1]
        r=y-4152074.2
        if abs(r)<810:# and val[0]<=max(x):
    #if np.allclose(y,4243929.72563): 
            x.append(val[0])
            elevation.append(val[2])
    return x, elevation
#-----------------------------------------------------
# initial stage 
# Formula for solitary wave in 1D 
 
def stage(x):
    a=0.0*x
    #a=gaussian_filtered(x)
    y = 0.0*x
    wd=abs(-6000)  # water depth
    wh=2.5    # wave height
    xm=908432.19999999995
    k=10*sqrt(3.0*wh/(4.0*wd))/wd
    for i in range(len(x)):
            th=tanh(k*(x[i]-xm))**2
            y[i]= wh*(1.0-th)#+a[i]
    return y
# --------------------------------------------------------


# Initial momentum is zero

def initial_momentum(x):
    return 1.0


# Length of channel (m)

# Define cells for finite volume and their size
# Define the number of cells
x,elev=elevation()
print x
boundary = { (0,0) : 'left',  (len(x)-2, 1) : 'right'} 
print boundary
domain = Domain(x,boundary)
print domain.boundary
print type(x), len(x), type(elev),len(elev)


# Set initial values of quantities - default to zero
domain.set_quantity('elevation', elev, location='unique_vertices')
#domain.add_quantity('stage',stage)
domain.set_quantity('stage',stage)
domain.set_quantity('xmomentum',initial_momentum)

# Set boundry type, order, timestepping method and limiter
print 'Available boundary tags: ', domain.get_boundary_tags()

Bd = Dirichlet_boundary([0, 0, 0, 0, 0]) #(w,uh,z,h,u)
Br = Reflective_boundary(domain)
Bt = Transmissive_boundary(domain)
domain.set_boundary({'left': Br , 'right' : Br})
#domain.set_boundary({'exterior': Br})
print 'Available boundary tags: ', domain.get_boundary_tags()
domain.order = 2
domain.set_timestepping_method(1)
domain.set_CFL(1.0)
domain.set_limiter("vanleer")
#domain.h0=0.0001


# Start timer
t0 = time.time()

# Set final time and yield time for simulation
finaltime = 1.0
yieldstep = 1.0

#===================================================================
# Time loop
#===================================================================
for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):

    w  = domain.quantities['stage'].vertex_values

    
    
    domain.write_time()


import pylab
pylab.ion()

x = domain.get_vertices().flatten()

z  = domain.quantities['elevation'].vertex_values.flatten()
w  = domain.quantities['stage'].vertex_values.flatten()
h  = domain.quantities['height'].vertex_values.flatten()
u  = domain.quantities['velocity'].vertex_values.flatten()
uh = domain.quantities['xmomentum'].vertex_values.flatten()

print x.shape
print z.shape
print w.shape
print u.shape


#print w
#print z
#print h
#print u
#print uh
#print x

#-------------------------------
# Top plot
#-------------------------------
plot1 = pylab.subplot(211)

zplot = pylab.plot(x, z, 'k')
wplot = pylab.plot(x, w, 'k')

#plot1.set_ylim(stage_lim)
#pylab.xlabel('Position')
pylab.ylabel('Bed/Stage')

#-------------------------------
# Bottom Plot
#-------------------------------
plot2 = pylab.subplot(212)

uplot = pylab.plot(x, u, 'k')
pylab.draw()

#plot2.set_ylim(width_lim)
#pylab.xlabel('Position')
#pylab.ylabel('Velocity')


#pylab.savefig('figure4.pdf')
pylab.show()