source: anuga_core/source/anuga/examples/bedslope_bsod_test.py @ 4150

Last change on this file since 4150 was 3514, checked in by duncan, 19 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

File size: 3.3 KB
Line 
1"""
2
3Test example to demonstrate the "Black Screen of Death"
4
5* Simple mesh set up on square with square internal polygon
6* Set elevation using geospatial data
7* Duplicate point in internal polygon (or bounding polygon) can
8  cause the "Black Screen of Death" only when the elevation is set
9  using the geospatial data.
10
11Jane Sexton 2006
12
13"""
14
15######################
16# Module imports
17#
18from anuga.pmesh.mesh_interface import create_mesh_from_regions
19from anuga.pyvolution.pmesh2domain import pmesh_to_domain_instance
20from caching import cache
21from anuga.pyvolution.shallow_water import Domain, Reflective_boundary,\
22     Dirichlet_boundary, Time_boundary, Transmissive_boundary
23from anuga.geospatial_data.geospatial_data import *
24
25
26# bounding polygon
27c0 = [331000, 6255000]
28c1 = [345000, 6255000]
29c2 = [345000, 6267000]
30c3 = [331000, 6267000]
31bound_poly = [c0, c1, c2, c3]
32
33# set up internal region
34p0 = [335000, 6260000]
35p1 = [340000, 6260000]
36p2 = [340000, 6265000]
37p3 = [335000, 6265000]
38interior_polygon = [p0, p1, p2, p3]
39# duplicate point, this should fail, i.e. conjugate gradient doesn't converge
40interior_polygon = [p0, p1, p2, p3, p0] 
41
42interior_region_area = 2500000000 #25000
43interior_regions = [[interior_polygon, interior_region_area]]
44
45meshname = 'bedslope_bsod.tsh'
46
47exterior_region_area = 10000000000000000 #100000
48# create a mesh
49_ = cache(create_mesh_from_regions,
50          bound_poly,
51          {'boundary_tags': {'bottom': [0], 'right': [1],
52                             'top': [2], 'left': [3]},
53           'maximum_triangle_area': exterior_region_area,
54           'filename': meshname,           
55           'interior_regions': interior_regions},
56          evaluate = True,
57          verbose = True)
58
59#Create shallow water domain
60domain = pmesh_to_domain_instance(meshname, Domain,
61                                  use_cache = True,
62                                  verbose = True)
63
64domain.set_name('bedslope_bsod')
65domain.set_datadir('.')                      #Use current directory for output
66domain.set_quantities_to_be_stored('stage')  #See shallow_water.py
67
68# set up geospatial data object
69test_pts = []
70test_elev = []
71a = [340000, 6264000]
72b = [340000, 6266000]
73c = [335000, 6259000]
74d = [335000, 6266000]
75e = [335000, 6258000]
76f = [335000, 6259000]
77g = [341000, 6260000]
78h = [341000, 6256000]
79test_pts = [a, b, c, d, e, f, g, h]
80test_elev = [1.0, 4.0, 3.0, 0.1, 5, -100.0, -200, -15]
81G = Geospatial_data(test_pts, test_elev)
82
83# set geospatial data object to elevation
84domain.set_quantity('elevation', G)
85# if set elevation as a constant and duplicate a point in the interior region
86# the black screen of death will NOT be seen
87#domain.set_quantity('elevation', 1)
88domain.set_quantity('friction', 0.0)
89domain.set_quantity('stage', 0.0)
90
91
92######################
93# Boundary conditions
94from math import pi, sin
95Br = Reflective_boundary(domain)
96Bt = Transmissive_boundary(domain)
97Bd = Dirichlet_boundary([0.2,0.,0.])
98Bw = Time_boundary(domain=domain,
99                   f=lambda t: [(-1000*sin(t*2*pi)), 0.0, 0.0])
100
101print 'Tags are ', domain.get_boundary_tags()
102
103domain.set_boundary({'left': Bd, 'right': Bw, 'top': Br, 'bottom': Br})
104
105######################
106#Evolution
107domain.check_integrity()
108
109for t in domain.evolve(yieldstep = 5, finaltime = 50.0):
110    domain.write_time()
Note: See TracBrowser for help on using the repository browser.