Changeset 5672

Timestamp:

Aug 21, 2008, 4:32:14 PM (16 years ago)

Author:

ole

Message:

Tested for bug where default_boundry was repeatedly added to by Field_boundary.
Fixed bug and updated manual accordingly.

This closes ticket:293

Files:

• anuga_core/documentation/user_manual/anuga_user_manual.tex (modified) (1 diff)
• anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py (modified) (2 diffs)
• anuga_core/source/anuga/shallow_water/shallow_water_domain.py (modified) (1 diff)
• anuga_core/source/anuga/shallow_water/test_data_manager.py (modified) (2 diffs)
• anuga_work/development/boxingday08/project.py (modified) (1 diff)
• anuga_work/development/boxingday08/run_boxingday_polyline.py (modified) (3 diffs)

anuga_core/documentation/user_manual/anuga_user_manual.tex

@@ -2254 +2254 @@
 appropriate segments for each conserved quantity.
 
-Optional argument \code{default_boundary} can be used to specify another boundary object to be used in case model time exceeds the time availabel in the file used by \code{File\_boundary}.
-The default_boundary could be a simple Dirichlet condition or even another File\_boundary 
+Optional argument \code{default\_boundary} can be used to specify another boundary object to be used in case model time exceeds the time availabel in the file used by \code{File\_boundary}.
+The \code{default\_boundary} could be a simple Dirichlet condition or 
+even another \code{File\_boundary} 
 typically using data pertaining to another time interval.  
 \end{classdesc}
 
-
+\begin{classdesc}{Field\_boundary}{Boundary}
+Module: \module{shallow\_water.shallow\_water\_domain}
+
+This method works in the same way as \code{File\_boundary} except that it 
+allows for the value of stage to be offset by a constant specified in the 
+keyword argument \code{mean\_stage}.
+
+This functionality allows for models to be run for a range of tides using
+the same boundary information (from .sts, .sww or .tms files). The tidal value
+for each run would then be specified in the keyword argument 
+\code{mean\_stage}. 
+If \code{mean\_stage} = 0.0, \code{Field\_boundary} and \code{File\_boundary} 
+behave identically. 
+
+
+Note that if the optional argument \code{default\_boundary} is specified
+it's stage value will be adjusted by \code{mean\_stage} just like the values
+obtained from the file. 
+
+See \code{File\_boundary} for further details.
+\end{classdesc}
 
 %%%

anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

@@ -332 +332 @@
                     # Pass control to default boundary
                     res = self.default_boundary.evaluate(vol_id, edge_id)
-                                    
+                    
+                    # Ensure that result cannot be manipulated
+                    # This is a real danger in case the 
+                    # default_boundary is a Dirichlet type 
+                    # for instance. 
+                    res = res.copy() 
+                    
                     if self.default_boundary_invoked is False:
                         # Issue warning the first time
@@ -347 +353 @@
                     
 
-                    
-                
             if res == NAN:
                 x,y=self.midpoint_coordinates[i,:]

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -1145 +1145 @@
         vol_id and edge_id are ignored
         """
-
+        
         # Evaluate file boundary
         q = self.file_boundary.evaluate(vol_id, edge_id)

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -7801 +7801 @@
         from anuga.pmesh.mesh_interface import create_mesh_from_regions
 
-        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]]
+        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],
+                          [6.02,97.02],[6.00,97.02]]
         tide = 0.37
         time_step_count = 5
         time_step = 2
-        lat_long_points =bounding_polygon[0:3]
+        lat_long_points = bounding_polygon[0:3]
         n=len(lat_long_points)
         first_tstep=ones(n,Int)
@@ -7936 +7937 @@
         os.remove(meshname)
 
-        
+
+        
+                
+        
+    def test_field_boundary_stsI_beyond_model_time(self):
+        """test_field_boundary(self):
+        
+        Test that field_boundary can work when model time
+        exceeds available data whilst adjusting mean_stage.
+        
+        """
+        
+        # Don't do warnings in unit test
+        import warnings
+        warnings.simplefilter('ignore')
+        
+        from anuga.shallow_water import Domain
+        from anuga.shallow_water import Reflective_boundary
+        from anuga.shallow_water import Dirichlet_boundary
+        from anuga.shallow_water import File_boundary
+        from anuga.pmesh.mesh_interface import create_mesh_from_regions
+
+        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],
+                          [6.02,97.02],[6.00,97.02]]
+        tide = 0.37
+        time_step_count = 5
+        time_step = 2
+        lat_long_points = bounding_polygon[0:3]
+        n=len(lat_long_points)
+        first_tstep=ones(n,Int)
+        last_tstep=(time_step_count)*ones(n,Int)
+
+        h = 20        
+        w = 2
+        u = 10
+        v = -10
+        gauge_depth=h*ones(n,Float)
+        ha=w*ones((n,time_step_count),Float)
+        ua=u*ones((n,time_step_count),Float)
+        va=v*ones((n,time_step_count),Float)
+        base_name, files = self.write_mux2(lat_long_points,
+                                           time_step_count, time_step,
+                                           first_tstep, last_tstep,
+                                           depth=gauge_depth,
+                                           ha=ha,
+                                           ua=ua,
+                                           va=va)
+
+        sts_file=base_name
+        urs2sts(base_name,
+                sts_file,
+                mean_stage=0.0, # Deliberately let Field_boundary do the adjustment
+                verbose=False)
+        self.delete_mux(files)
+
+        #print 'start create mesh from regions'
+        for i in range(len(bounding_polygon)):
+            zone,\
+            bounding_polygon[i][0],\
+            bounding_polygon[i][1]=redfearn(bounding_polygon[i][0],
+                                            bounding_polygon[i][1])
+                                            
+        extent_res=1000000
+        meshname = 'urs_test_mesh' + '.tsh'
+        interior_regions=None
+        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
+        create_mesh_from_regions(bounding_polygon,
+                                 boundary_tags=boundary_tags,
+                                 maximum_triangle_area=extent_res,
+                                 filename=meshname,
+                                 interior_regions=interior_regions,
+                                 verbose=False)
+        
+        domain_fbound = Domain(meshname)
+        domain_fbound.set_quantity('stage', tide)
+        
+        Br = Reflective_boundary(domain_fbound)
+        Bd = Dirichlet_boundary([w+tide, u*(w+h), v*(w+h)])
+        Bdefault = Dirichlet_boundary([w, u*(w+h), v*(w+h)])        
+                
+        Bf = Field_boundary(sts_file+'.sts', 
+                           domain_fbound, 
+                           mean_stage=tide, # Field boundary to adjust for tide
+                           boundary_polygon=bounding_polygon,
+                           default_boundary=Bdefault) # Condition to be used 
+                                                      # if model time exceeds 
+                                                      # available data
+
+        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
+        
+        # Check boundary object evaluates as it should
+        for i, ((vol_id, edge_id), B) in enumerate(domain_fbound.boundary_objects):
+            if B is Bf:
+            
+                qf = B.evaluate(vol_id, edge_id)  # Field boundary
+                qd = Bd.evaluate(vol_id, edge_id) # Dirichlet boundary
+                
+                msg = 'Got %s, should have been %s' %(qf, qd)
+                assert allclose(qf, qd), msg 
+                
+        # Evolve
+        data_finaltime = time_step*(time_step_count-1)
+        finaltime = data_finaltime + 10 # Let model time exceed available data
+        yieldstep = time_step
+        temp_fbound=zeros(int(finaltime/yieldstep)+1, Float)
+
+        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,
+                                                   finaltime=finaltime, 
+                                                   skip_initial_step=False)):
+                                                   
+            D = domain_fbound
+            temp_fbound[i]=D.quantities['stage'].centroid_values[2]
+
+            # Check that file boundary object has populated 
+            # boundary array correctly  
+            # FIXME (Ole): Do this for the other tests too!
+            for j, val in enumerate(D.get_quantity('stage').boundary_values):
+            
+                (vol_id, edge_id), B = D.boundary_objects[j]
+                if isinstance(B, Field_boundary):
+                    msg = 'Got %f should have been %f' %(val, w+tide)
+                    assert allclose(val, w + tide), msg
+
 
     def test_file_boundary_stsII(self):

anuga_work/development/boxingday08/project.py

@@ -99 +99 @@
 patong_bay_polygon = read_polygon(dir+sep+'patong_bay_polygon.csv')
     
-tide=0.45
+tide=0.15
 
 ######################################

anuga_work/development/boxingday08/run_boxingday_polyline.py

@@ -54 +54 @@
 timestamp = strftime('%Y%m%d_%H%M%S',localtime())
 basename = scenario[:4] + '_polyline'
+tide = project.tide
 
 
@@ -130 +131 @@
 print domain.statistics()
 
-domain.set_name(basename+timestamp)
+domain.set_name(basename+timestamp+'_'+str(tide))
 #domain.set_name(basename)
 domain.set_datadir(scenario)
@@ -146 +147 @@
 #------------------------------------------------------------------------------
      
-tide = project.tide
 domain.set_quantity('stage', tide)
 domain.set_quantity('friction', 0.01)