Changeset 3722
- Timestamp:
- Oct 10, 2006, 11:47:18 AM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_validation/automated_validation_tests/okushiri_tank_validation/validate_okushiri.py
r3718 r3722 1 """ Validation of the AnuGA implementation of the shallow water wave equation.1 """Automated validation of ANUGA. 2 2 3 This script sets up LWRU2 benchmark with initial condition stated 3 This script runs the benchmark on a regular grid and verifies that 4 the predicted timeseries at selected gauges are close enought to the observed 5 timeseries provided from the experiment. 4 6 5 See also 7 See anuga_validation/okushiri_2005 for scripts that produce simulations at 8 greater resolution. 6 9 7 http://www.cee.cornell.edu/longwave/index.cfm?page=benchmark&problem=28 9 Depth at western boundary is d = 13.5 cm10 10 """ 11 11 … … 13 13 from Numeric import array, zeros, Float, allclose 14 14 15 from anuga.utilities.numerical_tools import ensure_numeric 15 16 from anuga.shallow_water import Domain 16 17 from anuga.shallow_water import Reflective_boundary 17 18 from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary 19 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 18 20 from anuga.abstract_2d_finite_volumes.util import file_function 19 21 20 22 import project 21 23 24 finaltime = 22.5 25 timestep = 0.05 22 26 23 27 #------------------------- 24 28 # Create Domain 25 29 #------------------------- 30 #N = 50 31 N = 10 32 points, vertices, boundary = rectangular_cross(N, N/5*3, 33 len1=5.448, len2=3.402) 34 domain = Domain(points, vertices, boundary) 26 35 27 use_variable_mesh = True #Use large variable mesh generated by create_mesh.py 28 #use_variable_mesh = False #Use small structured mesh for speed 29 30 if use_variable_mesh is True: 31 print 'Creating domain from', project.mesh_filename 32 33 domain = Domain(project.mesh_filename, use_cache=True, verbose=True) 34 else: 35 print 'Creating regular from regular mesh' 36 N = 150 37 points, vertices, boundary = rectangular_cross(N, N/5*3, 38 len1=5.448, len2=3.402) 39 domain = Domain(points, vertices, boundary) 40 41 domain.set_name('lwru2') 36 domain.set_name('okushiri_automated_validation') 42 37 domain.set_default_order(2) 43 38 domain.set_minimum_storable_height(0.001) 44 39 domain.check_integrity() 45 print domain.statistics()46 40 47 41 … … 65 59 domain, 66 60 verbose=True) 67 68 61 Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function) 69 70 if use_variable_mesh is True: 71 domain.set_boundary({'wave': Bts, 'wall': Br}) 72 else: 73 domain.set_boundary({'left': Bts, 'right': Br, 'bottom': Br, 'top': Br}) 62 domain.set_boundary({'left': Bts, 'right': Br, 'bottom': Br, 'top': Br}) 74 63 75 64 76 65 #------------------------- 77 # Evolve through time 66 # Set up validation data 67 #------------------------- 68 gauge_locations = [[0.000, 1.696]] #Boundary 69 gauge_locations += [[4.521, 1.196], [4.521, 1.696], [4.521, 2.196]] #Ch 5-7-9 70 71 #Boundary input wave 72 filename = project.boundary_filename[:-4] + '.tms' 73 print 'Reading', filename 74 from Scientific.IO.NetCDF import NetCDFFile 75 from Numeric import allclose 76 fid = NetCDFFile(filename, 'r')# 77 78 input_time = fid.variables['time'][:] 79 input_stage = fid.variables['stage'][:] 80 81 82 # reference gauge timeseries 83 reference_time = [] 84 ch5 = [] 85 ch7 = [] 86 ch9 = [] 87 fid = open('output_ch5-7-9.txt') 88 lines = fid.readlines() 89 fid.close() 90 #for i, line in enumerate(lines[1:]): 91 for line in lines[1:]: 92 fields = line.split() 93 94 reference_time.append(float(fields[0])) 95 ch5.append(float(fields[1])/100) #cm2m 96 ch7.append(float(fields[2])/100) #cm2m 97 ch9.append(float(fields[3])/100) #cm2m 98 99 if fields[0] == str(finaltime): break 100 101 102 assert reference_time[0] == 0.0 103 assert reference_time[-1] == finaltime 104 105 reference_gauge_values = [ensure_numeric(input_stage), 106 ensure_numeric(ch5), 107 ensure_numeric(ch7), 108 ensure_numeric(ch9)] 109 110 111 112 predicted_gauge_values = [[], [],[],[]] 113 114 115 #------------------------- 116 # Evolve through time 78 117 #------------------------- 79 118 import time 80 119 t0 = time.time() 81 120 82 w_max = 0 83 for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5):121 for j, t in enumerate(domain.evolve(yieldstep = timestep, 122 finaltime = finaltime)): 84 123 domain.write_time() 85 124 86 w = domain.get_maximum_inundation_elevation() 87 x, y = domain.get_maximum_inundation_location() 88 print ' Coastline elevation = %.2f at (x,y)=(%.2f, %.2f)' %(w, x, y) 89 print 125 # Record gauge values 126 stage = domain.get_quantity('stage') 127 for i, (x, y) in enumerate(gauge_locations): 128 w = stage.get_values(interpolation_points=[[x,y]]) 129 predicted_gauge_values[i].append(w) 130 131 print 'difference', x, y, w - reference_gauge_values[i][j] 90 132 91 if w > w_max:92 w_max = w93 x_max = x94 y_max = y95 133 134 #------------------------- 135 # Validate 136 #------------------------- 137 for i, (x, y) in enumerate(gauge_locations): 138 predicted_gauge_values[i] = ensure_numeric(predicted_gauge_values[i]) 96 139 97 print '**********************************************'98 print 'Max coastline elevation = %.2f at (%.2f, %.2f)' %(w_max, x_max, y_max)99 print '**********************************************'100 140 101 141 sqsum = sum((predicted_gauge_values[i] - reference_gauge_values)**2) 142 denom = sum(reference_gauge_values**2) 143 144 print i, sqsum, sqsum/denom 145 146 msg 'Rms error is %f' %(sqsum/denom) 147 assert sqsum/denom < 0.01, msg 148 149 try: 150 from pylab import * 151 except: 152 pass 153 else: 154 ion() 155 hold(False) 156 plot(reference_time, reference_gauge_values[i], 'r.-', 157 reference_time, predicted_gauge_values[i], 'k-') 158 title('Gauge %s' %name) 159 xlabel('time(s)') 160 ylabel('stage (m)') 161 legend(('Observed', 'Modelled'), shadow=True, loc='upper left') 162 savefig(name, dpi = 300) 163 164 raw_input('Next') 165 166 167 try: 168 from pylab import * 169 except: 170 pass 171 else: 172 show() 173 102 174 103 175 print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracChangeset
for help on using the changeset viewer.