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

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

Bridgette starting work on Benchmark problem 4

File size: 4.6 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 = 10.
33width = 10.
34depth = 10.
35bounding_polygon = [[0,width],[0,0],[length,0],[length,width]]
36meshname = 'slide.msh'
37from caching import cache
38_ = cache(create_mesh_from_regions,
39          bounding_polygon,
40           {'boundary_tags': {'top': [0], 'left': [1], 
41                              'bottom': [2], 'right': [3]},
42           'maximum_triangle_area': remainder_res,
43           'filename': meshname},
44           #'interior_regions': interior_regions},
45          verbose = True, evaluate=False)
46print 'created mesh'
47
48#-------------------------------------------------------------------------------                                 
49# Setup computational domain
50#-------------------------------------------------------------------------------                                 
51domain = Domain(meshname, use_cache = True, verbose = True)
52
53print 'Number of triangles = ', len(domain)
54print 'The extent is ', domain.get_extent()
55print domain.statistics()
56
57domain.set_name('slide')
58domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
59domain.set_minimum_storable_height(0.01)
60
61#-------------------------------------------------------------------------------                                 
62# Setup initial conditions
63#-------------------------------------------------------------------------------
64from math import tan
65def topography(x,y):
66    z = -depth + tan(0.5)*x
67    return z
68domain.set_quantity('stage', 0.0)
69domain.set_quantity('friction', 0.0) 
70domain.set_quantity('elevation', topography,
71                    use_cache = True,
72                    verbose = True,
73                    alpha = 0.1
74                    )
75
76#-------------------------------------------------------------------------------
77# Set up scenario (tsunami_source is a callable object used with set_quantity)
78#-------------------------------------------------------------------------------
79from smf import slide_tsunami
80
81tsunami_source = slide_tsunami(length=2.,
82                               width=1.,
83                               depth=1.,
84                               slope=45,
85                               thickness=0.1, 
86                               x0=length/2, 
87                               y0=width/2, 
88                               alpha=0.,
89                                         #dx=0.001,
90                               domain=domain)
91                               #verbose=True)
92
93#-------------------------------------------------------------------------------                                 
94# Setup boundary conditions
95#-------------------------------------------------------------------------------
96print 'Available boundary tags', domain.get_boundary_tags()
97
98Br = Reflective_boundary(domain)
99Bd = Dirichlet_boundary([tide,0,0])
100
101domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd} )
102
103
104#-------------------------------------------------------------------------------                                 
105# Evolve system through time
106#-------------------------------------------------------------------------------
107import time
108from Numeric import allclose
109from anuga.abstract_2d_finite_volumes.quantity import Quantity
110
111t0 = time.time()
112
113for t in domain.evolve(yieldstep = 1, finaltime = 10): 
114    domain.write_time()
115    stagestep = domain.get_quantity('stage') 
116
117    if allclose(t, 1):
118        slide = Quantity(domain)
119        slide.set_values(tsunami_source)
120        domain.set_quantity('stage', slide + stagestep)
121   
122print 'That took %.2f seconds' %(time.time()-t0)
123
124print 'finished'
Note: See TracBrowser for help on using the repository browser.