# source:anuga_work/development/landslide_4/benchmark_problem4.py@5572

Last change on this file since 5572 was 5572, checked in by bridgette, 16 years ago

update script for bridgettes use

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