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

Last change on this file since 3567 was 3567, checked in by ole, 18 years ago

More heavy handed restructuring towards the agreed formats.
File size: 12.6 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 anuga.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                if  self.tri_full_flag[k] == 1 :
221                    try:
222                        optimal_timestep = min(optimal_timestep, radii[k]/max_speed)
223                    except ZeroDivisionError:
224                        pass
225
226            #Normalise by area and store for when all conserved
227            #quantities get updated
228            flux /= areas[k]
229            Stage.explicit_update[k] = flux[0]
230
231            timestep = min(timestep, optimal_timestep)
232
233        self.timestep = timestep
234
235    def compute_fluxes_weave(self):
236        """Compute all fluxes and the timestep suitable for all volumes
237        in domain.
238
239        Compute total flux for each conserved quantity using "flux_function"
240
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 anuga.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
278        tri_full_flag   = self.tri_full_flag
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
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                    }
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
354        logger.debug('Trying to weave advection.compute_fluxes')
355        weave.inline(code, ['stage_edge','stage_bdry','stage_update',
356                             'neighbours','neighbour_edges','normals',
357                             'areas','radii','edgelengths','tri_full_flag',
358                             'huge_timestep',
359                             'timestep','v1','v2','N'],
360                             type_converters = converters.blitz, compiler='gcc');
361
362        self.timestep = timestep[0]
363
364
365
366    def evolve(self,
367               yieldstep = None,
368               finaltime = None,
369               duration = None,
370               skip_initial_step = False):
371
372        """Specialisation of basic evolve method from parent class
373        """
374
375        #Initialise real time viz if requested
376        if self.visualise is True and self.time == 0.0:
377            #pass
378            #import realtime_visualisation_new as visualise
379            #visualise.create_surface(self)
380            self.initialise_visualiser()
381            self.visualiser.setup_all()
382            self.visualiser.update_timer()
383
384
385        #Call basic machinery from parent class
386        for t in Generic_domain.evolve(self,
387                                       yieldstep=yieldstep,
388                                       finaltime=finaltime,
389                                       duration=duration,
390                                       skip_initial_step=skip_initial_step):
391
392
393
394
395            #Real time viz
396            if self.visualise is True:
397                #pass
398                self.visualiser.update_all()
399                self.visualiser.update_timer()
400
401
402            #Pass control on to outer loop for more specific actions
403            yield(t)
Note: See TracBrowser for help on using the repository browser.