source: anuga_core/source/anuga/damage_modelling/test_inundation_damage.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: 14.0 KB
Line 
1#!/usr/bin/env python
2#
3
4import unittest
5import tempfile
6import os
7import time
8import csv
9
10#from anuga.damage.inundation_damage import _calc_collapse_structures
11from inundation_damage import *
12from anuga.geospatial_data.geospatial_data import Geospatial_data
13from anuga.pmesh.mesh import Mesh
14from anuga.coordinate_transforms.geo_reference import Geo_reference
15from anuga.shallow_water import Domain, Transmissive_boundary
16from anuga.utilities.numerical_tools import mean
17from anuga.shallow_water.data_manager import get_dataobject
18
19from Numeric import zeros, Float, allclose
20
21
22def elevation_function(x, y):
23    return -x
24
25class Test_inundation_damage(unittest.TestCase):
26    def setUp(self):
27        #print "****set up****"
28        # Create an sww file
29       
30        # Set up an sww that has a geo ref.
31        # have it cover an area in Australia.  'gong maybe
32        #Don't have many triangles though!
33       
34        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
35        #Zone:   56   
36        #Easting:  222908.705  Northing: 6233785.284
37        #Latitude:   -34  0 ' 0.00000 ''  Longitude: 150  0 ' 0.00000 ''
38        #Grid Convergence:  -1  40 ' 43.13 ''  Point Scale: 1.00054660
39
40        #geo-ref
41        #Zone:   56   
42        #Easting:  220000  Northing: 6230000
43
44
45        #have  a big area covered.
46
47        mesh_file = tempfile.mktemp(".tsh")
48       
49        points_lat_long = [[-33,152],[-35,152],[-35,150],[-33,150]]
50       
51        spat = Geospatial_data(data_points=points_lat_long,
52                               points_are_lats_longs=True)
53        points_ab = spat.get_data_points( absolute = True)
54
55       
56        geo =  Geo_reference(56,400000,6000000)
57        spat.set_geo_reference(geo)
58
59        m = Mesh()
60        m.add_vertices(spat)
61        m.auto_segment()
62        m.generate_mesh(verbose=False)
63        m.export_mesh_file(mesh_file)
64
65       
66        #Create shallow water domain
67        domain = Domain(mesh_file)
68
69        os.remove(mesh_file)
70       
71        domain.default_order=2
72        domain.beta_h = 0
73
74
75        #Set some field values
76        #domain.set_quantity('stage', 1.0)
77        domain.set_quantity('elevation', -0.5)
78        domain.set_quantity('friction', 0.03)
79
80
81        ######################
82        # Boundary conditions
83        B = Transmissive_boundary(domain)
84        domain.set_boundary( {'exterior': B})
85
86
87        ######################
88        #Initial condition - with jumps
89
90        bed = domain.quantities['elevation'].vertex_values
91        stage = zeros(bed.shape, Float)
92
93        h = 0.3
94        for i in range(stage.shape[0]):
95            if i % 2 == 0:
96                stage[i,:] = bed[i,:] + h
97            else:
98                stage[i,:] = bed[i,:]
99
100        domain.set_quantity('stage', stage)
101        domain.set_quantity('xmomentum', stage*22.0)
102        domain.set_quantity('ymomentum', stage*55.0)
103
104        domain.distribute_to_vertices_and_edges()
105
106
107        self.domain = domain
108
109        C = domain.get_vertex_coordinates()
110        self.X = C[:,0:6:2].copy()
111        self.Y = C[:,1:6:2].copy()
112
113        self.F = bed
114
115       
116        #sww_file = tempfile.mktemp("")
117        self.domain.set_name('datatest' + str(time.time()))
118        self.domain.format = 'sww'
119        self.domain.smooth = True
120        self.domain.reduction = mean
121
122        sww = get_dataobject(self.domain)
123        sww.store_connectivity()
124        sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
125        self.domain.time = 2.
126        sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
127        self.sww = sww # so it can be deleted
128       
129        #Create a csv file
130        self.csv_file = tempfile.mktemp(".csv")
131        fd = open(self.csv_file,'wb')
132        writer = csv.writer(fd)
133        writer.writerow(['LONGITUDE','LATITUDE','STR_VALUE','C_VALUE','ROOF_TYPE','WALLS', 'SHORE_DIST'])
134        writer.writerow(['151.5','-34','199770','130000','Metal','Timber',20.])
135        writer.writerow(['151','-34.5','150000','76000','Metal','Double Brick',200.])
136        writer.writerow(['151','-34.25','150000','76000','Metal','Brick Veneer',200.])
137        fd.close()
138
139       
140    def tearDown(self):
141        #print "***** tearDown  ********"
142
143        # FIXME (Ole): Sometimes this fails - is the file open or is it sometimes not created?
144        try:
145            # Sometimes this fails - don't know why.
146            # Seems to be that the file is not created, since after it
147            # fails there are no sww files in the anuga directory
148            os.remove(self.sww.filename)
149        except OSError:
150            pass
151        os.remove(self.csv_file)
152
153   
154    def test_inundation_damage(self):
155
156        # Note, this isn't testing the results,
157        # just that is all runs
158        sww_file = self.domain.get_name() + "." + self.domain.format
159        #print "sww_file",sww_file
160        inundation_damage(sww_file, self.csv_file, verbose=False)
161
162   
163    def test_inundation_damage2(self):
164
165        # create mesh
166        mesh_file = tempfile.mktemp(".tsh")   
167        points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
168        m = Mesh()
169        m.add_vertices(points)
170        m.auto_segment()
171        m.generate_mesh(verbose=False)
172        m.export_mesh_file(mesh_file)
173       
174        #Create shallow water domain
175        domain = Domain(mesh_file)
176        os.remove(mesh_file)
177       
178        domain.default_order=2
179        domain.beta_h = 0
180
181        #Set some field values
182        domain.set_quantity('elevation', elevation_function)
183        domain.set_quantity('friction', 0.03)
184        domain.set_quantity('xmomentum', 22.0)
185        domain.set_quantity('ymomentum', 55.0)
186
187        ######################
188        # Boundary conditions
189        B = Transmissive_boundary(domain)
190        domain.set_boundary( {'exterior': B})
191
192        # This call mangles the stage values.
193        domain.distribute_to_vertices_and_edges()
194        domain.set_quantity('stage', 0.3)
195
196        #sww_file = tempfile.mktemp("")
197        domain.set_name('datatest' + str(time.time()))
198        domain.format = 'sww'
199        domain.smooth = True
200        domain.reduction = mean
201
202        sww = get_dataobject(domain)
203        sww.store_connectivity()
204        sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
205        domain.set_quantity('stage', -0.3)
206        domain.time = 2.
207        sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
208       
209        #Create a csv file
210        csv_file = tempfile.mktemp(".csv")
211        fd = open(csv_file,'wb')
212        writer = csv.writer(fd)
213        writer.writerow(['x','y','STR_VALUE','C_VALUE','ROOF_TYPE','WALLS', 'SHORE_DIST'])
214        writer.writerow([5.5,0.5,'10','130000','Metal','Timber',20])
215        writer.writerow([4.5,1.0,'150','76000','Metal','Double Brick',20])
216        writer.writerow([0.1,1.5,'100','76000','Metal','Brick Veneer',300])
217        writer.writerow([6.1,1.5,'100','76000','Metal','Brick Veneer',300])
218        fd.close()
219
220        sww_file = domain.get_name() + "." + domain.format
221        #print "sww_file",sww_file
222        inundation_damage(sww_file, csv_file, verbose=False)
223
224        csv_handle = Exposure_csv(csv_file)
225        struct_loss = csv_handle.get_column(EventDamageModel.STRUCT_LOSS_TITLE)
226        #print "struct_loss",struct_loss
227        struct_loss = [float(x) for x in struct_loss]
228        assert allclose(struct_loss,[10,150,16.9,0])
229        depth = csv_handle.get_column(EventDamageModel.MAX_DEPTH_TITLE)
230        #print "depth",depth
231        depth = [float(x) for x in depth]
232        assert allclose(depth,[5.5,4.5,0.1,-0.3])
233        os.remove(sww.filename)
234        os.remove(csv_file)
235       
236    def ztest_add_depth_and_momentum2csv(self):
237        sww_file = self.domain.get_name() + "." + self.domain.format
238        #print "sww_file",sww_file
239       
240        out_csv = tempfile.mktemp(".csv")
241        print "out_csv",out_csv
242        add_depth_and_momentum2csv(sww_file, self.csv_file,
243                                   out_csv, verbose=False)
244       
245    def test_calc_damage_percentages(self):
246        max_depths = [-0.3, 0.0, 1.0,-0.3, 0.0, 1.0,-0.3, 0.0, 1.0]
247        shore_distances = [100, 100, 100, 100, 100, 100, 100, 100, 100]
248        walls = ['Double Brick',
249                 'Double Brick',
250                 'Double Brick',
251                 'Timber',
252                 'Timber',
253                 'Timber',
254                 'Brick Veneer',
255                 'Brick Veneer',
256                 'Brick Veneer']
257        struct_costs = [10,
258                        10,
259                        10,
260                        10,
261                        10,
262                        10,
263                        1,
264                        1,
265                        1]
266        content_costs = [100,
267                        100,
268                        100,
269                        100,
270                        100,
271                        100,
272                        10,
273                        10,
274                        10]
275
276        edm = EventDamageModel(max_depths, shore_distances, walls,
277                               struct_costs, content_costs)
278        edm.calc_damage_percentages()
279        assert allclose(edm.struct_damage,[0.0,0.016,0.572,
280                                            0.0,0.016,0.618,
281                                            0.0,0.016,0.618])
282        assert allclose(edm.contents_damage,[0.0,0.013,0.970,
283                                             0.0,0.013,0.970,
284                                             0.0,0.013,0.970])
285        edm.calc_cost()
286        assert allclose(edm.struct_loss,[0.0,.16,5.72,
287                                            0.0,.16,6.18,
288                                            0.0,0.016,0.618])
289        assert allclose(edm.contents_loss,[0.0,1.3,97,
290                                             0.0,1.3,97,
291                                             0.0,0.13,9.7])
292       
293       
294    def test_calc_collapse_structures(self):
295        edm = EventDamageModel([0.0]*17, [0.0]*17, [0.0]*17,
296                               [0.0]*17, [0.0]*17)
297        edm.struct_damage = zeros(17,Float) 
298        edm.contents_damage = zeros(17,Float) 
299        collapse_probability = {0.4:[0], #0
300                                0.6:[1], #1
301                                0.5:[2], #1
302                                0.25:[3,4], #1
303                                0.1:[5,6,7,8], #0
304                                0.2:[9,10,11,12,13,14,15,16]} #2
305        edm._calc_collapse_structures(collapse_probability, verbose_csv=True)
306
307        self.failUnless( edm.struct_damage[0]  == 0.0 and
308                         edm.contents_damage[0]  == 0.0,
309                        'Error!')
310        self.failUnless( edm.struct_damage[1]  == 1.0 and
311                         edm.contents_damage[1]  == 1.0,
312                        'Error!')
313        self.failUnless( edm.struct_damage[2]  == 1.0 and
314                         edm.contents_damage[2]  == 1.0,
315                        'Error!')
316        self.failUnless( edm.struct_damage[3]+ edm.struct_damage[4] == 1.0 and
317                         edm.contents_damage[3] + edm.contents_damage[4] ==1.0,
318                        'Error!')
319        sum_struct = 0.0
320        sum_contents = 0.0
321        for i in [5,6,7,8]:
322            sum_struct += edm.struct_damage[i]
323            sum_contents += edm.contents_damage[i]
324        print "", 
325        self.failUnless( sum_struct == 0.0 and sum_contents  == 0.0,
326                        'Error!')
327        sum_struct = 0.0
328        sum_contents = 0.0
329        for i in [9,10,11,12,13,14,15,16]:
330            sum_struct += edm.struct_damage[i]
331            sum_contents += edm.contents_damage[i]
332        self.failUnless( sum_struct == 2.0 and sum_contents  == 2.0,
333                        'Error!')
334       
335    def test_calc_collapse_probability(self):
336        depth =          [0.0, 0.5, 0.5  , 1.5, 2.5, 4.5, 10000, 2.0]
337        shore_distance = [0.0, 125, 250.1, 0.0, 150, 225, 10000, 251]
338        dummy = depth
339        edm = EventDamageModel(depth, shore_distance, dummy, dummy, dummy)
340        struct_coll_prob = edm.calc_collapse_probability()
341        answer = {0.05:[1,7],
342                  0.6:[3],
343                  0.4:[4],
344                  0.5:[5],
345                  0.45:[6]}
346        #print "struct_coll_prob",struct_coll_prob
347        #print "answer",answer
348
349        self.failUnless( struct_coll_prob ==  answer,
350                        'Error!')
351       
352       
353    def test_calc_damage_and_costs(self):
354                             
355        max_depths = [-0.3, 0.0, 1.0,-0.3, 0.0, 1.0,-0.3, 0.0, 10.0]
356        shore_distances = [100, 100, 100, 100, 100, 100, 100, 100, 100]
357        walls = ['Double Brick',
358                 'Double Brick',
359                 'Double Brick',
360                 'Timber',
361                 'Timber',
362                 'Timber',
363                 'Brick Veneer',
364                 'Brick Veneer',
365                 'Brick Veneer']
366        struct_costs = [10,
367                        10,
368                        10,
369                        10,
370                        10,
371                        10,
372                        1,
373                        1,
374                        1]
375        content_costs = [100,
376                        100,
377                        100,
378                        100,
379                        100,
380                        100,
381                        10,
382                        10,
383                        10]
384
385        edm = EventDamageModel(max_depths, shore_distances, walls,
386                               struct_costs, content_costs)
387        results_dic = edm.calc_damage_and_costs(verbose_csv=True)
388        #print "results_dic",results_dic
389#-------------------------------------------------------------
390if __name__ == "__main__":
391    #suite = unittest.makeSuite(Test_inundation_damage,'test_in_damage2')
392    suite = unittest.makeSuite(Test_inundation_damage,'test')
393    runner = unittest.TextTestRunner()
394    runner.run(suite)
395
Note: See TracBrowser for help on using the repository browser.