source: anuga_work/development/parallel/run_sw_lost_mass.py @ 5162

Last change on this file since 5162 was 3460, checked in by jack, 19 years ago

Moved parallel to development from inundation

File size: 2.6 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 os import sep
34sys.path.append('..'+sep+'pyvolution')
35
36from Numeric import array, zeros, Float
37# pmesh
38
39#from shallow_water import Domain
40
41from shallow_water import Domain
42
43# mesh partition routines
44
45from pmesh2domain import pmesh_to_domain_instance
46
47# read in the processor information
48
49
50#-------
51# Domain
52rect = zeros( 4, Float) # Buffer for results
53
54class Set_Stage:
55    """Set an initial condition with constant water height, for x0<x<x1
56    """
57
58    def __init__(self, x0=0.25, x1=0.5, h=1.0):
59        self.x0 = x0
60        self.x1 = x1
61        self.= h
62
63    def __call__(self, x, y):
64        return self.h*((x>self.x0)&(x<self.x1))
65
66
67
68    # read in the test files
69
70filename = 'test-100.tsh'
71
72nx = 1
73ny = 1
74
75domain = pmesh_to_domain_instance(filename, Domain)
76
77domain.set_quantity('stage', Set_Stage(200.0,300.0,5.0))
78rect = array(domain.xy_extent, Float)
79
80try:
81    domain.initialise_visualiser(rect=rect)
82    #domain.visualiser.coloring['stage'] = True
83    domain.visualiser.scale_z['stage'] = 0.2
84    domain.visualiser.scale_z['elevation'] = 0.05
85except:
86    print 'No visualiser'
87
88
89domain.default_order = 1
90
91#Boundaries
92from shallow_water import Transmissive_boundary, Reflective_boundary
93
94T = Transmissive_boundary(domain)
95R = Reflective_boundary(domain)
96domain.set_boundary( {'outflow': R, 'inflow': R, 'inner':R, 'exterior': R, 'open':R} )
97
98
99#domain.set_quantity('stage', quantities['stage'])
100#domain.set_quantity('elevation', quantities['elevation'])
101
102#domain.store = True
103#domain.filename = 'merimbula-%d' %domain.processor
104
105#---------
106# Evolution
107t0 = time.time()
108
109print 'No of elements %d'%(domain.number_of_elements)
110
111
112for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
113    #domain.write_time()
114    print '    Integral stage = ', domain.quantities['stage'].get_integral(),' Time = ',domain.time
115
116
Note: See TracBrowser for help on using the repository browser.