source: anuga_work/development/landslide_4/benchmark_problem4.py @ 5568

Last change on this file since 5568 was 5568, checked in by sexton, 14 years ago

incorporated definitions from web and Synolakis paper

File size: 4.9 KB
Line 
1"""Script for running the 3rd Long Wave Workshop Benchmark Problem 4.
2
3Bridgette Lewis and Jane Sexton, July 2008
4"""
5
6#-------------------------------------------------------------------------------
7# Import necessary modules
8#-------------------------------------------------------------------------------
9
10# Standard modules
11import os
12import time
13from shutil import copy
14from os.path import dirname, basename
15from os import mkdir, access, F_OK, sep
16import sys
17
18# Related major packages
19from anuga.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary
20from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
21
22#----------------------------------------------------------------------------
23# Create the triangular mesh based on the problem statement
24#-------------------------------------------------------------------------------
25
26from anuga.pmesh.mesh_interface import create_mesh_from_regions
27remainder_res = 250000.
28#local_res = 50000.
29#coast_res = 500.
30#interior_regions = [[project_slide.poly_syd1, local_res],
31#                    [project_slide.poly_coast, coast_res]]
32length = 104.
33width = 3.7
34
35bounding_polygon = [[0,0],[length,0],[length,width],[0,width]]
36meshname = 'slide.msh'
37remainder_res = 1.
38
39from caching import cache
40_ = cache(create_mesh_from_regions,
41          bounding_polygon,
42           {'boundary_tags': {'bottom': [0], 'right': [1], 
43                              'top': [2], 'left': [3]},
44           'maximum_triangle_area': remainder_res,
45           'filename': meshname},
46           #'interior_regions': interior_regions},
47          verbose = True, evaluate=False)
48print 'created mesh'
49
50#-------------------------------------------------------------------------------                                 
51# Setup computational domain
52#-------------------------------------------------------------------------------                                 
53domain = Domain(meshname, use_cache = True, verbose = True)
54
55print 'Number of triangles = ', len(domain)
56print 'The extent is ', domain.get_extent()
57print domain.statistics()
58
59domain.set_name('slide')
60domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
61domain.set_minimum_storable_height(0.01)
62
63#-------------------------------------------------------------------------------                                 
64# Setup initial conditions
65#-------------------------------------------------------------------------------
66from math import tan
67depth = 4.6
68b = 9.144 # from Synolakis paper
69def topography(x,y):
70    z = []
71    for xi in x:
72        if xi < b: 
73            z.append(-xi/2.)
74        else:
75            z.append(-depth)
76        #print 'hello', xi,b, z
77    return z
78domain.set_quantity('stage', 0.0)
79domain.set_quantity('friction', 0.0) 
80domain.set_quantity('elevation', topography,
81                    use_cache = True,
82                    verbose = True,
83                    alpha = 0.1
84                    )
85
86#-------------------------------------------------------------------------------
87# Set up scenario (tsunami_source is a callable object used with set_quantity)
88#-------------------------------------------------------------------------------
89from smf import slide_tsunami
90from math import atan, degrees
91
92a = .4572 # from Synolakis paper
93tsunami_source = slide_tsunami(length=0.91,
94                               width=0.61,
95                               depth=.1, # delta from Synolakis paper
96                               slope=degrees(atan(0.5)),
97                               thickness=a, 
98                               x0=0.75, 
99                               y0=width/2, 
100                               alpha=0.,
101                               #dx=0.001,
102                               domain=domain,
103                               verbose=True)
104
105#-------------------------------------------------------------------------------                                 
106# Setup boundary conditions
107#-------------------------------------------------------------------------------
108print 'Available boundary tags', domain.get_boundary_tags()
109
110Br = Reflective_boundary(domain)
111Bd = Dirichlet_boundary([0,0,0])
112
113domain.set_boundary( {'top': Bd,  'bottom': Bd, 'right': Bd, 'left': Bd} )
114
115
116#-------------------------------------------------------------------------------                                 
117# Evolve system through time
118#-------------------------------------------------------------------------------
119import time
120from Numeric import allclose
121from anuga.abstract_2d_finite_volumes.quantity import Quantity
122
123t0 = time.time()
124
125for t in domain.evolve(yieldstep = 1, finaltime = 500): 
126    domain.write_time()
127    stagestep = domain.get_quantity('stage') 
128
129    if allclose(t, 1):
130        slide = Quantity(domain)
131        slide.set_values(tsunami_source)
132        domain.set_quantity('stage', slide + stagestep)
133   
134print 'That took %.2f seconds' %(time.time()-t0)
135
136print 'finished'
Note: See TracBrowser for help on using the repository browser.