source: inundation/pyvolution/advection.py @ 2267

Last change on this file since 2267 was 2152, checked in by linda, 19 years ago

Added caching to the sw_merimbula files. Also updated some of the comments

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