Changeset 8486


Ignore:
Timestamp:
Jul 31, 2012, 8:37:12 PM (13 years ago)
Author:
steve
Message:

Changing calculation of boundary conditions

Location:
trunk/anuga_core
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/documentation/user_manual/demos/channel_variable.py

    r7376 r8486  
    77# Import necessary modules
    88#------------------------------------------------------------------------------
    9 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
    10 from anuga.shallow_water import Domain
    11 from anuga.shallow_water import Reflective_boundary
    12 from anuga.shallow_water import Dirichlet_boundary
    13 from anuga.shallow_water import Time_boundary
     9from anuga import rectangular_cross
     10from anuga import Domain
     11from anuga import Reflective_boundary
     12from anuga import Dirichlet_boundary
     13from anuga import Time_boundary
    1414
    1515#------------------------------------------------------------------------------

  • trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

    r8124 r8486  
    1 """boundary.py - Classes for implementing boundary conditions
     1
     2"""
     3boundary.py - Classes for implementing boundary conditions
    24"""
    35
    … …  
    2931
    3032
     33    def evaluate_segment(self, domain=None, segment_edges=None):
     34        """
     35        Evaluate boundary condition at edges of a domain in a list
     36        defined by segment_edges
     37
     38        Go through list of boundary objects and update boundary values
     39        for all conserved quantities on boundary.
     40        It is assumed that the ordering of conserved quantities is
     41        consistent between the domain and the boundary object, i.e.
     42        the jth element of vector q must correspond to the jth conserved
     43        quantity in domain.
     44        """
     45
     46        if segment_edges is None:
     47            return
     48        if domain is None:
     49            return
     50       
     51        for i in segment_edges:
     52            vol_id  = domain.boundary_cells[i]
     53            edge_id = domain.boundary_edges[i]
     54
     55
     56            q_bdry = self.evaluate(vol_id, edge_id)
     57
     58            if len(q_bdry) == len(domain.evolved_quantities):
     59                # conserved and evolved quantities are the same
     60                q_evol = q_bdry
     61            elif len(q_bdry) == len(domain.conserved_quantities):
     62                # boundary just returns conserved quantities
     63                # Need to calculate all the evolved quantities
     64                # Use default conversion
     65
     66                q_evol = domain.get_evolved_quantities(vol_id, edge = edge_id)
     67
     68                q_evol = domain.conserved_values_to_evolved_values \
     69                                                        (q_bdry, q_evol)
     70            else:
     71                msg = 'Boundary must return array of either conserved'
     72                msg += ' or evolved quantities'
     73                raise Exception(msg)
     74
     75            for j, name in enumerate(domain.evolved_quantities):
     76                Q = domain.quantities[name]
     77                Q.boundary_values[i] = q_evol[j]
     78
     79
    3180class Transmissive_boundary(Boundary):
    3281    """Transmissive boundary returns same conserved quantities as
    … …  
    62111
    63112
     113
     114    def evaluate_segment(self, domain, segment_edges):
     115
     116        if segment_edges is None:
     117            return
     118        if domain is None:
     119            return
     120
     121
     122        ids = segment_edges
     123        vol_ids  = domain.boundary_cells[ids]
     124        edge_ids = domain.boundary_edges[ids]
     125
     126
     127        if domain.centroid_transmissive_bc :
     128            for j, name in enumerate(domain.evolved_quantities):
     129                Q = domain.quantities[name]
     130                Q.boundary_values[ids] = Q.centroid_values[vol_ids]
     131        else:
     132            for j, name in enumerate(domain.evolved_quantities):
     133                Q = domain.quantities[name]
     134                Q.boundary_values[ids] = Q.edge_values[vol_ids,edge_ids]
     135
     136
     137
    64138class Dirichlet_boundary(Boundary):
    65139    """Dirichlet boundary returns constant values for the
    … …  
    68142
    69143
    70     def __init__(self, conserved_quantities=None):
     144    def __init__(self, dirichlet_values=None):
    71145        Boundary.__init__(self)
    72146
    73         if conserved_quantities is None:
    74             msg = 'Must specify one value for each conserved quantity'
     147        if dirichlet_values is None:
     148            msg = 'Must specify one value for each quantity'
    75149            raise Exception(msg)
    76150
    77         self.conserved_quantities=num.array(conserved_quantities, num.float)
     151        self.dirichlet_values=num.array(dirichlet_values, num.float)
     152
     153
    78154
    79155    def __repr__(self):
    80         return 'Dirichlet boundary (%s)' %self.conserved_quantities
     156        return 'Dirichlet boundary (%s)' %self.dirichlet_values
     157
    81158
    82159    def evaluate(self, vol_id=None, edge_id=None):
    83         return self.conserved_quantities
     160        return self.dirichlet_values
     161
     162    def evaluate_segment(self, domain, segment_edges):
     163
     164        if segment_edges is None:
     165            return
     166        if domain is None:
     167            return
     168
     169
     170        ids = segment_edges
     171
     172
     173        vol_ids  = domain.boundary_cells[ids]
     174        edge_ids = domain.boundary_edges[ids]
     175
     176        q_bdry = self.dirichlet_values
     177
     178        conserved_quantities = True
     179        if len(q_bdry) == len(domain.evolved_quantities):
     180            # enough dirichlet values to set evolved quantities
     181            conserved_quantities = False
     182
     183        #--------------------------------------------------
     184        # First populate all the boundary values with
     185        # interior edge values
     186        #--------------------------------------------------
     187        if  conserved_quantities:
     188            for j, name in enumerate(domain.evolved_quantities):
     189                Q = domain.quantities[name]
     190                Q.boundary_values[ids] = Q.edge_values[vol_ids,edge_ids]
     191
     192        #--------------------------------------------------
     193        # Now over write with constant Dirichlet value
     194        #--------------------------------------------------
     195        if conserved_quantities:
     196            quantities = domain.conserved_quantities
     197        else:
     198            quantities = domain.evolved_quantities
     199
     200        #-------------------------------------------------
     201        # Now update to Dirichlet values
     202        #-------------------------------------------------
     203        for j, name in enumerate(quantities):
     204            Q = domain.quantities[name]
     205            Q.boundary_values[ids] = q_bdry[j]
    84206
    85207
    … …  
    101223    """
    102224
    103     # FIXME (Ole): We should rename f to function to be consistent with
    104     # Transmissive_Momentum_Set_Stage_Boundary (cf posting by rrraman)
    105225    def __init__(self, domain=None,
    106                  f=None, # Should be removed and replaced by function below
     226                 #f=None, # Should be removed and replaced by function below
    107227                 function=None,
    108228                 default_boundary=None,
    … …  
    117237        if domain is None:
    118238            raise Exception('You must specify a domain to Time_boundary')
    119 
    120         # FIXME: Temporary code to deal with both f and function
    121         if function is not None and f is not None:
    122             raise Exception('Specify either function or f to Time_boundary')
    123            
    124         if function is None and f is None:
     239           
     240        if function is None:
    125241            raise Exception('You must specify a function to Time_boundary')
    126242           
    127         if f is None:
    128             f = function
    129         #####
    130243       
    131244        try:
    132             q = f(0.0)
     245            q = function(0.0)
    133246        except Exception, e:
    134247            msg = 'Function for time boundary could not be executed:\n%s' %e
    … …  
    152265        assert len(q) == d, msg
    153266
    154         self.f = f
     267        self.f = function
    155268        self.domain = domain
    156269
    … …  
    158271        return 'Time boundary'
    159272
     273    def get_time(self):
     274
     275        return self.domain.get_time()
     276
    160277    def evaluate(self, vol_id=None, edge_id=None):
    161         # FIXME (Ole): I think this should be get_time(), see ticket:306
     278
     279        return self.get_boundary_values()
     280
     281
     282    def evaluate_segment(self, domain, segment_edges):
     283
     284        if segment_edges is None:
     285            return
     286        if domain is None:
     287            return
     288
     289        ids = segment_edges
     290
     291        vol_ids  = domain.boundary_cells[ids]
     292        edge_ids = domain.boundary_edges[ids]
     293
     294        q_bdry = self.get_boundary_values()
     295
     296        conserved_quantities = True
     297        if len(q_bdry) == len(domain.evolved_quantities):
     298            # enough dirichlet values to set evolved quantities
     299            conserved_quantities = False
     300
     301        #--------------------------------------------------
     302        # First populate all the boundary values with
     303        # interior edge values
     304        #--------------------------------------------------
     305        if  conserved_quantities:
     306            for j, name in enumerate(domain.evolved_quantities):
     307                Q = domain.quantities[name]
     308                Q.boundary_values[ids] = Q.edge_values[vol_ids,edge_ids]
     309
     310        #--------------------------------------------------
     311        # Now over write with constant Dirichlet value
     312        #--------------------------------------------------
     313        if conserved_quantities:
     314            quantities = domain.conserved_quantities
     315        else:
     316            quantities = domain.evolved_quantities
     317
     318        #-------------------------------------------------
     319        # Now update to Dirichlet values
     320        #-------------------------------------------------
     321        for j, name in enumerate(quantities):
     322            Q = domain.quantities[name]
     323            Q.boundary_values[ids] = q_bdry[j]
     324           
     325
     326    def get_boundary_values(self, t=None):
     327
     328        if t is None:
     329            t = self.get_time()
     330           
    162331        try:
    163             res = self.f(self.domain.time)
     332            res = self.f(t)
    164333        except Modeltime_too_early, e:
    165334            raise Modeltime_too_early(e)
    … …  
    169338            else:
    170339                # Pass control to default boundary
    171                 res = self.default_boundary.evaluate(vol_id, edge_id)
     340                res = self.default_boundary
    172341               
    173342                # Ensure that result cannot be manipulated
    … …  
    192361
    193362        return res
     363
     364
     365
     366
    194367
    195368class Time_space_boundary(Boundary):

  • trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_domain.py

    r8456 r8486  
    2424
    2525import numpy as num
     26
     27
     28
     29
    2630
    2731class Generic_Domain:
    … …  
    13481352        self.update_ghosts()
    13491353
     1354
    13501355        # Initial update of vertex and edge values
    13511356        self.distribute_to_vertices_and_edges()
    13521357
     1358
    13531359        # Initial update boundary values
    13541360        self.update_boundary()
    … …  
    13571363        # Update extrema if necessary (for reporting)
    13581364        self.update_extrema()
     1365
     1366
    13591367
    13601368        # Or maybe restore from latest checkpoint
    … …  
    14091417            if self._order_ == 1:
    14101418                self.number_of_first_order_steps += 1
     1419
     1420           
    14111421
    14121422            # Yield results
    … …  
    18191829           
    18201830            B = self.boundary_map[tag]
    1821            
    1822             #if B is None:
    1823             #        log.critical('WARNING: Ignored boundary segment (None)')
    1824 
    1825             for i in self.tag_boundary_cells[tag]:
    1826                 vol_id  = self.boundary_cells[i]
    1827                 edge_id = self.boundary_edges[i]
    1828 
    1829                 if B is None:
    1830                     pass
    1831                 else:
    1832                     q_bdry = B.evaluate(vol_id, edge_id)
    1833 
    1834                     if len(q_bdry) == len(self.evolved_quantities):
    1835                         # conserved and evolved quantities are the same
    1836                         q_evol = q_bdry
    1837                     elif len(q_bdry) == len(self.conserved_quantities):
    1838                         # boundary just returns conserved quantities
    1839                         # Need to calculate all the evolved quantities
    1840                         # Use default conversion
    1841 
    1842                         q_evol = self.get_evolved_quantities(vol_id, edge = edge_id)
    1843 
    1844                         q_evol = self.conserved_values_to_evolved_values \
    1845                                                                 (q_bdry, q_evol)
    1846                     else:
    1847                         msg = 'Boundary must return array of either conserved'
    1848                         msg += ' or evolved quantities'
    1849                         raise Exception(msg)
    1850 
    1851                     for j, name in enumerate(self.evolved_quantities):
    1852                         Q = self.quantities[name]
    1853                         Q.boundary_values[i] = q_evol[j]
     1831
     1832            if B is None:
     1833                continue
     1834
     1835            boundary_segment_edges = self.tag_boundary_cells[tag]
     1836
     1837            B.evaluate_segment(self, boundary_segment_edges)
     1838       
    18541839
    18551840

  • trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

    r8455 r8486  
    13581358        self.y_gradient *= 0.0
    13591359
    1360     def get_integral(self):
     1360    def get_integral(self, full_only=False):
     1361
    13611362        """Compute the integral of quantity across entire domain."""
    13621363
    13631364       
    13641365        areas = self.domain.get_areas()
    1365         """
    1366         integral = 0
    1367         for k in range(len(self.domain)):
    1368             area = areas[k]
    1369             qc = self.centroid_values[k]
    1370             integral += qc*area
    1371         """
    1372         return num.sum(areas*self.centroid_values)
    1373 
    1374         #return integral
     1366
     1367        if full_only:
     1368            indices = num.where(self.domain.tri_full_flag ==1)[0]
     1369            return num.sum(areas[indices]*self.centroid_values[indices])
     1370        else:
     1371            return num.sum(areas*self.centroid_values)
     1372
    13751373
    13761374    def get_gradients(self):

  • trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py

    r8420 r8486  
    594594                              'Second': anuga.Transmissive_boundary(domain)} )
    595595
    596         try:
    597             domain.update_boundary()
    598         except:
    599             pass
    600         else:
    601             msg = 'Should have caught the lack of conserved_values_to_evolved_values member function'
    602             raise Exception, msg
    603 
     596#        try:
     597#            domain.update_boundary()
     598#        except:
     599#            pass
     600#        else:
     601#            msg = 'Should have caught the lack of conserved_values_to_evolved_values member function'
     602#            raise Exception, msg
     603
     604        domain.update_boundary()
    604605
    605606        def  conserved_values_to_evolved_values(q_cons, q_evol):

  • trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/test_gauge.py

    r8204 r8486  
    135135            line.append([float(row[0]),float(row[1]),float(row[2]),float(row[3]),
    136136                         float(row[4]),float(row[5]),float(row[6])])
    137 #            print 'assert line',line[i],'point1',point1_answers_array[i]
     137            #print 'assert line',line[i],'point1',point1_answers_array[i]
    138138            assert num.allclose(line[i], point1_answers_array[i])
    139139

  • trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/test_pmesh2domain.py

    r7737 r8486  
    6565
    6666         tags = {}
    67          b1 =  Dirichlet_boundary(conserved_quantities = num.array([0.0]))
    68          b2 =  Dirichlet_boundary(conserved_quantities = num.array([1.0]))
    69          b3 =  Dirichlet_boundary(conserved_quantities = num.array([2.0]))
     67         b1 =  Dirichlet_boundary(dirichlet_values = num.array([0.0]))
     68         b2 =  Dirichlet_boundary(dirichlet_values = num.array([1.0]))
     69         b3 =  Dirichlet_boundary(dirichlet_values = num.array([2.0]))
    7070         tags["1"] = b1
    7171         tags["2"] = b2
    … …  
    159159       
    160160         tags = {}
    161          b1 =  Dirichlet_boundary(conserved_quantities = num.array([0.0]))
    162          b2 =  Dirichlet_boundary(conserved_quantities = num.array([1.0]))
    163          b3 =  Dirichlet_boundary(conserved_quantities = num.array([2.0]))
     161         b1 =  Dirichlet_boundary(dirichlet_values = num.array([0.0]))
     162         b2 =  Dirichlet_boundary(dirichlet_values = num.array([1.0]))
     163         b3 =  Dirichlet_boundary(dirichlet_values = num.array([2.0]))
    164164         tags["1"] = b1
    165165         tags["2"] = b2

  • trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

    r8164 r8486  
    24442444if __name__ == "__main__":
    24452445    suite = unittest.makeSuite(Test_Quantity, 'test')
    2446     runner = unittest.TextTestRunner()
     2446    runner = unittest.TextTestRunner(verbosity=2)
    24472447    runner.run(suite)

  • trunk/anuga_core/source/anuga/operators/run_polygon_erosion.py

    r8485 r8486  
    103103
    104104#===============================================================================
    105 # Now to the standard model setup
     105# Now to the standard anuga model setup
    106106#===============================================================================
    107107
    … …  
    152152# Evolve system through time
    153153#-------------------------------------------------------------------------------
    154 for t in domain.evolve(yieldstep=0.2, finaltime=40.0):
     154for t in domain.evolve(yieldstep=0.2, finaltime=60.0):
    155155    domain.print_timestepping_statistics()
    156156    domain.print_operator_timestepping_statistics()

  • trunk/anuga_core/source/anuga/shallow_water/boundaries.py

    r8456 r8486  
    8585
    8686        return q
     87
     88
     89    def evaluate_segment_new(self, domain, segment_edges):
     90
     91        if segment_edges is None:
     92            return
     93        if domain is None:
     94            return
     95
     96        ids = segment_edges
     97
     98    def evaluate_segment(self, domain, segment_edges):
     99
     100        if segment_edges is None:
     101            return
     102        if domain is None:
     103            return
     104
     105
     106        ids = segment_edges
     107        vol_ids  = domain.boundary_cells[ids]
     108        edge_ids = domain.boundary_edges[ids]
     109
     110
     111        for j, name in enumerate(domain.evolved_quantities):
     112            Q = domain.quantities[name]
     113            Q.boundary_values[ids] = Q.edge_values[vol_ids,edge_ids]
     114
     115        #-----------------------------------------
     116        # Rotate Momentum and  velocity
     117        #-----------------------------------------
    87118
    88119

  • trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r8473 r8486  
    298298
    299299
     300
     301    def set_quantity(self, name, *args, **kwargs):
     302        """Set values for named quantity
     303
     304        We have to do something special for 'elevation'
     305        otherwise pass through to generic set_quantity
     306        """
     307
     308#        if name == 'elevation':
     309#            stage_c = self.get_quantity('stage').centroid_values
     310#            elev_c =  self.get_quantity('elevation').centroid_values
     311#            height_c = stage_c - elev_c
     312#            Generic_Domain.set_quantity(self, name, *args, **kwargs)
     313#            stage_c[:] = elev_c + height_c
     314#        else:
     315#            Generic_Domain.set_quantity(self, name, *args, **kwargs)
     316
     317        Generic_Domain.set_quantity(self, name, *args, **kwargs)
     318
     319
    300320    def set_store(self, flag=True):
    301321        """Set whether data saved to sww file.
    … …  
    10491069            if self.store is True:
    10501070                self.store_timestep()
     1071
     1072
     1073
    10511074
    10521075            # Pass control on to outer loop for more specific actions

  • trunk/anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r8476 r8486  
    28982898
    28992899
     2900
    29002901        if (number_of_boundaries[k] <= 1) {
    29012902            //==============================================
    … …  
    32813282            //ymom_vertex_values[k3 + 2] = ymom_centroid_values[k] + dqv[2];
    32823283        } // else [number_of_boundaries==2]
     3284
     3285
     3286
     3287
    32833288    } // for k=0 to number_of_elements-1
    32843289
    … …  
    33073312        }
    33083313    }
    3309    
     3314
     3315
    33103316    free(xmom_centroid_store);
    33113317    free(ymom_centroid_store);
    … …  
    41974203
    41984204
     4205
     4206
     4207
    41994208PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
    42004209    //
    … …  
    42564265    return PyArray_Return(r);
    42574266}
     4267
     4268
     4269PyObject *reflect_momentum_velocity(PyObject *self, PyObject *args) {
     4270  /*Rotate momentum and velcity for all edges in a segment
     4271   *
     4272    Python call:
     4273    reflect_momentum_velocity(domain, segment_ids)
     4274
     4275    Post conditions:
     4276        The edges of each triangle have values from a
     4277        limited linear reconstruction
     4278        based on centroid values
     4279
     4280  */
     4281
     4282    struct domain D;
     4283    PyObject *domain;
     4284    PyArrayObject *segment_ids;
     4285
     4286    long *seg_ids;
     4287
     4288    int err;
     4289
     4290    if (!PyArg_ParseTuple(args, "OO", &domain, &segment_ids)) {
     4291        report_python_error(AT, "could not parse input arguments");
     4292        return NULL;
     4293    }
     4294
     4295    // populate the C domain structure with pointers
     4296    // to the python domain data
     4297    get_python_domain(&D, domain);
     4298
     4299
     4300
     4301
     4302
     4303
     4304    //print_domain_struct(&D);
     4305
     4306    // Call underlying flux computation routine and update
     4307    // the explicit update arrays
     4308    err = _extrapolate_second_order_sw(&D);
     4309
     4310    if (err == -1) {
     4311        // Use error string set inside computational routine
     4312        return NULL;
     4313    }
     4314
     4315  return Py_BuildValue("");
     4316
     4317}// extrapolate_second-order_sw
     4318
     4319
     4320
    42584321
    42594322

  • trunk/anuga_core/source/anuga/shallow_water/test_data_manager.py

    r8411 r8486  
    784784        Br = Reflective_boundary(domain_time)
    785785        Bw = Time_boundary(domain=domain_time,
    786                          f=lambda t: [num.sin(t)+tide,3.*(20.+num.sin(t)+tide),2.*(20.+num.sin(t)+tide)])
     786                         function=lambda t: [num.sin(t)+tide,3.*(20.+num.sin(t)+tide),2.*(20.+num.sin(t)+tide)])
    787787        domain_time.set_boundary({'ocean': Bw,'otherocean': Br})
    788788       

  • trunk/anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r8472 r8486  
    730730        #--------------------------------------------------------------
    731731        # Time dependent boundary
    732         Bt = Time_boundary(domain=domain, f=lambda t: [t, 0.0, 0.0])
     732        Bt = Time_boundary(domain=domain, function=lambda t: [t, 0.0, 0.0])
    733733
    734734        Br = Reflective_boundary(domain)              # Reflective wall
    … …  
    66246624
    66256625         tags = {}
    6626          b1 =  Dirichlet_boundary(conserved_quantities = num.array([0.0]))
    6627          b2 =  Dirichlet_boundary(conserved_quantities = num.array([1.0]))
    6628          b3 =  Dirichlet_boundary(conserved_quantities = num.array([2.0]))
     6626         b1 =  Dirichlet_boundary(dirichlet_values = num.array([0.0]))
     6627         b2 =  Dirichlet_boundary(dirichlet_values = num.array([1.0]))
     6628         b3 =  Dirichlet_boundary(dirichlet_values = num.array([2.0]))
    66296629         tags["1"] = b1
    66306630         tags["2"] = b2

  • trunk/anuga_core/source/anuga/shallow_water/test_system.py

    r8068 r8486  
    5757        Bd = anuga.Dirichlet_boundary([tide,0.,0.]) # Constant boundary values
    5858        Bd = anuga.Time_boundary(domain=domain,     # Time dependent boundary 
    59                    f=lambda t: [t, 0.0, 0.0])
     59                   function=lambda t: [t, 0.0, 0.0])
    6060        domain.set_boundary({'exterior': Bd})
    6161        for t in domain.evolve(yieldstep = 10, finaltime = 20.0):
Note: See TracChangeset for help on using the changeset viewer.