source: anuga_core/source/anuga_parallel/run_sw_lost_mass.py @ 3846

Last change on this file since 3846 was 3846, checked in by ole, 17 years ago

Refactored references to domain.filename away.
Use

domain.set_name()
domain.get_name()

instead.

File size: 2.4 KB
Line 
1#!/usr/bin/env/python
2#########################################################
3#
4# Shows problem of conservation of mass. Due to material being
5# removed if too shallow. Can lead to problems if timestep
6# set very small.
7#
8#
9#
10#  Authors: Linda Stals, Steve Roberts and Matthew Hardy,
11# June 2005
12#
13#
14#
15#########################################################
16
17
18##############################################
19# Set a small yieldstep (say 0.005) to see
20# lose of conservation
21# Set it to a larger value (say 0.1) and have
22# conservation "restored"
23##############################################
24yieldstep = 0.005
25#yieldstep = 0.08
26#yieldstep = 0.05
27finaltime = 100.0
28
29import sys
30import time
31
32
33from Numeric import array, zeros, Float
34
35#from shallow_water import Domain
36
37from anuga.shallow_water import Domain
38
39# mesh partition routines
40
41from anuga.abstract_2d_finite_volumes.pmesh2domain\
42     import pmesh_to_domain_instance
43
44# read in the processor information
45
46
47#-------
48# Domain
49rect = zeros( 4, Float) # Buffer for results
50
51class Set_Stage:
52    """Set an initial condition with constant water height, for x0<x<x1
53    """
54
55    def __init__(self, x0=0.25, x1=0.5, h=1.0):
56        self.x0 = x0
57        self.x1 = x1
58        self.= h
59
60    def __call__(self, x, y):
61        return self.h*((x>self.x0)&(x<self.x1))
62
63
64
65    # read in the test files
66
67filename = 'test-100.tsh'
68
69nx = 1
70ny = 1
71
72domain = pmesh_to_domain_instance(filename, Domain)
73
74domain.set_quantity('stage', Set_Stage(200.0,300.0,5.0))
75rect = array(domain.xy_extent, Float)
76
77try:
78    domain.initialise_visualiser(rect=rect)
79    #domain.visualiser.coloring['stage'] = True
80    domain.visualiser.scale_z['stage'] = 0.2
81    domain.visualiser.scale_z['elevation'] = 0.05
82except:
83    print 'No visualiser'
84
85
86domain.default_order = 1
87
88#Boundaries
89from shallow_water import Transmissive_boundary, Reflective_boundary
90
91T = Transmissive_boundary(domain)
92R = Reflective_boundary(domain)
93domain.set_boundary( {'outflow': R, 'inflow': R, 'inner':R, 'exterior': R, 'open':R} )
94
95
96
97#---------
98# Evolution
99t0 = time.time()
100
101print 'No of elements %d'%(domain.number_of_elements)
102
103
104for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
105    #domain.write_time()
106    print '    Integral stage = ', domain.quantities['stage'].get_integral(),' Time = ',domain.time
107
108
Note: See TracBrowser for help on using the repository browser.