Changeset 4312


Ignore:
Timestamp:
Mar 22, 2007, 12:16:57 PM (17 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:
4 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

    r4293 r4312  
    180180    Note that the resulting solution history is not exactly the same as if
    181181    the models were coupled as there is no feedback into the source model.
     182       
    182183    """
    183184
     
    244245        self.F = file_function(filename, domain,
    245246                               interpolation_points=self.midpoint_coordinates,
    246                            time_thinning=time_thinning,
    247                            use_cache=use_cache,
    248                            verbose=verbose)
     247                               time_thinning=time_thinning,
     248                               use_cache=use_cache,
     249                               verbose=verbose)
    249250        self.domain = domain
    250251
     
    264265    def evaluate(self, vol_id=None, edge_id=None):
    265266        """Return linearly interpolated values based on domain.time
    266 
    267         vol_id and edge_id are ignored
    268267        """
    269268
  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r4306 r4312  
    323323              'verbose': verbose
    324324              }
     325   
    325326#    print'CACHING FROM UTIL.py for interpolation_function', use_cache
    326     from anuga.caching import myhash
     327#    from anuga.caching import myhash
    327328    if use_cache is True:
    328329        from caching import cache
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r4258 r4312  
    13921392
    13931393
    1394 
    1395 
     1394class Field_boundary(Boundary):
     1395    """Set boundary from given field represented in an sww file containing values
     1396    for stage, xmomentum and ymomentum.
     1397    Optionally, the user can specify mean_stage to offset the stage provided in the
     1398    sww file.
     1399
     1400    This function is a thin wrapper around the generic File_boundary.
     1401    """
     1402
     1403
     1404    def __init__(self, filename, domain,
     1405                 mean_stage=0.0,
     1406                 time_thinning=1,
     1407                 use_cache=False,
     1408                 verbose=False):
     1409        """Constructor
     1410
     1411        filename: Name of sww file
     1412        domain: pointer to shallow water domain for which the boundary applies
     1413        mean_stage: The mean water level which will be added to stage derived from the sww file
     1414        time_thinning:
     1415        use_cache:
     1416        verbose:
     1417       
     1418        """
     1419
     1420        # Create generic file_boundary object
     1421        self.file_boundary = File_boundary(filename, domain,
     1422                                           time_thinning=time_thinning,
     1423                                           use_cache=use_cache,
     1424                                           verbose=verbose)
     1425       
     1426        # Record information from File_boundary
     1427        self.F = self.file_boundary.F
     1428        self.domain = self.file_boundary.domain
     1429       
     1430        # Record mean stage
     1431        self.mean_stage = mean_stage
     1432
     1433
     1434    def __repr__(self):
     1435        return 'Field boundary'
     1436
     1437
     1438    def evaluate(self, vol_id=None, edge_id=None):
     1439        """Return linearly interpolated values based on domain.time
     1440
     1441        vol_id and edge_id are ignored
     1442        """
     1443
     1444        # Evaluate file boundary
     1445        q = self.file_boundary.evaluate(vol_id, edge_id)
     1446
     1447        # Adjust stage
     1448        for j, name in enumerate(self.domain.conserved_quantities):
     1449            if name == 'stage':
     1450                q[j] += self.mean_stage
     1451        return q
     1452
     1453   
    13961454
    13971455#########################
  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r4303 r4312  
    37183718        #Bf = Spatio_temporal_boundary(domain1.get_name() + '.' +\
    37193719        #                              domain1.format, domain2)
    3720         Bf = File_boundary(domain1.get_name() + '.' +\
    3721                            domain1.format, domain2)       
     3720        Bf = Field_boundary(domain1.get_name() + '.' +\
     3721                            domain1.format, domain2)       
    37223722        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
    37233723        domain2.check_integrity()
     
    38593859        #Bf = Spatio_temporal_boundary(domain1.get_name() + '.' + domain1.format,
    38603860        #                              domain2)
    3861         Bf = File_boundary(domain1.get_name() + '.' + domain1.format,
    3862                            domain2)       
     3861        Bf = Field_boundary(domain1.get_name() + '.' + domain1.format,
     3862                            domain2)       
    38633863        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
    38643864        domain2.check_integrity()
     
    39233923        assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3)
    39243924
     3925
     3926
     3927        #Cleanup
     3928        os.remove(domain1.get_name() + '.' + domain1.format)
     3929
     3930
     3931    def test_spatio_temporal_boundary_3(self):
     3932        """Test that boundary values can be read from file and interpolated
     3933        in both time and space.
     3934        This is a more basic test, verifying that boundary object
     3935        produces the expected results
     3936
     3937        This tests adjusting using mean_stage
     3938
     3939        """
     3940
     3941        import time
     3942
     3943        mean_stage = 5.2 # Adjust stage by this amount in boundary
     3944
     3945        #Create sww file of simple propagation from left to right
     3946        #through rectangular domain
     3947
     3948        from mesh_factory import rectangular
     3949
     3950        #Create basic mesh
     3951        points, vertices, boundary = rectangular(3, 3)
     3952
     3953        #Create shallow water domain
     3954        domain1 = Domain(points, vertices, boundary)
     3955
     3956        domain1.reduction = mean
     3957        domain1.smooth = True #To mimic MOST output
     3958
     3959        domain1.default_order = 2
     3960        domain1.store = True
     3961        domain1.set_datadir('.')
     3962        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
     3963
     3964        #FIXME: This is extremely important!
     3965        #How can we test if they weren't stored?
     3966        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
     3967
     3968
     3969        #Bed-slope and friction at vertices (and interpolated elsewhere)
     3970        domain1.set_quantity('elevation', 0)
     3971        domain1.set_quantity('friction', 0)
     3972
     3973        # Boundary conditions
     3974        Br = Reflective_boundary(domain1)
     3975        Bd = Dirichlet_boundary([0.3,0,0])
     3976        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
     3977        #Initial condition
     3978        domain1.set_quantity('stage', 0)
     3979        domain1.check_integrity()
     3980
     3981        finaltime = 5
     3982        #Evolution  (full domain - large steps)
     3983        for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
     3984            pass
     3985            #domain1.write_time()
     3986
     3987
     3988        #Create an triangle shaped domain (coinciding with some
     3989        #coordinates from domain 1),
     3990        #formed from the lower and right hand  boundaries and
     3991        #the sw-ne diagonal
     3992        #from domain 1. Call it domain2
     3993
     3994        points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
     3995                   [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
     3996                   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
     3997
     3998        vertices = [ [1,2,0],
     3999                     [3,4,1], [2,1,4], [4,5,2],
     4000                     [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]
     4001
     4002        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
     4003                     (4,2):'right', (6,2):'right', (8,2):'right',
     4004                     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}
     4005
     4006        domain2 = Domain(points, vertices, boundary)
     4007
     4008        domain2.reduction = domain1.reduction
     4009        domain2.smooth = False
     4010        domain2.default_order = 2
     4011
     4012        #Bed-slope and friction at vertices (and interpolated elsewhere)
     4013        domain2.set_quantity('elevation', 0)
     4014        domain2.set_quantity('friction', 0)
     4015        domain2.set_quantity('stage', 0)
     4016
     4017
     4018        #Read results for specific timesteps t=1 and t=2
     4019        from Scientific.IO.NetCDF import NetCDFFile
     4020        fid = NetCDFFile(domain1.get_name() + '.' + domain1.format)
     4021
     4022        x = fid.variables['x'][:]
     4023        y = fid.variables['y'][:]
     4024        s1 = fid.variables['stage'][1,:]
     4025        s2 = fid.variables['stage'][2,:]
     4026        fid.close()
     4027
     4028        from Numeric import take, reshape, concatenate
     4029        shp = (len(x), 1)
     4030        points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)
     4031        #The diagonal points of domain 1 are 0, 5, 10, 15
     4032
     4033        #print points[0], points[5], points[10], points[15]
     4034        assert allclose( take(points, [0,5,10,15]),
     4035                         [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]])
     4036
     4037
     4038        # Boundary conditions
     4039        Br = Reflective_boundary(domain2)
     4040        #Bf = Spatio_temporal_boundary(domain1.get_name() + '.' + domain1.format,
     4041        #                              domain2)
     4042        Bf = Field_boundary(domain1.get_name() + '.' + domain1.format,
     4043                            domain2, mean_stage=mean_stage)
     4044       
     4045        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
     4046        domain2.check_integrity()
     4047
     4048        #Test that interpolation points are the mid points of the all boundary
     4049        #segments
     4050
     4051        boundary_midpoints = [[1.0/6, 0], [1.0/2, 0], [5.0/6,0],
     4052                              [1.0, 1.0/6], [1.0, 1.0/2], [1.0, 5.0/6],
     4053                              [1.0/6, 1.0/6], [0.5, 0.5], [5.0/6, 5.0/6]]
     4054
     4055        boundary_midpoints.sort()
     4056        R = Bf.F.interpolation_points.tolist()
     4057        R.sort()
     4058        assert allclose(boundary_midpoints, R)
     4059
     4060        #Check spatially interpolated output at time == 1
     4061        domain2.time = 1
     4062
     4063        #First diagonal midpoint
     4064        R0 = Bf.evaluate(0,0)
     4065        assert allclose(R0[0], (s1[0] + s1[5])/2 + mean_stage)
     4066
     4067        #Second diagonal midpoint
     4068        R0 = Bf.evaluate(3,0)
     4069        assert allclose(R0[0], (s1[5] + s1[10])/2 + mean_stage)
     4070
     4071        #First diagonal midpoint
     4072        R0 = Bf.evaluate(8,0)
     4073        assert allclose(R0[0], (s1[10] + s1[15])/2 + mean_stage)
     4074
     4075        #Check spatially interpolated output at time == 2
     4076        domain2.time = 2
     4077
     4078        #First diagonal midpoint
     4079        R0 = Bf.evaluate(0,0)
     4080        assert allclose(R0[0], (s2[0] + s2[5])/2 + mean_stage)
     4081
     4082        #Second diagonal midpoint
     4083        R0 = Bf.evaluate(3,0)
     4084        assert allclose(R0[0], (s2[5] + s2[10])/2 + mean_stage)
     4085
     4086        #First diagonal midpoint
     4087        R0 = Bf.evaluate(8,0)
     4088        assert allclose(R0[0], (s2[10] + s2[15])/2 + mean_stage)
     4089
     4090
     4091        #Now check temporal interpolation
     4092
     4093        domain2.time = 1 + 2.0/3
     4094
     4095        #First diagonal midpoint
     4096        R0 = Bf.evaluate(0,0)
     4097        assert allclose(R0[0], ((s1[0] + s1[5])/2 + 2.0*(s2[0] + s2[5])/2)/3 + mean_stage)
     4098
     4099        #Second diagonal midpoint
     4100        R0 = Bf.evaluate(3,0)
     4101        assert allclose(R0[0], ((s1[5] + s1[10])/2 + 2.0*(s2[5] + s2[10])/2)/3 + mean_stage)
     4102
     4103        #First diagonal midpoint
     4104        R0 = Bf.evaluate(8,0)
     4105        assert allclose(R0[0], ((s1[10] + s1[15])/2 + 2.0*(s2[10] + s2[15])/2)/3 + mean_stage)
    39254106
    39264107
Note: See TracChangeset for help on using the changeset viewer.