source: trunk/anuga_work/development/2010-projects/anuga_1d/sww/sww_domain.py @ 8251

Last change on this file since 8251 was 8251, checked in by steve, 12 years ago

Tracking down NaN in 1d anuga code

File size: 16.4 KB
Line 
1"""Class Domain -
21D interval domains for finite-volume computations of
3the shallow water wave equation.
4
5This module contains a specialisation of class Domain from module domain.py
6consisting of methods specific to the Shallow Water Wave Equation
7
8
9U_t + E_x = S
10
11where
12
13U = [w, uh]
14E = [uh, u^2h + gh^2/2]
15S represents source terms forcing the system
16(e.g. gravity, friction, wind stress, ...)
17
18and _t, _x, _y denote the derivative with respect to t, x and y respectiely.
19
20The quantities are
21
22symbol    variable name    explanation
23x         x                horizontal distance from origin [m]
24z         elevation        elevation of bed on which flow is modelled [m]
25h         height           water height above z [m]
26w         stage            absolute water level, w = z+h [m]
27u                          speed in the x direction [m/s]
28uh        xmomentum        momentum in the x direction [m^2/s]
29
30eta                        mannings friction coefficient [to appear]
31nu                         wind stress coefficient [to appear]
32
33The conserved quantities are w, uh
34
35For details see e.g.
36Christopher Zoppou and Stephen Roberts,
37Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
38Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999
39
40
41John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
42Geoscience Australia, 2006
43"""
44
45import numpy
46
47from anuga_1d.base.generic_domain import Generic_domain
48from sww_boundary_conditions import *
49from sww_forcing_terms import *
50
51#Shallow water domain
52class Domain(Generic_domain):
53
54    def __init__(self, coordinates, boundary = None, tagged_elements = None):
55
56        conserved_quantities = ['stage', 'xmomentum']
57        evolved_quantities   = ['stage', 'xmomentum', 'elevation', 'height', 'velocity']
58        other_quantities     = ['friction']
59        Generic_domain.__init__(self,
60                                coordinates          = coordinates,
61                                boundary             = boundary,
62                                conserved_quantities = conserved_quantities,
63                                evolved_quantities   = evolved_quantities,
64                                other_quantities     = other_quantities,
65                                tagged_elements      = tagged_elements)
66       
67        from anuga_1d.config import minimum_allowed_height, g, h0
68        self.minimum_allowed_height = minimum_allowed_height
69        self.g = g
70        self.h0 = h0
71
72        #forcing terms
73        self.forcing_terms.append(gravity)
74        #self.forcing_terms.append(manning_friction)
75       
76       
77        #Stored output
78        self.store = True
79        self.format = 'sww'
80        self.smooth = True
81
82        #Evolve parametrs
83        self.cfl = 1.0
84       
85        #Reduction operation for get_vertex_values
86        from anuga_1d.base.util import mean
87        self.reduction = mean
88        #self.reduction = min  #Looks better near steep slopes
89
90        self.quantities_to_be_stored = ['stage','xmomentum']
91
92        self.__doc__ = 'sww_domain'
93
94        self.set_spatial_order(2)
95        self.set_timestepping_method('rk2')
96        self.set_CFL(1.0)
97        self.set_limiter("minmod_kurganov")
98        self.set_beta(1.5)
99
100
101    def set_quantities_to_be_stored(self, q):
102        """Specify which quantities will be stored in the sww file.
103
104        q must be either:
105          - the name of a quantity
106          - a list of quantity names
107          - None
108
109        In the two first cases, the named quantities will be stored at each
110        yieldstep
111        (This is in addition to the quantities elevation and friction) 
112        If q is None, storage will be switched off altogether.
113        """
114
115
116        if q is None:
117            self.quantities_to_be_stored = []           
118            self.store = False
119            return
120
121        if isinstance(q, basestring):
122            q = [q] # Turn argument into a list
123
124        #Check correcness   
125        for quantity_name in q:
126            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
127            assert quantity_name in self.conserved_quantities, msg
128       
129        self.quantities_to_be_stored = q
130       
131
132
133    def check_integrity(self):
134        Generic_Domain.check_integrity(self)
135        #Check that we are solving the shallow water wave equation
136
137        msg = 'First conserved quantity must be "stage"'
138        assert self.conserved_quantities[0] == 'stage', msg
139        msg = 'Second conserved quantity must be "xmomentum"'
140        assert self.conserved_quantities[1] == 'xmomentum', msg
141
142
143
144    def compute_fluxes(self):
145        #Call correct module function
146        #(either from this module or C-extension)
147
148        import sys
149
150
151        timestep = float(sys.maxint)
152
153        stage    = self.quantities['stage']
154        xmom     = self.quantities['xmomentum']
155        bed      = self.quantities['elevation']
156        height   = self.quantities['height']
157        velocity = self.quantities['velocity']
158
159
160        #from anuga_1d.sww.sww_comp_flux_ext import compute_fluxes_ext_short as comp_flux_ext
161        #self.flux_timestep = comp_flux_ext(timestep,self,stage,xmom,bed)
162
163        from anuga_1d.sww.sww_vel_comp_flux_ext import compute_fluxes_vel_ext as comp_flux_ext
164        self.flux_timestep = comp_flux_ext(timestep,self,stage,xmom,bed,height,velocity)
165
166
167
168    def distribute_to_vertices_and_edges(self):
169        #Call correct module function
170        #(either from this module or C-extension)
171        distribute_to_vertices_and_edges(self)
172
173       
174    def evolve(self, yieldstep = None, finaltime = None, duration = None,
175               skip_initial_step = False):
176        """Specialisation of basic evolve method from parent class
177        """
178
179        #Call basic machinery from parent class
180        for t in Generic_domain.evolve(self, yieldstep, finaltime,duration,
181                                       skip_initial_step):
182
183            #Pass control on to outer loop for more specific actions
184            yield(t)
185
186
187
188
189
190#=============== End of Shallow Water Domain ===============================
191
192
193# Module functions for gradient limiting (distribute_to_vertices_and_edges)
194
195def distribute_to_vertices_and_edges(domain):
196    """Distribution from centroids to vertices specific to the
197    shallow water wave
198    equation.
199
200    It will ensure that h (w-z) is always non-negative even in the
201    presence of steep bed-slopes by taking a weighted average between shallow
202    and deep cases.
203
204    In addition, all conserved quantities get distributed as per either a
205    constant (order==1) or a piecewise linear function (order==2).
206
207    FIXME: more explanation about removal of artificial variability etc
208
209    Precondition:
210      All quantities defined at centroids and bed elevation defined at
211      vertices.
212
213    Postcondition
214      Conserved quantities defined at vertices
215
216    """
217
218    #from config import optimised_gradient_limiter
219
220    #Remove very thin layers of water
221    #protect_against_infinitesimal_and_negative_heights(domain) 
222
223    import sys
224    import numpy as np
225    from anuga_1d.config import epsilon, h0
226       
227    N = domain.number_of_elements
228
229    #Shortcuts
230    Stage      = domain.quantities['stage']
231    Xmom       = domain.quantities['xmomentum']
232    Bed        = domain.quantities['elevation']
233    Height     = domain.quantities['height']
234    Velocity   = domain.quantities['velocity']
235
236    #Arrays   
237    w_C   = Stage.centroid_values   
238    uh_C  = Xmom.centroid_values   
239    z_C   = Bed.centroid_values
240    h_C   = Height.centroid_values
241    u_C   = Velocity.centroid_values
242
243    w_V = domain.quantities['stage'].vertex_values
244    z_V   = domain.quantities['elevation'].vertex_values
245    h_V     = domain.quantities['height'].vertex_values
246    u_V     = domain.quantities['velocity'].vertex_values
247    uh_V  = domain.quantities['xmomentum'].vertex_values
248
249
250#    print 'w_C', np.any(np.isnan(w_C))
251#    print 'uh_C', np.any(np.isnan(uh_C))
252#    print 'z_C', np.any(np.isnan(z_C))
253#    print 'h_C', np.any(np.isnan(h_C))
254#    print 'u_C', np.any(np.isnan(u_C))
255#
256#    print 'w_V', np.any(np.isnan(w_V))
257#    print 'uh_V', np.any(np.isnan(uh_V))
258#    print 'z_V', np.any(np.isnan(z_V))
259#    print 'h_V', np.any(np.isnan(h_V))
260#    print 'u_V', np.any(np.isnan(u_V))
261       
262    #Calculate height (and fix negatives)better be non-negative!
263    w_C[:] = numpy.maximum(w_C, z_C)
264    h_C[:] = w_C - z_C
265
266    u_C[:] = numpy.where(h_C <= 1.0e-14, 0.0, uh_C/h_C)
267   
268
269    #print 'domain.order', domain.order
270   
271    for name in [ 'velocity', 'stage' ]:
272        Q = domain.quantities[name]
273        if domain.order == 1:
274            Q.extrapolate_first_order()
275        elif domain.order == 2:
276            #print "add extrapolate second order to shallow water"
277            #if name != 'height':
278            Q.extrapolate_second_order()
279            #Q.limit()
280        else:
281            raise 'Unknown order'
282
283
284#    print 'w_C', np.any(np.isnan(w_C))
285#    print 'uh_C', np.any(np.isnan(uh_C))
286#    print 'z_C', np.any(np.isnan(z_C))
287#    print 'h_C', np.any(np.isnan(h_C))
288#    print 'u_C', np.any(np.isnan(u_C))
289
290
291
292    h_V[:,:]    = w_V - z_V
293
294    #print 'any' ,numpy.any( h_V[:,0] < 0.0)
295    #print 'any' ,numpy.any( h_V[:,1] < 0.0)
296
297
298    h_0 = numpy.where(h_V[:,0] < 0.0, 0.0, h_V[:,0])
299    h_1 = numpy.where(h_V[:,0] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,1])
300
301    h_V[:,0] = h_0
302    h_V[:,1] = h_1
303
304
305    h_0 = numpy.where(h_V[:,1] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,0])
306    h_1 = numpy.where(h_V[:,1] < 0.0, 0.0, h_V[:,1])
307
308    h_V[:,0] = h_0
309    h_V[:,1] = h_1
310   
311    #print 'any' ,numpy.any( h_V[:,0] < 0.0)
312    #print 'any' ,numpy.any( h_V[:,1] < 0.0)
313
314    #h00 = 1e-12
315    #print h00
316
317    #h_V[:,:] = numpy.where (h_V <= h00, 0.0, h_V)
318    #u_V[:,:] = numpy.where (h_V <= h00, 0.0, u_V)
319
320#    # protect from edge values going negative
321#    h_V[:,1] = numpy.where(h_V[:,0] < 0.0 , h_V[:,1]-h_V[:,0], h_V[:,1])
322#    h_V[:,0] = numpy.where(h_V[:,0] < 0.0 , 0.0, h_V[:,0])
323#
324#    h_V[:,0] = numpy.where(h_V[:,1] < 0.0 , h_V[:,0]-h_V[:,1], h_V[:,0])
325#    h_V[:,1] = numpy.where(h_V[:,1] < 0.0 , 0.0, h_V[:,1])
326#
327#
328    w_V[:,:] = z_V + h_V
329    #bed_V[:,:] = stage_V - h_V
330    uh_V[:,:] = u_V * h_V
331
332
333#    print 'w_V', np.any(np.isnan(w_V))
334#    print 'uh_V', np.any(np.isnan(uh_V))
335#    print 'z_V', np.any(np.isnan(z_V))
336#    print 'h_V', np.any(np.isnan(h_V))
337#    print 'u_V', np.any(np.isnan(u_V))
338   
339    return
340#
341
342
343
344
345
346
347
348#
349def protect_against_infinitesimal_and_negative_heights(domain):
350    """Protect against infinitesimal heights and associated high velocities
351    """
352
353    #Shortcuts
354    wc    = domain.quantities['stage'].centroid_values
355    zc    = domain.quantities['elevation'].centroid_values
356    xmomc = domain.quantities['xmomentum'].centroid_values
357    hc = wc - zc  #Water depths at centroids
358
359    zv = domain.quantities['elevation'].vertex_values
360    wv = domain.quantities['stage'].vertex_values
361    hv = wv-zv
362    xmomv = domain.quantities['xmomentum'].vertex_values
363    #remove the above two lines and corresponding code below
364
365    #Update
366    #FIXME replace with numpy.where
367    for k in range(domain.number_of_elements):
368
369        if hc[k] < domain.minimum_allowed_height:
370            #Control stage
371            if hc[k] < domain.epsilon:
372                wc[k] = zc[k] # Contain 'lost mass' error
373                wv[k,0] = zv[k,0]
374                wv[k,1] = zv[k,1]
375               
376            xmomc[k] = 0.0
377           
378        #N = domain.number_of_elements
379        #if (k == 0) | (k==N-1):
380        #    wc[k] = zc[k] # Contain 'lost mass' error
381        #    wv[k,0] = zv[k,0]
382        #    wv[k,1] = zv[k,1]
383   
384def h_limiter(domain):
385    """Limit slopes for each volume to eliminate artificial variance
386    introduced by e.g. second order extrapolator
387
388    limit on h = w-z
389
390    This limiter depends on two quantities (w,z) so it resides within
391    this module rather than within quantity.py
392    """
393
394    N = domain.number_of_elements
395    beta_h = domain.beta_h
396
397    #Shortcuts
398    wc = domain.quantities['stage'].centroid_values
399    zc = domain.quantities['elevation'].centroid_values
400    hc = wc - zc
401
402    wv = domain.quantities['stage'].vertex_values
403    zv = domain.quantities['elevation'].vertex_values
404    hv = wv-zv
405
406    hvbar = zeros(hv.shape, numpy.float) #h-limited values
407
408    #Find min and max of this and neighbour's centroid values
409    hmax = zeros(hc.shape, numpy.float)
410    hmin = zeros(hc.shape, numpy.float)
411
412    for k in range(N):
413        hmax[k] = hmin[k] = hc[k]
414        #for i in range(3):
415        for i in range(2):   
416            n = domain.neighbours[k,i]
417            if n >= 0:
418                hn = hc[n] #Neighbour's centroid value
419
420                hmin[k] = min(hmin[k], hn)
421                hmax[k] = max(hmax[k], hn)
422
423
424    #Diffences between centroids and maxima/minima
425    dhmax = hmax - hc
426    dhmin = hmin - hc
427
428    #Deltas between vertex and centroid values
429    dh = zeros(hv.shape, numpy.float)
430    #for i in range(3):
431    for i in range(2):
432        dh[:,i] = hv[:,i] - hc
433
434    #Phi limiter
435    for k in range(N):
436
437        #Find the gradient limiter (phi) across vertices
438        phi = 1.0
439        #for i in range(3):
440        for i in range(2):
441            r = 1.0
442            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
443            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
444
445            phi = min( min(r*beta_h, 1), phi )
446
447        #Then update using phi limiter
448        #for i in range(3):
449        for i in range(2):
450            hvbar[k,i] = hc[k] + phi*dh[k,i]
451
452    return hvbar
453
454def balance_deep_and_shallow(domain):
455    """Compute linear combination between stage as computed by
456    gradient-limiters limiting using w, and stage computed by
457    gradient-limiters limiting using h (h-limiter).
458    The former takes precedence when heights are large compared to the
459    bed slope while the latter takes precedence when heights are
460    relatively small.  Anything in between is computed as a balanced
461    linear combination in order to avoid numpyal disturbances which
462    would otherwise appear as a result of hard switching between
463    modes.
464
465    The h-limiter is always applied irrespective of the order.
466    """
467
468    #Shortcuts
469    wc = domain.quantities['stage'].centroid_values
470    zc = domain.quantities['elevation'].centroid_values
471    hc = wc - zc
472
473    wv = domain.quantities['stage'].vertex_values
474    zv = domain.quantities['elevation'].vertex_values
475    hv = wv-zv
476
477    #Limit h
478    hvbar = h_limiter(domain)
479
480    for k in range(domain.number_of_elements):
481        #Compute maximal variation in bed elevation
482        #  This quantitiy is
483        #    dz = max_i abs(z_i - z_c)
484        #  and it is independent of dimension
485        #  In the 1d case zc = (z0+z1)/2
486        #  In the 2d case zc = (z0+z1+z2)/3
487
488        dz = max(abs(zv[k,0]-zc[k]),
489                 abs(zv[k,1]-zc[k]))#,
490#                 abs(zv[k,2]-zc[k]))
491
492
493        hmin = min( hv[k,:] )
494
495        #Create alpha in [0,1], where alpha==0 means using the h-limited
496        #stage and alpha==1 means using the w-limited stage as
497        #computed by the gradient limiter (both 1st or 2nd order)
498
499        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
500        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
501
502        if dz > 0.0:
503            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
504        else:
505            #Flat bed
506            alpha = 1.0
507
508        alpha  = 0.0
509        #Let
510        #
511        #  wvi be the w-limited stage (wvi = zvi + hvi)
512        #  wvi- be the h-limited state (wvi- = zvi + hvi-)
513        #
514        #
515        #where i=0,1,2 denotes the vertex ids
516        #
517        #Weighted balance between w-limited and h-limited stage is
518        #
519        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
520        #
521        #It follows that the updated wvi is
522        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
523        #
524        # Momentum is balanced between constant and limited
525
526
527        #for i in range(3):
528        #    wv[k,i] = zv[k,i] + hvbar[k,i]
529
530        #return
531
532        if alpha < 1:
533
534            #for i in range(3):
535            for i in range(2):
536                wv[k,i] = zv[k,i] + (1.0-alpha)*hvbar[k,i] + alpha*hv[k,i]
537
538            #Momentums at centroids
539            xmomc = domain.quantities['xmomentum'].centroid_values
540#            ymomc = domain.quantities['ymomentum'].centroid_values
541
542            #Momentums at vertices
543            xmomv = domain.quantities['xmomentum'].vertex_values
544#           ymomv = domain.quantities['ymomentum'].vertex_values
545
546            # Update momentum as a linear combination of
547            # xmomc and ymomc (shallow) and momentum
548            # from extrapolator xmomv and ymomv (deep).
549            xmomv[k,:] = (1.0-alpha)*xmomc[k] + alpha*xmomv[k,:]
550#            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
551
552
553
Note: See TracBrowser for help on using the repository browser.