source: inundation/pyvolution/advection.py @ 3354

Last change on this file since 3354 was 3167, checked in by linda, 19 years ago

Fixed the if self.tri_full_flag[k] == 1 in advection.py

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