# source:anuga_core/source/obsolete_code/_test_sww2domain2.py@7585

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

Moved obsolete code and examples to their appropriate locations

File size: 2.8 KB
Line
1"""Simple example of shallow water wave equation using Pyvolution
2
3Water driven by linear slope and Dirichlet boundary
4
5"""
6
7######################
8# Module imports
9#
10from mesh_factory import rectangular
11from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
12     Constant_height, Time_boundary, Transmissive_boundary
13from Numeric import array
14
15#Create basic mesh
16points, vertices, boundary = rectangular(10, 10)
17
18#Create shallow water domain
19domain = Domain(points, vertices, boundary)
20domain.smooth = False
21domain.visualise = False
22domain.store = True
23domain.filename = 'bedslope'
24domain.default_order=2
25domain.quantities_to_be_stored=['stage']
26
27#######################
28#Bed-slope and friction
29domain.set_quantity('elevation', lambda x,y: -x/3)
30domain.set_quantity('friction', 0.1)
31
32######################
33# Boundary conditions
34from math import sin, pi
35Br = Reflective_boundary(domain)
36Bt = Transmissive_boundary(domain)
37Bd = Dirichlet_boundary([0.2,0.,0.])
38Bw = Time_boundary(domain=domain,
39                   f=lambda t: [(0.1*sin(t*2*pi)), 0.0, 0.0])
40
41domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
42
43######################
44#Initial condition
45h = 0.05
46elevation = domain.quantities['elevation'].vertex_values
47domain.set_quantity('stage', elevation + h)
48#elevation = domain.get_quantity('elevation',location='unique vertices')
49#domain.set_quantity('stage', elevation + h,location='unique vertices')
50
51domain.check_integrity()
52######################
53#Evolution
54for t in domain.evolve(yieldstep = 1, finaltime = 2.0):
55    pass
56    #domain.write_time()
57
58
59##NOW TEST IT!!!
60
61from anuga.pyvolution.data_manager import sww2domain
62from Numeric import allclose
63
65
66try:
67    domain2 = sww2domain(filename,verbose=False)
68    assert True == False
69except:
70    filler = 0
71    domain2 = sww2domain(filename,fail_if_NaN=False,NaN_filler = filler,verbose=False)
72
73bits = ['xllcorner','yllcorner','vertex_coordinates','time','starttime']
74
75for quantity in ['elevation']+domain.quantities_to_be_stored:
76    bits.append('get_quantity("%s")'%quantity)
77
78for bit in bits:
79#    print 'testing that domain.'+bit+' has been restored'
80    assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
81
82#print max(max(domain2.get_quantity('xmomentum')))
83#print min(min(domain2.get_quantity('xmomentum')))
84#print max(max(domain2.get_quantity('ymomentum')))
85#print min(min(domain2.get_quantity('ymomentum')))
86
87assert max(max(domain2.get_quantity('xmomentum')))==filler
88assert min(min(domain2.get_quantity('xmomentum')))==filler
89assert max(max(domain2.get_quantity('ymomentum')))==filler
90assert min(min(domain2.get_quantity('ymomentum')))==filler
91
92#print 'passed'
93
94#cleanup
95#import os