Changeset 2679

Timestamp:

Apr 7, 2006, 3:49:48 PM (18 years ago)

Author:

duncan

Message:

checking in so I can run some tests

Location:

inundation/pyvolution

Files:

• test_util.py (modified) (7 diffs)
• util.py (modified) (1 diff)

inundation/pyvolution/test_util.py

@@ -5 +5 @@
 from Numeric import zeros, array, allclose, Float
 from math import sqrt, pi
+import tempfile
 
 from util import *
@@ -10 +11 @@
 from data_manager import timefile2netcdf
 
+from utilities.numerical_tools import INF
 
 def test_function(x, y):
@@ -498 +500 @@
 
 
-    def test_spatio_temporal_file_function_time(self):
+    def qtest_spatio_temporal_file_function_time(self):
         """Test that File function interpolates correctly
         between given times.
@@ -531 +533 @@
         points, vertices, boundary =\
                 rectangular(4, 4, 15, 30, origin = (0, -20))
-
+        print "points", points
 
         #print 'Number of elements', len(vertices)
@@ -571 +573 @@
 
         interpolation_points = [[0,-20], [1,0], [0,1], [1.1, 3.14], [10,-12.5]]
-
-
-
+      
         #Deliberately set domain.starttime to too early
         domain.starttime = start - 1
@@ -610 +610 @@
                     q1 = F(t+60, point_id=id)
 
-
+                if q0 == INF:
+                    actual = q0
+                else:
+                    actual = (k*q1 + (6-k)*q0)/6
                 q = F(t, point_id=id)
                 #print i, k, t, q
                 #print ' ', q0
                 #print ' ', q1
-                #print 's', (k*q1 + (6-k)*q0)/6
+                print "q",q
+                print "actual", actual
                 #print
-                assert allclose(q, (k*q1 + (6-k)*q0)/6)
+                if q0 == INF:
+                     self.failUnless( q == actual, 'Fail!')
+                else:
+                    assert allclose(q, actual)
 
 
@@ -667 +674 @@
 
 
+
+    def xtest_spatio_temporal_file_function_time(self):
+        # Test that File function interpolates correctly
+        # When some points are outside the mesh
+
+        import os, time
+        from config import time_format
+        from Numeric import sin, pi, exp
+        from mesh_factory import rectangular
+        from shallow_water import Domain
+        import data_manager 
+        from pmesh.mesh_interface import create_mesh_from_regions
+        finaltime = 1200
+        
+        filename = tempfile.mktemp()
+        #print "filename",filename 
+        filename = 'test_file_function'
+
+        meshfilename = tempfile.mktemp(".tsh")
+
+        boundary_tags = {'walls':[0,1],'bom':[2]}
+        
+        polygon_absolute = [[0,-20],[10,-20],[10,15],[-20,15]]
+        
+        create_mesh_from_regions(polygon_absolute,
+                                 boundary_tags,
+                                 10000000,
+                                 filename=meshfilename)
+        domain = Domain(mesh_filename=meshfilename)
+        domain.smooth = False
+        domain.default_order = 2
+        domain.set_datadir('.')
+        domain.set_name(filename)
+        domain.store = True
+        domain.format = 'sww'   #Native netcdf visualisation format
+
+        #print points
+        start = time.mktime(time.strptime('2000', '%Y'))
+        domain.starttime = start
+        
+
+        #Store structure
+        domain.initialise_storage()
+
+        #Compute artificial time steps and store
+        dt = 60  #One minute intervals
+        t = 0.0
+        while t <= finaltime:
+            #Compute quantities
+            f1 = lambda x,y: 3*x - y**2 + 2*t + 4
+            domain.set_quantity('stage', f1)
+
+            f2 = lambda x,y: x+y+t**2
+            domain.set_quantity('xmomentum', f2)
+
+            f3 = lambda x,y: x**2 + y**2 * sin(t*pi/600)
+            domain.set_quantity('ymomentum', f3)
+
+            #Store and advance time
+            domain.time = t
+            domain.store_timestep(domain.conserved_quantities)
+            t += dt
+
+        interpolation_points = [[1,0]]
+        interpolation_points = [[100,1000]]
+        
+        interpolation_points = [[0,-20], [1,0], [0,1], [1.1, 3.14], [10,-12.5],
+                                [78787,78787],[7878,3432]]
+            
+        #Deliberately set domain.starttime to too early
+        domain.starttime = start - 1
+
+        #Create file function
+        F = file_function(filename + '.sww', domain,
+                          quantities = domain.conserved_quantities,
+                          interpolation_points = interpolation_points)
+
+        #Check that FF updates fixes domain starttime
+        assert allclose(domain.starttime, start)
+
+        #Check that domain.starttime isn't updated if later
+        domain.starttime = start + 1
+        F = file_function(filename + '.sww', domain,
+                          quantities = domain.conserved_quantities,
+                          interpolation_points = interpolation_points)
+        assert allclose(domain.starttime, start+1)
+        domain.starttime = start
+
+
+        #Check linear interpolation in time
+        # checking points inside and outside the mesh
+        F = file_function(filename + '.sww', domain,
+                          quantities = domain.conserved_quantities,
+                          interpolation_points = interpolation_points)
+        
+        for id in range(len(interpolation_points)):
+            x = interpolation_points[id][0]
+            y = interpolation_points[id][1]
+
+            for i in range(20):
+                t = i*10
+                k = i%6
+
+                if k == 0:
+                    q0 = F(t, point_id=id)
+                    q1 = F(t+60, point_id=id)
+
+                if q0 == INF:
+                    actual = q0
+                else:
+                    actual = (k*q1 + (6-k)*q0)/6
+                q = F(t, point_id=id)
+                #print i, k, t, q
+                #print ' ', q0
+                #print ' ', q1
+                #print "q",q
+                #print "actual", actual
+                #print
+                if q0 == INF:
+                     self.failUnless( q == actual, 'Fail!')
+                else:
+                    assert allclose(q, actual)
+
+        # now lets check points inside the mesh
+        interpolation_points = [[0,-20], [1,0], [0,1], [1.1, 3.14]] #, [10,-12.5]] - this point doesn't work WHY?
+        interpolation_points = [[10,-12.5]]
+            
+        print "len(interpolation_points)",len(interpolation_points) 
+        F = file_function(filename + '.sww', domain,
+                          quantities = domain.conserved_quantities,
+                          interpolation_points = interpolation_points)
+
+        domain.starttime = start
+
+
+        #Check linear interpolation in time
+        F = file_function(filename + '.sww', domain,
+                          quantities = domain.conserved_quantities,
+                          interpolation_points = interpolation_points)                
+        for id in range(len(interpolation_points)):
+            x = interpolation_points[id][0]
+            y = interpolation_points[id][1]
+
+            for i in range(20):
+                t = i*10
+                k = i%6
+
+                if k == 0:
+                    q0 = F(t, point_id=id)
+                    q1 = F(t+60, point_id=id)
+
+                if q0 == INF:
+                    actual = q0
+                else:
+                    actual = (k*q1 + (6-k)*q0)/6
+                q = F(t, point_id=id)
+                print "############"
+                print "id, x, y ", id, x, y #k, t, q
+                print "t", t
+                #print ' ', q0
+                #print ' ', q1
+                print "q",q
+                print "actual", actual
+                #print
+                if q0 == INF:
+                     self.failUnless( q == actual, 'Fail!')
+                else:
+                    assert allclose(q, actual)
+
+
+        #Another check of linear interpolation in time
+        for id in range(len(interpolation_points)):
+            q60 = F(60, point_id=id)
+            q120 = F(120, point_id=id)
+
+            t = 90 #Halfway between 60 and 120
+            q = F(t, point_id=id)
+            assert allclose( (q120+q60)/2, q )
+
+            t = 100 #Two thirds of the way between between 60 and 120
+            q = F(t, point_id=id)
+            assert allclose(q60/3 + 2*q120/3, q)
+
+
+
+        #Check that domain.starttime isn't updated if later than file starttime but earlier
+        #than file end time
+        delta = 23
+        domain.starttime = start + delta
+        F = file_function(filename + '.sww', domain,
+                          quantities = domain.conserved_quantities,
+                          interpolation_points = interpolation_points)
+        assert allclose(domain.starttime, start+delta)
+
+
+
+
+        #Now try interpolation with delta offset
+        for id in range(len(interpolation_points)):            
+            x = interpolation_points[id][0]
+            y = interpolation_points[id][1]
+
+            for i in range(20):
+                t = i*10
+                k = i%6
+
+                if k == 0:
+                    q0 = F(t-delta, point_id=id)
+                    q1 = F(t+60-delta, point_id=id)
+
+                q = F(t-delta, point_id=id)
+                assert allclose(q, (k*q1 + (6-k)*q0)/6)
+
+
+        os.remove(filename + '.sww')
 
     def test_file_function_time_with_domain(self):

inundation/pyvolution/util.py

@@ -307 +307 @@
 
     from least_squares import Interpolation_function
+    #from fit_interpolate.interpolate import Interpolation_function
 
     if not spatial: