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

Last change on this file since 4967 was 4967, checked in by steve, 16 years ago

Playing with f2py for adding C extensions

File size: 8.9 KB
Line 
1"""Class Domain -
22D triangular domains for finite-volume computations of
3the advection equation.
4
5This module contains a specialisation of class Domain from module domain.py
6consisting of methods specific to the advection equantion
7
8The equation is
9
10  u_t + (v_1 u)_x + (v_2 u)_y = 0
11
12There is only one conserved quantity, the stage u
13
14The advection equation is a very simple specialisation of the generic
15domain and may serve as an instructive example or a test of other
16components such as visualisation.
17
18Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
19Geoscience Australia, 2004
20"""
21
22
23#import logging, logging.config
24#logger = logging.getLogger('advection')
25#logger.setLevel(logging.WARNING)
26#
27#try:
28#    logging.config.fileConfig('log.ini')
29#except:
30#    pass
31
32
33from anuga.abstract_2d_finite_volumes.domain import *
34Generic_domain = Domain # Rename
35
36class Domain(Generic_domain):
37
38    def __init__(self,
39                 coordinates,
40                 vertices,
41                 boundary = None,
42                 tagged_elements = None,
43                 geo_reference = None,
44                 use_inscribed_circle=False,
45                 velocity = None,
46                 full_send_dict=None,
47                 ghost_recv_dict=None,
48                 processor=0,
49                 numproc=1
50                 ):
51
52        conserved_quantities = ['stage']
53        other_quantities = []
54        Generic_domain.__init__(self,
55                                source=coordinates,
56                                triangles=vertices,
57                                boundary=boundary,
58                                conserved_quantities=conserved_quantities,
59                                other_quantities=other_quantities,
60                                tagged_elements=tagged_elements,
61                                geo_reference=geo_reference,
62                                use_inscribed_circle=use_inscribed_circle,
63                                full_send_dict=full_send_dict,
64                                ghost_recv_dict=ghost_recv_dict,
65                                processor=processor,
66                                numproc=numproc)
67
68
69        if velocity is None:
70            self.velocity = [1,0]
71        else:
72            self.velocity = velocity
73
74        #Only first is implemented for advection
75        self.default_order = self.order = 1
76
77        self.smooth = True
78
79    def check_integrity(self):
80        Generic_domain.check_integrity(self)
81
82        msg = 'Conserved quantity must be "stage"'
83        assert self.conserved_quantities[0] == 'stage', msg
84
85
86    def flux_function(self, normal, ql, qr, zl=None, zr=None):
87        """Compute outward flux as inner product between velocity
88        vector v=(v_1, v_2) and normal vector n.
89
90        if <n,v> > 0 flux direction is outward bound and its magnitude is
91        determined by the quantity inside volume: ql.
92        Otherwise it is inbound and magnitude is determined by the
93        quantity outside the volume: qr.
94        """
95
96        v1 = self.velocity[0]
97        v2 = self.velocity[1]
98
99
100        normal_velocity = v1*normal[0] + v2*normal[1]
101
102        if normal_velocity < 0:
103            flux = qr * normal_velocity
104        else:
105            flux = ql * normal_velocity
106
107        max_speed = abs(normal_velocity)
108        return flux, max_speed
109
110    def compute_fluxes(self):
111
112       
113        self.compute_fluxes_ext()
114##         try:
115##             self.compute_fluxes_ext()
116##         except:
117##             self.compute_fluxes_python()
118 
119
120
121    def compute_fluxes_python(self):
122        """Compute all fluxes and the timestep suitable for all volumes
123        in domain.
124
125        Compute total flux for each conserved quantity using "flux_function"
126
127        Fluxes across each edge are scaled by edgelengths and summed up
128        Resulting flux is then scaled by area and stored in
129        domain.explicit_update
130
131        The maximal allowable speed computed by the flux_function
132        for each volume
133        is converted to a timestep that must not be exceeded. The minimum of
134        those is computed as the next overall timestep.
135
136        Post conditions:
137        domain.explicit_update is reset to computed flux values
138        domain.timestep is set to the largest step satisfying all volumes.
139        """
140
141        import sys
142        from Numeric import zeros, Float
143        from anuga.config import max_timestep
144
145        N = len(self)
146
147        neighbours = self.neighbours
148        neighbour_edges = self.neighbour_edges
149        normals = self.normals
150
151        areas = self.areas
152        radii = self.radii
153        edgelengths = self.edgelengths
154
155        timestep = max_timestep #FIXME: Get rid of this
156
157        #Shortcuts
158        Stage = self.quantities['stage']
159
160        #Arrays
161        stage = Stage.edge_values
162
163        stage_bdry = Stage.boundary_values
164
165        flux = zeros(1, Float) #Work array for summing up fluxes
166
167        #Loop
168        for k in range(N):
169            optimal_timestep = float(sys.maxint)
170
171            flux[:] = 0.  #Reset work array
172            for i in range(3):
173                #Quantities inside volume facing neighbour i
174                ql = stage[k, i]
175
176                #Quantities at neighbour on nearest face
177                n = neighbours[k,i]
178                if n < 0:
179                    m = -n-1 #Convert neg flag to index
180                    qr = stage_bdry[m]
181                else:
182                    m = neighbour_edges[k,i]
183                    qr = stage[n, m]
184
185
186                #Outward pointing normal vector
187                normal = normals[k, 2*i:2*i+2]
188
189                #Flux computation using provided function
190                edgeflux, max_speed = self.flux_function(normal, ql, qr)
191                flux -= edgeflux * edgelengths[k,i]
192
193                #Update optimal_timestep
194                if  self.tri_full_flag[k] == 1 :
195                    try:
196                        optimal_timestep = min(optimal_timestep, radii[k]/max_speed)
197                    except ZeroDivisionError:
198                        pass
199
200            #Normalise by area and store for when all conserved
201            #quantities get updated
202            flux /= areas[k]
203            Stage.explicit_update[k] = flux[0]
204
205            timestep = min(timestep, optimal_timestep)
206
207        self.timestep = timestep
208
209
210    def compute_fluxes_ext(self):
211        """Compute all fluxes and the timestep suitable for all volumes
212        in domain.
213
214        Compute total flux for each conserved quantity using "flux_function"
215
216        Fluxes across each edge are scaled by edgelengths and summed up
217        Resulting flux is then scaled by area and stored in
218        domain.explicit_update
219
220        The maximal allowable speed computed by the flux_function
221        for each volume
222        is converted to a timestep that must not be exceeded. The minimum of
223        those is computed as the next overall timestep.
224
225        Post conditions:
226        domain.explicit_update is reset to computed flux values
227        domain.timestep is set to the largest step satisfying all volumes.
228        """
229
230        import sys
231        from Numeric import zeros, Float
232        from anuga.config import max_timestep
233
234        N = len(self)
235
236        timestep = max_timestep
237
238        #Shortcuts
239        Stage = self.quantities['stage']
240
241        #Arrays
242        neighbours      = self.neighbours
243        neighbour_edges = self.neighbour_edges
244        normals         = self.normals
245        areas           = self.areas
246        radii           = self.radii
247        edgelengths     = self.edgelengths
248        tri_full_flag   = self.tri_full_flag
249
250        stage_edge      = Stage.edge_values
251        stage_bdry      = Stage.boundary_values
252        stage_update    = Stage.explicit_update
253
254        huge_timestep = float(sys.maxint)
255
256        v = self.velocity
257
258        from advection_ext import compute_fluxes
259
260
261        print 'N = ',N
262        print 'timestep = ',timestep
263        print 'huge_timestep = ',huge_timestep
264               
265        timestep = compute_fluxes(stage_edge,stage_bdry,stage_update,
266                                  neighbours,neighbour_edges,normals,
267                                  areas,radii,edgelengths,
268                                  tri_full_flag,
269                                  huge_timestep,max_timestep,v,N)
270       
271        print 'timestep out2 =',timestep
272
273        self.timestep = timestep
274
275
276    def evolve(self,
277               yieldstep = None,
278               finaltime = None,
279               duration = None,
280               skip_initial_step = False):
281
282        """Specialisation of basic evolve method from parent class
283        """
284
285        #Call basic machinery from parent class
286        for t in Generic_domain.evolve(self,
287                                       yieldstep=yieldstep,
288                                       finaltime=finaltime,
289                                       duration=duration,
290                                       skip_initial_step=skip_initial_step):
291
292            #Pass control on to outer loop for more specific actions
293            yield(t)
Note: See TracBrowser for help on using the repository browser.