source: trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py @ 9016

"""Runup example from the manual, slightly modified
"""
#---------
#Import Modules
#--------
import anuga

import numpy

from math import sin, pi, exp
#from anuga.shallow_water.shallow_water_domain import Domain as Domain
#from anuga.shallow_water_balanced2.swb2_domain import Domain as Domain
#path.append('/home/gareth/storage/anuga_clean/anuga_jan12/trunk/anuga_work/development/gareth/balanced_basic')
#from swb2_domain import *
#from balanced_basic import *
#from balanced_dev import *
from bal_and import *
#from anuga_tsunami import *

#---------
#Setup computational domain
#---------
points, vertices, boundary = anuga.rectangular_cross(20,20, len1=1., len2=1.)

domain=Domain(points,vertices,boundary)    # Create Domain
domain.set_name('runup_sinusoid_v2')                         # Output to file runup.sww
#domain.set_timestepping_method('euler')
#domain.set_flow_algorithm('tsunami')
#------------------
# Define topography
#------------------

#### Pathological
#scale_me=100.0
#boundary_ws=-0.2#1999
#init_ws=-0.2
#bumpiness=50. # Higher = shorter wavelength oscillations in topography
#tstep=0.2
#lasttime=20.1

### Sensible
scale_me=1.0
boundary_ws=-0.1
init_ws=-0.2
bumpiness=50. # Higher = shorter wavelength oscillations in topography
tstep=0.2
lasttime=150.

#domain.minimum_allowed_height=domain.minimum_allowed_height*scale_me # Seems needed to make the algorithms behave

def topography(x,y):
        return (-x/2.0 +0.05*numpy.sin((x+y)*bumpiness))*scale_me

def stagefun(x,y):
    stge=init_ws*scale_me# +0.1*(x>0.9)*scale_me 
    #topo=topography(x,y) 
    return stge#*(stge>topo) + (topo)*(stge<=topo)

domain.set_quantity('elevation',topography)     # Use function for elevation
domain.get_quantity('elevation').smooth_vertex_values() 
domain.set_quantity('friction',0.00)             # Constant friction

print '#-#', domain.minimum_allowed_height
#domain.minimum_allowed_height=0.001
#def frict_change(x,y):
#       return 0.2*(x>0.5)+0.1*(x<=0.5)
#
#domain.set_quantity('friction',frict_change)

domain.set_quantity('stage', stagefun)              # Constant negative initial stage
domain.get_quantity('stage').smooth_vertex_values()

#print domain.forcing_terms
#assert 0==1
# Experiment with rain.
# rainin = anuga.shallow_water.forcing.Rainfall(domain, rate=0.001) #, center=(0.,0.), radius=1000. )
# domain.forcing_terms.append(rainin)

#--------------------------
# Setup boundary conditions
#--------------------------
Br=anuga.Reflective_boundary(domain)            # Solid reflective wall
Bt=anuga.Transmissive_boundary(domain)          # Continue all values of boundary -- not used in this example
Bd=anuga.Dirichlet_boundary([boundary_ws*scale_me,0.,0.])       # Constant boundary values -- not used in this example
def waveform(t):
    return boundary_ws*scale_me
    #return -0.2*scale_me #-0.1 #(0.0*sin(t*2*pi)-0.1)*exp(-t)-0.1
def waveform2(t):
    return [ (0.1*sin(2*pi*t)-0.3)*exp(-t), 0., 0.]

Bt2=anuga.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain,waveform)
Bw=anuga.Time_boundary(domain=domain, waveform2) # Time varying boundary -- get rid of the 0.0 to do a runup.

#----------------------------------------------
# Associate boundary tags with boundary objects
#----------------------------------------------
domain.set_boundary({'left': Br, 'right': Bt2, 'top': Br, 'bottom':Br})

#------------------------------
#Evolve the system through time
#------------------------------

for t in domain.evolve(yieldstep=tstep,finaltime=lasttime):
#for t in domain.evolve(yieldstep=0.2,finaltime=60.):
    print domain.timestepping_statistics()
    #print domain.boundary_flux_integral
    xx = domain.quantities['xmomentum'].centroid_values
    yy = domain.quantities['ymomentum'].centroid_values
    dd = domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values 
    dd_raw=1.0*dd
    dd = dd*(dd>1.0e-03)+1.0e-03*(dd<=1.0e-03)
    vv = ( (xx/dd)**2 + (yy/dd)**2)**0.5
    vv = vv*(dd>1.0e-03)
    print 'Peak velocity is: ', vv.max(), vv.argmax(), dd[vv.argmax()]
    print 'Volume is', sum(dd_raw*domain.areas)    
    #print 'Volume less flux int', sum(dd_raw*domain.areas) - domain.boundary_flux_integral


vel_final_inds=(vv>1.0e-01).nonzero()
print 'vel_final_inds', vel_final_inds

print 'Finished'