# source:trunk/anuga_core/source/anuga_parallel/run_sw_lost_mass_simpler.py@8022

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

Removed all references to pyvolution in parallel code

File size: 1.9 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, Matthew Hardy, Ole Nielsen
11#  July 2005
12#
13#
14#########################################################
15
16
17##############################################
18# Change min_depth and see how larger values aggravate the problem
19yieldstep = 0.1
20finaltime = 10.0
21#min_depth = 1.0e-4
22min_depth = 1.0e-2
23
24
25from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
26from anuga.shallow_water import Domain, Reflective_boundary
27
28
29
30#Create shallow water domain
31points, vertices, boundary = rectangular(50, 50, len1=500, len2=500)
32domain = Domain(points, vertices, boundary)
33domain.smooth = False
34domain.visualise = True
35domain.default_order = 1
36domain.minimum_allowed_height = min_depth
37print 'Extent', domain.get_extent()
38
39# Set initial condition
40class Set_IC:
41    """Set an initial condition with a constant value, for x0<x<x1
42    """
43
44    def __init__(self, x0=0.25, x1=0.5, h=1.0):
45        self.x0 = x0
46        self.x1 = x1
47        self.= h
48
49    def __call__(self, x, y):
50        return self.h*((x>self.x0)&(x<self.x1))
51
52
53domain.set_quantity('stage', Set_IC(200.0,300.0,5.0))
54
55try:
56    domain.initialise_visualiser()
57except:
58    print 'No visualiser'
59else:
60    domain.visualiser.scale_z['stage'] = 0.2
61    domain.visualiser.scale_z['elevation'] = 0.05
62
63
64#Boundaries
65R = Reflective_boundary(domain)
66domain.set_boundary( {'left': R, 'right': R, 'top':R, 'bottom': R} )
67
68
69# Evolution
70print 'Minimal allowed water height = ', domain.minimum_allowed_height
71for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
72    print 'Integral stage = ', domain.quantities['stage'].get_integral(),' Time = ',domain.time
73
74
Note: See TracBrowser for help on using the repository browser.