Changeset 3722

Timestamp:

Oct 10, 2006, 11:47:18 AM (16 years ago)

Author:

ole

Message:

Work on automated okushiri island validation

File:

• anuga_validation/automated_validation_tests/okushiri_tank_validation/validate_okushiri.py (modified) (3 diffs)

anuga_validation/automated_validation_tests/okushiri_tank_validation/validate_okushiri.py

@@ -1 @@
-"""Validation of the AnuGA implementation of the shallow water wave equation.
+"""Automated validation of ANUGA.
 
-This script sets up LWRU2 benchmark with initial condition stated
+This script runs the benchmark on a regular grid and verifies that
+the predicted timeseries at selected gauges are close enought to the observed
+timeseries provided from the experiment.
 
-See also
+See anuga_validation/okushiri_2005 for scripts that produce simulations at
+greater resolution.
 
-http://www.cee.cornell.edu/longwave/index.cfm?page=benchmark&problem=2
-
-Depth at western boundary is d = 13.5 cm
 """
 
@@ -13 +13 @@
 from Numeric import array, zeros, Float, allclose
 
+from anuga.utilities.numerical_tools import ensure_numeric
 from anuga.shallow_water import Domain
 from anuga.shallow_water import Reflective_boundary
 from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
 from anuga.abstract_2d_finite_volumes.util import file_function
 
 import project
 
+finaltime = 22.5
+timestep = 0.05
 
 #-------------------------
 # Create Domain
 #-------------------------
+#N = 50
+N = 10
+points, vertices, boundary = rectangular_cross(N, N/5*3,
+                                               len1=5.448, len2=3.402)
+domain = Domain(points, vertices, boundary)
 
-use_variable_mesh = True #Use large variable mesh generated by create_mesh.py
-#use_variable_mesh = False #Use small structured mesh for speed
-
-if use_variable_mesh is True:
-    print 'Creating domain from', project.mesh_filename
-
-    domain = Domain(project.mesh_filename, use_cache=True, verbose=True)
-else:
-    print 'Creating regular from regular mesh'
-    N = 150
-    points, vertices, boundary = rectangular_cross(N, N/5*3,
-                                                   len1=5.448, len2=3.402)
-    domain = Domain(points, vertices, boundary)
-
-domain.set_name('lwru2')
+domain.set_name('okushiri_automated_validation')
 domain.set_default_order(2)
 domain.set_minimum_storable_height(0.001)
 domain.check_integrity()
-print domain.statistics()
 
 
@@ -65 +59 @@
                          domain,
                          verbose=True)
-
 Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
-
-if use_variable_mesh is True:
-    domain.set_boundary({'wave': Bts, 'wall': Br})
-else:
-    domain.set_boundary({'left': Bts, 'right': Br, 'bottom': Br, 'top': Br})
+domain.set_boundary({'left': Bts, 'right': Br, 'bottom': Br, 'top': Br})
 
 
 #-------------------------
-# Evolve through time
+# Set up validation data
+#-------------------------
+gauge_locations = [[0.000, 1.696]]  #Boundary
+gauge_locations += [[4.521, 1.196], [4.521, 1.696], [4.521, 2.196]] #Ch 5-7-9
+
+#Boundary input wave
+filename = project.boundary_filename[:-4] + '.tms'
+print 'Reading', filename
+from Scientific.IO.NetCDF import NetCDFFile
+from Numeric import allclose
+fid = NetCDFFile(filename, 'r')#
+
+input_time = fid.variables['time'][:]
+input_stage = fid.variables['stage'][:]
+
+
+# reference gauge timeseries
+reference_time = []
+ch5 = []
+ch7 = []
+ch9 = []
+fid = open('output_ch5-7-9.txt')
+lines = fid.readlines()
+fid.close()
+#for i, line in enumerate(lines[1:]):
+for line in lines[1:]:
+    fields = line.split()
+
+    reference_time.append(float(fields[0]))
+    ch5.append(float(fields[1])/100)   #cm2m
+    ch7.append(float(fields[2])/100)   #cm2m
+    ch9.append(float(fields[3])/100)   #cm2m
+
+    if fields[0] == str(finaltime): break    
+    
+
+assert reference_time[0] == 0.0
+assert reference_time[-1] == finaltime
+
+reference_gauge_values = [ensure_numeric(input_stage),
+                          ensure_numeric(ch5),
+                          ensure_numeric(ch7),
+                          ensure_numeric(ch9)] 
+
+
+
+predicted_gauge_values = [[], [],[],[]]
+
+
+#-------------------------
+# Evolve through time 
 #-------------------------
 import time
 t0 = time.time()
 
-w_max = 0
-for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5):
+for j, t in enumerate(domain.evolve(yieldstep = timestep,
+                                    finaltime = finaltime)):
     domain.write_time()
 
-    w = domain.get_maximum_inundation_elevation()
-    x, y = domain.get_maximum_inundation_location()
-    print '  Coastline elevation = %.2f at (x,y)=(%.2f, %.2f)' %(w, x, y)
-    print
+    # Record gauge values
+    stage = domain.get_quantity('stage')
+    for i, (x, y) in enumerate(gauge_locations):
+        w = stage.get_values(interpolation_points=[[x,y]])
+        predicted_gauge_values[i].append(w)
+
+        print 'difference', x, y, w - reference_gauge_values[i][j]
     
-    if w > w_max:
-        w_max = w
-        x_max = x
-        y_max = y
 
+#-------------------------
+# Validate
+#-------------------------
+for i, (x, y) in enumerate(gauge_locations):        
+    predicted_gauge_values[i] = ensure_numeric(predicted_gauge_values[i])
 
-print '**********************************************'
-print 'Max coastline elevation = %.2f at (%.2f, %.2f)' %(w_max, x_max, y_max)
-print '**********************************************'
     
-   
+    sqsum = sum((predicted_gauge_values[i] - reference_gauge_values)**2)
+    denom = sum(reference_gauge_values**2)
+
+    print i, sqsum, sqsum/denom
+
+    msg 'Rms error is %f' %(sqsum/denom)
+    assert sqsum/denom < 0.01, msg
+
+    try:
+        from pylab import *
+    except:
+        pass
+    else:
+        ion()
+        hold(False)
+        plot(reference_time, reference_gauge_values[i], 'r.-',
+             reference_time, predicted_gauge_values[i], 'k-')
+        title('Gauge %s' %name)
+        xlabel('time(s)')
+        ylabel('stage (m)')    
+        legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
+        savefig(name, dpi = 300)
+        
+        raw_input('Next')
+        
+
+try:
+    from pylab import *
+except:
+    pass
+else:
+    show()
+    
 
 print 'That took %.2f seconds' %(time.time()-t0)