source: inundation/pyvolution/advection.py @ 2813

Last change on this file since 2813 was 2813, checked in by steve, 18 years ago

Moving ghosts into domain.py

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