source: anuga_validation/circular_island/simple_circular/circular.py @ 5317

"""Simple water flow example using ANUGA

Water flowing down a channel with more complex topography
"""

#------------------------------------------------------------------------------
# Import necessary modules
#------------------------------------------------------------------------------
from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
from anuga.shallow_water import Domain
from anuga.shallow_water import Reflective_boundary
from anuga.shallow_water import Dirichlet_boundary
from anuga.shallow_water import Time_boundary
from anuga.pmesh.mesh_interface import create_mesh_from_regions
from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
from math import tan, sqrt, sin, pi

#------------------------------------------------------------------------------
# Setup computational domain
#------------------------------------------------------------------------------
length = 30.
width = 25.
Cx = 12.96                  # centre of island on the x axis
Cy = 13.8                   # centre of island on the y axis
dx = dy = 1                 # Resolution: Length of subdivisions on both axes
water_depth = 0.32          # Can be 0.32 or 0.42

#boundary
poly_domain = [[0,0],[length,0],[length,width],[0,width]]
meshname = 'test.msh'

# exporting asc grid
xmin = 0
xmax = length   
ymin = 0
ymax = width

#Create interior region
Dome = [[(Cx)-4,(Cy)-4],[(Cx)+4,(Cy)-4,],
        [(Cx)+4,(Cy)+4],[(Cx)-4,(Cy)+4]]

remainder_res=1
Dome_res = .01

interior_dome = [[Dome, Dome_res]]

#Create mesh
create_mesh_from_regions(poly_domain,
                         boundary_tags={'wavemaker': [0], 'right': [1],
                                        'top': [2], 'left': [3]},
                         maximum_triangle_area = remainder_res,
                         filename=meshname,
                         interior_regions = interior_dome,
                         use_cache=False,
                         verbose=True)

domain = Domain(meshname, use_cache=False, verbose = True)
domain.set_name('circular')                  # Output name
print domain.statistics()

#------------------------------------------------------------------------------
# Setup initial conditions
#------------------------------------------------------------------------------
def topography(x,y):
    """Complex topography defined by a function of vectors x and y
    """
    
                                    
    z= 0*x                      # defining z for all values other than the if statements
    r= 3.6                      # radius, provided in document
    angle = 14                  # angle, provided in document
    h= r*tan(angle/57.2957795)  # finding height of cone if not truncated
    

    N = len(x)
    for i in range(N):
       
        #truncated top
        if (x[i]-Cx)**2 + (y[i]-Cy)**2 <1.1**2:
            z[i] += 0.625

        # cone
        if (x[i]-Cx)**2 + (y[i]-Cy)**2 <r**2 and (x[i]-Cx)**2 + (y[i]-Cy)**2 >1.1**2:
           z[i] = -(sqrt(((x[i]-Cx)**2+(y[i]-Cy)**2)/((r/h)**2))-h)
    return z   

domain.set_quantity('elevation', topography, verbose=True)  # Use function for elevation
domain.set_quantity('friction', 0.01)         # Constant friction 
domain.set_quantity('stage',water_depth)   # Dry initial condition


#------------------------------------------------------------------------------
# Setup boundary conditions
#------------------------------------------------------------------------------
# Create boundary function from timeseries provided in file

##boundary_filename="ts2cnew1_input_20_80sec_new"
##prepare_timeboundary(boundary_filename+'.txt')
##
##function = file_function(boundary_filename+'.tms',
##                         domain, verbose=True)
def wave_form(t):
    return 0.1*sin(2*pi*t/50.)

# Create and assign boundary objects
Bw = Dirichlet_boundary([water_depth, 0, 0])          #wall
Bt = Transmissive_Momentum_Set_Stage_boundary(domain, wave_form) #wavemaker
domain.set_boundary({'left': Bw, 'right': Bw, 'top': Bw, 'wavemaker': Bt})




#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------
for t in domain.evolve(yieldstep = 0.2, finaltime = 100.0):
    print domain.timestepping_statistics()