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

Last change on this file since 5847 was 5847, checked in by steve, 15 years ago

Changed parallel_api so that global mesh only needs to
be constructed on processor 0

File size: 9.2 KB
Line 
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
68        import Numeric
69        if velocity is None:
70            self.velocity = Numeric.array([1,0],'d')
71        else:
72            self.velocity = Numeric.array(velocity,'d')
73
74        #Only first is implemented for advection
75        self.set_default_order(1)
76        self.set_beta(1.0)
77       
78        self.smooth = True
79
80    def check_integrity(self):
81        Generic_domain.check_integrity(self)
82
83        msg = 'Conserved quantity must be "stage"'
84        assert self.conserved_quantities[0] == 'stage', msg
85
86
87    def distribute_to_vertices_and_edges(self):
88        """Extrapolate conserved quantities from centroid to
89        vertices and edge-midpoints for each volume
90
91        Default implementation is straight first order,
92        i.e. constant values throughout each element and
93        no reference to non-conserved quantities.
94        """
95
96        for name in self.conserved_quantities:
97            Q = self.quantities[name]
98            if self._order_ == 1:
99                Q.extrapolate_first_order()
100            elif self._order_ == 2:
101                Q.extrapolate_second_order_and_limit_by_edge()
102                #Q.limit()
103            else:
104                raise 'Unknown order'
105            #Q.interpolate_from_vertices_to_edges()
106
107
108
109
110    def flux_function(self, normal, ql, qr, zl=None, zr=None):
111        """Compute outward flux as inner product between velocity
112        vector v=(v_1, v_2) and normal vector n.
113
114        if <n,v> > 0 flux direction is outward bound and its magnitude is
115        determined by the quantity inside volume: ql.
116        Otherwise it is inbound and magnitude is determined by the
117        quantity outside the volume: qr.
118        """
119
120        v1 = self.velocity[0]
121        v2 = self.velocity[1]
122
123
124        normal_velocity = v1*normal[0] + v2*normal[1]
125
126        if normal_velocity < 0:
127            flux = qr * normal_velocity
128        else:
129            flux = ql * normal_velocity
130
131        max_speed = abs(normal_velocity)
132        return flux, max_speed
133
134
135
136    def compute_fluxes(self):
137        """Compute all fluxes and the timestep suitable for all volumes
138        in domain.
139
140        Compute total flux for each conserved quantity using "flux_function"
141
142        Fluxes across each edge are scaled by edgelengths and summed up
143        Resulting flux is then scaled by area and stored in
144        domain.explicit_update
145
146        The maximal allowable speed computed by the flux_function
147        for each volume
148        is converted to a timestep that must not be exceeded. The minimum of
149        those is computed as the next overall timestep.
150
151        Post conditions:
152        domain.explicit_update is reset to computed flux values
153        domain.timestep is set to the largest step satisfying all volumes.
154        """
155
156        import sys
157        from Numeric import zeros, Float
158        from anuga.config import max_timestep
159
160
161        huge_timestep = float(sys.maxint)
162        Stage = self.quantities['stage']
163
164        """
165        print "======================================"
166        print "BEFORE compute_fluxes"
167        print "stage_update",Stage.explicit_update
168        print "stage_edge",Stage.edge_values
169        print "stage_bdry",Stage.boundary_values
170        print "neighbours",self.neighbours
171        print "neighbour_edges",self.neighbour_edges
172        print "normals",self.normals
173        print "areas",self.areas
174        print "radii",self.radii
175        print "edgelengths",self.edgelengths
176        print "tri_full_flag",self.tri_full_flag
177        print "huge_timestep",huge_timestep
178        print "max_timestep",max_timestep
179        print "velocity",self.velocity
180        """
181
182        import advection_ext           
183        self.flux_timestep = advection_ext.compute_fluxes(self, Stage, huge_timestep, max_timestep)
184
185
186
187##     def evolve(self,
188##                yieldstep = None,
189##                finaltime = None,
190##                duration = None,
191##                skip_initial_step = False):
192
193##         """Specialisation of basic evolve method from parent class
194##         """
195
196##         #Call basic machinery from parent class
197##         for t in Generic_domain.evolve(self,
198##                                        yieldstep=yieldstep,
199##                                        finaltime=finaltime,
200##                                        duration=duration,
201##                                        skip_initial_step=skip_initial_step):
202
203##             #Pass control on to outer loop for more specific actions
204##             yield(t)
205
206
207
208
209    def compute_fluxes_python(self):
210        """Compute all fluxes and the timestep suitable for all volumes
211        in domain.
212
213        Compute total flux for each conserved quantity using "flux_function"
214
215        Fluxes across each edge are scaled by edgelengths and summed up
216        Resulting flux is then scaled by area and stored in
217        domain.explicit_update
218
219        The maximal allowable speed computed by the flux_function
220        for each volume
221        is converted to a timestep that must not be exceeded. The minimum of
222        those is computed as the next overall timestep.
223
224        Post conditions:
225        domain.explicit_update is reset to computed flux values
226        domain.timestep is set to the largest step satisfying all volumes.
227        """
228
229        import sys
230        from Numeric import zeros, Float
231        from anuga.config import max_timestep
232
233        N = len(self)
234
235        neighbours = self.neighbours
236        neighbour_edges = self.neighbour_edges
237        normals = self.normals
238
239        areas = self.areas
240        radii = self.radii
241        edgelengths = self.edgelengths
242
243        timestep = max_timestep #FIXME: Get rid of this
244
245        #Shortcuts
246        Stage = self.quantities['stage']
247
248        #Arrays
249        stage = Stage.edge_values
250
251        stage_bdry = Stage.boundary_values
252
253        flux = zeros(1, Float) #Work array for summing up fluxes
254
255        #Loop
256        for k in range(N):
257            optimal_timestep = float(sys.maxint)
258
259            flux[:] = 0.  #Reset work array
260            for i in range(3):
261                #Quantities inside volume facing neighbour i
262                ql = stage[k, i]
263
264                #Quantities at neighbour on nearest face
265                n = neighbours[k,i]
266                if n < 0:
267                    m = -n-1 #Convert neg flag to index
268                    qr = stage_bdry[m]
269                else:
270                    m = neighbour_edges[k,i]
271                    qr = stage[n, m]
272
273
274                #Outward pointing normal vector
275                normal = normals[k, 2*i:2*i+2]
276
277                #Flux computation using provided function
278                edgeflux, max_speed = self.flux_function(normal, ql, qr)
279                flux -= edgeflux * edgelengths[k,i]
280
281                #Update optimal_timestep
282                if  self.tri_full_flag[k] == 1 :
283                    try:
284                        optimal_timestep = min(optimal_timestep, radii[k]/max_speed)
285                    except ZeroDivisionError:
286                        pass
287
288            #Normalise by area and store for when all conserved
289            #quantities get updated
290            flux /= areas[k]
291            Stage.explicit_update[k] = flux[0]
292
293            timestep = min(timestep, optimal_timestep)
294
295        self.timestep = timestep
296
Note: See TracBrowser for help on using the repository browser.