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

Last change on this file since 4026 was 4026, checked in by jack, 17 years ago

Moved the vpython visualiser to obsolete_code and cleared out things that call it from other code.

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