Changeset 4312

Timestamp:

Mar 22, 2007, 12:16:57 PM (18 years ago)

Author:

ole

Message:

Created new boundary class: Field_boundary designed to replace the generic File_boundary. Field_boundary is specific to the shallow water wave equation and allows mean_stage to be specified adjusting the stage derived from the sww. This should address ticket:132.
On test was added, two updated.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/generic_boundary_conditions.py (modified) (3 diffs)
• abstract_2d_finite_volumes/util.py (modified) (1 diff)
• shallow_water/shallow_water_domain.py (modified) (1 diff)
• shallow_water/test_shallow_water_domain.py (modified) (3 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

@@ -180 +180 @@
     Note that the resulting solution history is not exactly the same as if
     the models were coupled as there is no feedback into the source model.
+       
     """
 
@@ -244 +245 @@
         self.F = file_function(filename, domain,
                                interpolation_points=self.midpoint_coordinates,
-                           time_thinning=time_thinning,
-                           use_cache=use_cache, 
-                           verbose=verbose)
+                               time_thinning=time_thinning,
+                               use_cache=use_cache, 
+                               verbose=verbose)
         self.domain = domain
 
@@ -264 +265 @@
     def evaluate(self, vol_id=None, edge_id=None):
         """Return linearly interpolated values based on domain.time
-
-        vol_id and edge_id are ignored
         """
 

anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

@@ -323 +323 @@
               'verbose': verbose
               }
+    
 #    print'CACHING FROM UTIL.py for interpolation_function', use_cache
-    from anuga.caching import myhash
+#    from anuga.caching import myhash
     if use_cache is True:
         from caching import cache

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -1392 +1392 @@
 
 
-
-
+class Field_boundary(Boundary):
+    """Set boundary from given field represented in an sww file containing values
+    for stage, xmomentum and ymomentum.
+    Optionally, the user can specify mean_stage to offset the stage provided in the
+    sww file.
+
+    This function is a thin wrapper around the generic File_boundary.
+    """
+
+
+    def __init__(self, filename, domain,
+                 mean_stage=0.0,
+                 time_thinning=1, 
+                 use_cache=False,
+                 verbose=False):
+        """Constructor
+
+        filename: Name of sww file
+        domain: pointer to shallow water domain for which the boundary applies
+        mean_stage: The mean water level which will be added to stage derived from the sww file
+        time_thinning:
+        use_cache:
+        verbose:
+        
+        """
+
+        # Create generic file_boundary object
+        self.file_boundary = File_boundary(filename, domain,
+                                           time_thinning=time_thinning,
+                                           use_cache=use_cache,
+                                           verbose=verbose)
+        
+        # Record information from File_boundary
+        self.F = self.file_boundary.F
+        self.domain = self.file_boundary.domain
+        
+        # Record mean stage
+        self.mean_stage = mean_stage
+
+
+    def __repr__(self):
+        return 'Field boundary'
+
+
+    def evaluate(self, vol_id=None, edge_id=None):
+        """Return linearly interpolated values based on domain.time
+
+        vol_id and edge_id are ignored
+        """
+
+        # Evaluate file boundary
+        q = self.file_boundary.evaluate(vol_id, edge_id)
+
+        # Adjust stage
+        for j, name in enumerate(self.domain.conserved_quantities):
+            if name == 'stage':
+                q[j] += self.mean_stage
+        return q
+
+    
 
 #########################

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -3718 +3718 @@
         #Bf = Spatio_temporal_boundary(domain1.get_name() + '.' +\
         #                              domain1.format, domain2)
-        Bf = File_boundary(domain1.get_name() + '.' +\
-                           domain1.format, domain2)        
+        Bf = Field_boundary(domain1.get_name() + '.' +\
+                            domain1.format, domain2)        
         domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
         domain2.check_integrity()
@@ -3859 +3859 @@
         #Bf = Spatio_temporal_boundary(domain1.get_name() + '.' + domain1.format,
         #                              domain2)
-        Bf = File_boundary(domain1.get_name() + '.' + domain1.format,
-                           domain2)        
+        Bf = Field_boundary(domain1.get_name() + '.' + domain1.format,
+                            domain2)        
         domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
         domain2.check_integrity()
@@ -3923 +3923 @@
         assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3)
 
+
+
+        #Cleanup
+        os.remove(domain1.get_name() + '.' + domain1.format)
+
+
+    def test_spatio_temporal_boundary_3(self):
+        """Test that boundary values can be read from file and interpolated
+        in both time and space.
+        This is a more basic test, verifying that boundary object
+        produces the expected results
+
+        This tests adjusting using mean_stage
+
+        """
+
+        import time
+
+        mean_stage = 5.2 # Adjust stage by this amount in boundary
+
+        #Create sww file of simple propagation from left to right
+        #through rectangular domain
+
+        from mesh_factory import rectangular
+
+        #Create basic mesh
+        points, vertices, boundary = rectangular(3, 3)
+
+        #Create shallow water domain
+        domain1 = Domain(points, vertices, boundary)
+
+        domain1.reduction = mean
+        domain1.smooth = True #To mimic MOST output
+
+        domain1.default_order = 2
+        domain1.store = True
+        domain1.set_datadir('.')
+        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
+
+        #FIXME: This is extremely important!
+        #How can we test if they weren't stored?
+        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
+
+
+        #Bed-slope and friction at vertices (and interpolated elsewhere)
+        domain1.set_quantity('elevation', 0)
+        domain1.set_quantity('friction', 0)
+
+        # Boundary conditions
+        Br = Reflective_boundary(domain1)
+        Bd = Dirichlet_boundary([0.3,0,0])
+        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
+        #Initial condition
+        domain1.set_quantity('stage', 0)
+        domain1.check_integrity()
+
+        finaltime = 5
+        #Evolution  (full domain - large steps)
+        for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
+            pass
+            #domain1.write_time()
+
+
+        #Create an triangle shaped domain (coinciding with some
+        #coordinates from domain 1),
+        #formed from the lower and right hand  boundaries and
+        #the sw-ne diagonal
+        #from domain 1. Call it domain2
+
+        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
+                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
+                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
+
+        vertices = [ [1,2,0],
+                     [3,4,1], [2,1,4], [4,5,2],
+                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
+
+        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
+                     (4,2):'right', (6,2):'right', (8,2):'right',
+                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
+
+        domain2 = Domain(points, vertices, boundary)
+
+        domain2.reduction = domain1.reduction
+        domain2.smooth = False
+        domain2.default_order = 2
+
+        #Bed-slope and friction at vertices (and interpolated elsewhere)
+        domain2.set_quantity('elevation', 0)
+        domain2.set_quantity('friction', 0)
+        domain2.set_quantity('stage', 0)
+
+
+        #Read results for specific timesteps t=1 and t=2
+        from Scientific.IO.NetCDF import NetCDFFile
+        fid = NetCDFFile(domain1.get_name() + '.' + domain1.format)
+
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        s1 = fid.variables['stage'][1,:]
+        s2 = fid.variables['stage'][2,:]
+        fid.close()
+
+        from Numeric import take, reshape, concatenate
+        shp = (len(x), 1)
+        points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)
+        #The diagonal points of domain 1 are 0, 5, 10, 15
+
+        #print points[0], points[5], points[10], points[15]
+        assert allclose( take(points, [0,5,10,15]),
+                         [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]])
+
+
+        # Boundary conditions
+        Br = Reflective_boundary(domain2)
+        #Bf = Spatio_temporal_boundary(domain1.get_name() + '.' + domain1.format,
+        #                              domain2)
+        Bf = Field_boundary(domain1.get_name() + '.' + domain1.format,
+                            domain2, mean_stage=mean_stage)
+        
+        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
+        domain2.check_integrity()
+
+        #Test that interpolation points are the mid points of the all boundary
+        #segments
+
+        boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0],
+                              [1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6],
+                              [1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]]
+
+        boundary_midpoints.sort()
+        R = Bf.F.interpolation_points.tolist()
+        R.sort()
+        assert allclose(boundary_midpoints, R)
+
+        #Check spatially interpolated output at time == 1
+        domain2.time = 1
+
+        #First diagonal midpoint
+        R0 = Bf.evaluate(0,0)
+        assert allclose(R0[0], (s1[0] + s1[5])/2 + mean_stage)
+
+        #Second diagonal midpoint
+        R0 = Bf.evaluate(3,0)
+        assert allclose(R0[0], (s1[5] + s1[10])/2 + mean_stage)
+
+        #First diagonal midpoint
+        R0 = Bf.evaluate(8,0)
+        assert allclose(R0[0], (s1[10] + s1[15])/2 + mean_stage)
+
+        #Check spatially interpolated output at time == 2
+        domain2.time = 2
+
+        #First diagonal midpoint
+        R0 = Bf.evaluate(0,0)
+        assert allclose(R0[0], (s2[0] + s2[5])/2 + mean_stage)
+
+        #Second diagonal midpoint
+        R0 = Bf.evaluate(3,0)
+        assert allclose(R0[0], (s2[5] + s2[10])/2 + mean_stage)
+
+        #First diagonal midpoint
+        R0 = Bf.evaluate(8,0)
+        assert allclose(R0[0], (s2[10] + s2[15])/2 + mean_stage)
+
+
+        #Now check temporal interpolation
+
+        domain2.time = 1 + 2.0/3
+
+        #First diagonal midpoint
+        R0 = Bf.evaluate(0,0)
+        assert allclose(R0[0], ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3 + mean_stage)
+
+        #Second diagonal midpoint
+        R0 = Bf.evaluate(3,0)
+        assert allclose(R0[0], ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3 + mean_stage)
+
+        #First diagonal midpoint
+        R0 = Bf.evaluate(8,0)
+        assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3 + mean_stage)