source: inundation/pyvolution/advection.py @ 2131

Last change on this file since 2131 was 2050, checked in by steve, 19 years ago

Fixed bug in visualisation of advection domain

File size: 11.2 KB
Line 
1import sys
2from os import sep
3sys.path.append('..'+sep+'pyvolution')
4
5"""Class Domain -
62D triangular domains for finite-volume computations of
7the advection equation.
8
9This module contains a specialisation of class Domain from module domain.py
10consisting of methods specific to the advection equantion
11
12The equation is
13
14  u_t + (v_1 u)_x + (v_2 u)_y = 0
15
16There is only one conserved quantity, the stage u
17
18The advection equation is a very simple specialisation of the generic
19domain and may serve as an instructive example or a test of other
20components such as visualisation.
21
22Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
23Geoscience Australia, 2004
24"""
25
26
27import logging, logging.config
28logger = logging.getLogger('advection')
29logger.setLevel(logging.WARNING)
30
31try:
32    logging.config.fileConfig('log.ini')
33except:
34    pass
35
36
37from domain import *
38Generic_domain = Domain #Rename
39
40class Domain(Generic_domain):
41
42    def __init__(self, coordinates, vertices, boundary = None,
43                 tagged_elements = None, geo_reference = None,
44                 use_inscribed_circle=False, velocity = None):
45
46        conserved_quantities = ['stage']
47        other_quantities = []
48        Generic_domain.__init__(self, coordinates, vertices, boundary,
49                                conserved_quantities, other_quantities,
50                                tagged_elements, geo_reference,
51                                use_inscribed_circle)
52
53
54        if velocity is None:
55            self.velocity = [1,0]
56        else:
57            self.velocity = velocity
58
59        #Only first is implemented for advection
60        self.default_order = self.order = 1
61
62        #Realtime visualisation
63        self.visualiser = None
64        self.visualise  = False
65        self.visualise_color_stage = False
66        self.visualise_stage_range = 1.0
67        self.visualise_timer = True
68        self.visualise_range_z = None
69
70
71        self.smooth = True
72
73    def initialise_visualiser(self,scale_z=1.0,rect=None):
74        #Realtime visualisation
75        if self.visualiser is None:
76            from realtime_visualisation_new import Visualiser
77            self.visualiser = Visualiser(self,scale_z,rect)
78        self.visualise = True
79
80
81    def check_integrity(self):
82        Generic_domain.check_integrity(self)
83
84        msg = 'Conserved quantity must be "stage"'
85        assert self.conserved_quantities[0] == 'stage', msg
86
87
88    def flux_function(self, normal, ql, qr, zl=None, zr=None):
89        """Compute outward flux as inner product between velocity
90        vector v=(v_1, v_2) and normal vector n.
91
92        if <n,v> > 0 flux direction is outward bound and its magnitude is
93        determined by the quantity inside volume: ql.
94        Otherwise it is inbound and magnitude is determined by the
95        quantity outside the volume: qr.
96        """
97
98        v1 = self.velocity[0]
99        v2 = self.velocity[1]
100
101
102        normal_velocity = v1*normal[0] + v2*normal[1]
103
104        if normal_velocity < 0:
105            flux = qr * normal_velocity
106        else:
107            flux = ql * normal_velocity
108
109        max_speed = abs(normal_velocity)
110        return flux, max_speed
111
112    def compute_fluxes(self):
113
114        try:
115            import weave
116            self.weave_available = True
117        except:
118            self.weave_available = False
119
120        if self.weave_available:
121            self.compute_fluxes_weave()
122        else:
123            self.compute_fluxes_python()
124
125
126
127    def compute_fluxes_python(self):
128        """Compute all fluxes and the timestep suitable for all volumes
129        in domain.
130
131        Compute total flux for each conserved quantity using "flux_function"
132
133        Fluxes across each edge are scaled by edgelengths and summed up
134        Resulting flux is then scaled by area and stored in
135        domain.explicit_update
136
137        The maximal allowable speed computed by the flux_function
138        for each volume
139        is converted to a timestep that must not be exceeded. The minimum of
140        those is computed as the next overall timestep.
141
142        Post conditions:
143        domain.explicit_update is reset to computed flux values
144        domain.timestep is set to the largest step satisfying all volumes.
145        """
146
147        import sys
148        from Numeric import zeros, Float
149        from config import max_timestep
150
151        N = self.number_of_elements
152
153        neighbours = self.neighbours
154        neighbour_edges = self.neighbour_edges
155        normals = self.normals
156
157        areas = self.areas
158        radii = self.radii
159        edgelengths = self.edgelengths
160
161        timestep = max_timestep #FIXME: Get rid of this
162
163        #Shortcuts
164        Stage = self.quantities['stage']
165
166        #Arrays
167        stage = Stage.edge_values
168
169        stage_bdry = Stage.boundary_values
170
171        flux = zeros(1, Float) #Work array for summing up fluxes
172
173        #Loop
174        for k in range(N):
175            optimal_timestep = float(sys.maxint)
176
177            flux[:] = 0.  #Reset work array
178            for i in range(3):
179                #Quantities inside volume facing neighbour i
180                ql = stage[k, i]
181
182                #Quantities at neighbour on nearest face
183                n = neighbours[k,i]
184                if n < 0:
185                    m = -n-1 #Convert neg flag to index
186                    qr = stage_bdry[m]
187                else:
188                    m = neighbour_edges[k,i]
189                    qr = stage[n, m]
190
191
192                #Outward pointing normal vector
193                normal = normals[k, 2*i:2*i+2]
194
195                #Flux computation using provided function
196                edgeflux, max_speed = self.flux_function(normal, ql, qr)
197                flux -= edgeflux * edgelengths[k,i]
198
199                #Update optimal_timestep
200                try:
201                    optimal_timestep = min(optimal_timestep, radii[k]/max_speed)
202                except ZeroDivisionError:
203                    pass
204
205            #Normalise by area and store for when all conserved
206            #quantities get updated
207            flux /= areas[k]
208            Stage.explicit_update[k] = flux[0]
209
210            timestep = min(timestep, optimal_timestep)
211
212        self.timestep = timestep
213
214    def compute_fluxes_weave(self):
215        """Compute all fluxes and the timestep suitable for all volumes
216        in domain.
217
218        Compute total flux for each conserved quantity using "flux_function"
219
220        Fluxes across each edge are scaled by edgelengths and summed up
221        Resulting flux is then scaled by area and stored in
222        domain.explicit_update
223
224        The maximal allowable speed computed by the flux_function
225        for each volume
226        is converted to a timestep that must not be exceeded. The minimum of
227        those is computed as the next overall timestep.
228
229        Post conditions:
230        domain.explicit_update is reset to computed flux values
231        domain.timestep is set to the largest step satisfying all volumes.
232        """
233
234        import sys
235        from Numeric import zeros, Float
236        from config import max_timestep
237
238        import weave
239        from weave import converters
240
241        N = self.number_of_elements
242
243
244        timestep    = zeros( 1, Float);
245        timestep[0] = float(max_timestep) #FIXME: Get rid of this
246
247        #Shortcuts
248        Stage = self.quantities['stage']
249
250        #Arrays
251        neighbours      = self.neighbours
252        neighbour_edges = self.neighbour_edges
253        normals         = self.normals
254        areas           = self.areas
255        radii           = self.radii
256        edgelengths     = self.edgelengths
257
258        stage_edge      = Stage.edge_values
259        stage_bdry      = Stage.boundary_values
260        stage_update    = Stage.explicit_update
261
262        huge_timestep = float(sys.maxint)
263
264        v1 = self.velocity[0]
265        v2 = self.velocity[1]
266
267        code = """
268        //Loop
269
270        double qr,ql;
271        int m,n;
272        double normal[2];
273        double normal_velocity;
274        double flux, edgeflux;
275        double max_speed;
276        double optimal_timestep;
277        for (int k=0; k<N; k++){
278
279            optimal_timestep = huge_timestep;
280            flux = 0.0;  //Reset work array
281            for (int i=0; i<3; i++){
282                //Quantities inside volume facing neighbour i
283                ql = stage_edge(k, i);
284
285                //Quantities at neighbour on nearest face
286                n = neighbours(k,i);
287                if (n < 0) {
288                    m = -n-1; //Convert neg flag to index
289                    qr = stage_bdry(m);
290                } else {
291                    m = neighbour_edges(k,i);
292                    qr = stage_edge(n, m);
293                }
294
295
296                //Outward pointing normal vector
297                for (int j=0; j<2; j++){
298                    normal[j] = normals(k, 2*i+j);
299                }
300
301
302                //Flux computation using provided function
303                normal_velocity = v1*normal[0] + v2*normal[1];
304
305                if (normal_velocity < 0) {
306                    edgeflux = qr * normal_velocity;
307                } else {
308                    edgeflux = ql * normal_velocity;
309                }
310
311                max_speed = fabs(normal_velocity);
312                flux = flux - edgeflux * edgelengths(k,i);
313
314                //Update optimal_timestep
315                if (max_speed != 0.0) {
316                    optimal_timestep = (optimal_timestep>radii(k)/max_speed) ? radii(k)/max_speed : optimal_timestep;
317                }
318
319            }
320
321            //Normalise by area and store for when all conserved
322            //quantities get updated
323            stage_update(k) = flux/areas(k);
324
325            timestep(0) = (timestep(0)>optimal_timestep) ? optimal_timestep : timestep(0);
326
327        }
328        """
329
330        logger.debug('Trying to weave advection.compute_fluxes')
331        weave.inline(code, ['stage_edge','stage_bdry','stage_update',
332                             'neighbours','neighbour_edges','normals',
333                             'areas','radii','edgelengths','huge_timestep',
334                             'timestep','v1','v2','N'],
335                             type_converters = converters.blitz, compiler='gcc');
336
337        self.timestep = timestep[0]
338
339
340    def evolve(self, yieldstep = None, finaltime = None):
341        """Specialisation of basic evolve method from parent class
342        """
343
344        #Initialise real time viz if requested
345        if self.visualise is True and self.time == 0.0:
346            #pass
347            #import realtime_visualisation_new as visualise
348            #visualise.create_surface(self)
349            self.initialise_visualiser()
350            self.visualiser.setup_all()
351            self.visualiser.update_timer()
352
353
354        #Call basic machinery from parent class
355        for t in Generic_domain.evolve(self, yieldstep, finaltime):
356
357
358            #Real time viz
359            if self.visualise is True:
360                #pass
361                self.visualiser.update_all()
362                self.visualiser.update_timer()
363
364
365            #Pass control on to outer loop for more specific actions
366            yield(t)
Note: See TracBrowser for help on using the repository browser.