source: anuga_core/source/anuga/advection/advection.py @ 5242

Last change on this file since 5242 was 5242, checked in by steve, 16 years ago
File size: 8.5 KB
RevLine 
[4967]1"""Class Domain -
22D triangular domains for finite-volume computations of
3the advection equation.
4
5This module contains a specialisation of class Domain from module domain.py
6consisting of methods specific to the advection equantion
7
8The equation is
9
10  u_t + (v_1 u)_x + (v_2 u)_y = 0
11
12There is only one conserved quantity, the stage u
13
14The advection equation is a very simple specialisation of the generic
15domain and may serve as an instructive example or a test of other
16components such as visualisation.
17
18Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
19Geoscience Australia, 2004
20"""
21
22
23#import logging, logging.config
24#logger = logging.getLogger('advection')
25#logger.setLevel(logging.WARNING)
26#
27#try:
28#    logging.config.fileConfig('log.ini')
29#except:
30#    pass
31
32
33from anuga.abstract_2d_finite_volumes.domain import *
34Generic_domain = Domain # Rename
35
36class Domain(Generic_domain):
37
38    def __init__(self,
39                 coordinates,
40                 vertices,
41                 boundary = None,
42                 tagged_elements = None,
43                 geo_reference = None,
44                 use_inscribed_circle=False,
45                 velocity = None,
46                 full_send_dict=None,
47                 ghost_recv_dict=None,
48                 processor=0,
49                 numproc=1
50                 ):
51
52        conserved_quantities = ['stage']
53        other_quantities = []
54        Generic_domain.__init__(self,
55                                source=coordinates,
56                                triangles=vertices,
57                                boundary=boundary,
58                                conserved_quantities=conserved_quantities,
59                                other_quantities=other_quantities,
60                                tagged_elements=tagged_elements,
61                                geo_reference=geo_reference,
62                                use_inscribed_circle=use_inscribed_circle,
63                                full_send_dict=full_send_dict,
64                                ghost_recv_dict=ghost_recv_dict,
65                                processor=processor,
66                                numproc=numproc)
67
[4978]68        import Numeric
[4967]69        if velocity is None:
[4978]70            self.velocity = Numeric.array([1,0],'d')
[4967]71        else:
[4978]72            self.velocity = Numeric.array(velocity,'d')
[4967]73
74        #Only first is implemented for advection
75        self.default_order = self.order = 1
76
77        self.smooth = True
78
79    def check_integrity(self):
80        Generic_domain.check_integrity(self)
81
82        msg = 'Conserved quantity must be "stage"'
83        assert self.conserved_quantities[0] == 'stage', msg
84
85
86    def flux_function(self, normal, ql, qr, zl=None, zr=None):
87        """Compute outward flux as inner product between velocity
88        vector v=(v_1, v_2) and normal vector n.
89
90        if <n,v> > 0 flux direction is outward bound and its magnitude is
91        determined by the quantity inside volume: ql.
92        Otherwise it is inbound and magnitude is determined by the
93        quantity outside the volume: qr.
94        """
95
96        v1 = self.velocity[0]
97        v2 = self.velocity[1]
98
99
100        normal_velocity = v1*normal[0] + v2*normal[1]
101
102        if normal_velocity < 0:
103            flux = qr * normal_velocity
104        else:
105            flux = ql * normal_velocity
106
107        max_speed = abs(normal_velocity)
108        return flux, max_speed
109
[4978]110
111
[4967]112    def compute_fluxes(self):
[4978]113        """Compute all fluxes and the timestep suitable for all volumes
114        in domain.
[4967]115
[4978]116        Compute total flux for each conserved quantity using "flux_function"
[4967]117
[4978]118        Fluxes across each edge are scaled by edgelengths and summed up
119        Resulting flux is then scaled by area and stored in
120        domain.explicit_update
[4967]121
[4978]122        The maximal allowable speed computed by the flux_function
123        for each volume
124        is converted to a timestep that must not be exceeded. The minimum of
125        those is computed as the next overall timestep.
126
127        Post conditions:
128        domain.explicit_update is reset to computed flux values
129        domain.timestep is set to the largest step satisfying all volumes.
130        """
131
132        import sys
133        from Numeric import zeros, Float
134        from anuga.config import max_timestep
135
136
137        huge_timestep = float(sys.maxint)
138        Stage = self.quantities['stage']
139
140        """
141        print "======================================"
142        print "BEFORE compute_fluxes"
143        print "stage_update",Stage.explicit_update
144        print "stage_edge",Stage.edge_values
145        print "stage_bdry",Stage.boundary_values
146        print "neighbours",self.neighbours
147        print "neighbour_edges",self.neighbour_edges
148        print "normals",self.normals
149        print "areas",self.areas
150        print "radii",self.radii
151        print "edgelengths",self.edgelengths
152        print "tri_full_flag",self.tri_full_flag
153        print "huge_timestep",huge_timestep
154        print "max_timestep",max_timestep
155        print "velocity",self.velocity
156        """
157
158        import advection_ext           
[5242]159        self.flux_timestep = advection_ext.compute_fluxes(self, Stage, huge_timestep, max_timestep)
[4978]160
161
162
163    def evolve(self,
164               yieldstep = None,
165               finaltime = None,
166               duration = None,
167               skip_initial_step = False):
168
169        """Specialisation of basic evolve method from parent class
170        """
171
172        #Call basic machinery from parent class
173        for t in Generic_domain.evolve(self,
174                                       yieldstep=yieldstep,
175                                       finaltime=finaltime,
176                                       duration=duration,
177                                       skip_initial_step=skip_initial_step):
178
179            #Pass control on to outer loop for more specific actions
180            yield(t)
181
182
183
184
[4967]185    def compute_fluxes_python(self):
186        """Compute all fluxes and the timestep suitable for all volumes
187        in domain.
188
189        Compute total flux for each conserved quantity using "flux_function"
190
191        Fluxes across each edge are scaled by edgelengths and summed up
192        Resulting flux is then scaled by area and stored in
193        domain.explicit_update
194
195        The maximal allowable speed computed by the flux_function
196        for each volume
197        is converted to a timestep that must not be exceeded. The minimum of
198        those is computed as the next overall timestep.
199
200        Post conditions:
201        domain.explicit_update is reset to computed flux values
202        domain.timestep is set to the largest step satisfying all volumes.
203        """
204
205        import sys
206        from Numeric import zeros, Float
207        from anuga.config import max_timestep
208
209        N = len(self)
210
211        neighbours = self.neighbours
212        neighbour_edges = self.neighbour_edges
213        normals = self.normals
214
215        areas = self.areas
216        radii = self.radii
217        edgelengths = self.edgelengths
218
219        timestep = max_timestep #FIXME: Get rid of this
220
221        #Shortcuts
222        Stage = self.quantities['stage']
223
224        #Arrays
225        stage = Stage.edge_values
226
227        stage_bdry = Stage.boundary_values
228
229        flux = zeros(1, Float) #Work array for summing up fluxes
230
231        #Loop
232        for k in range(N):
233            optimal_timestep = float(sys.maxint)
234
235            flux[:] = 0.  #Reset work array
236            for i in range(3):
237                #Quantities inside volume facing neighbour i
238                ql = stage[k, i]
239
240                #Quantities at neighbour on nearest face
241                n = neighbours[k,i]
242                if n < 0:
243                    m = -n-1 #Convert neg flag to index
244                    qr = stage_bdry[m]
245                else:
246                    m = neighbour_edges[k,i]
247                    qr = stage[n, m]
248
249
250                #Outward pointing normal vector
251                normal = normals[k, 2*i:2*i+2]
252
253                #Flux computation using provided function
254                edgeflux, max_speed = self.flux_function(normal, ql, qr)
255                flux -= edgeflux * edgelengths[k,i]
256
257                #Update optimal_timestep
258                if  self.tri_full_flag[k] == 1 :
259                    try:
260                        optimal_timestep = min(optimal_timestep, radii[k]/max_speed)
261                    except ZeroDivisionError:
262                        pass
263
264            #Normalise by area and store for when all conserved
265            #quantities get updated
266            flux /= areas[k]
267            Stage.explicit_update[k] = flux[0]
268
269            timestep = min(timestep, optimal_timestep)
270
271        self.timestep = timestep
272
Note: See TracBrowser for help on using the repository browser.