Context Navigation

← Previous Change
Next Change →

shallow_water.py

Timestamp:

Aug 23, 2004, 2:45:42 PM (21 years ago)

Author:

ole

Message:

Finished and tested gradient and limiters
Added forcing terms (not tested yet)

File:

: 1 edited

inundation/ga/storm_surge/pyvolution/shallow_water.py (modified) (9 diffs)

Legend:

: Unmodified
: Added
: Removed

inundation/ga/storm_surge/pyvolution/shallow_water.py

-                      r198
+                      r205
         compute_fluxes(self)
     def first_order_limiter(self):
+    def distribute_to_vertices_and_edges(self):
         #Call correct module function
+        first_order_limiter(self)
+        #(either from this module or C-extension)
+        distribute_to_vertices_and_edges(self)
 ####################################
 …
         Xmom.explicit_update[k] = flux[1]
         Ymom.explicit_update[k] = flux[2]
         timestep = min(timestep, optimal_timestep)
     domain.timestep = timestep
 ####################################
+# Functions for gradient limiting (distribute_to_vertices_and_edges)
+# Module functions for gradient limiting (distribute_to_vertices_and_edges)
+def distribute_to_vertices_and_edges(domain):
+    protect_against_infinitesimal_heights_centroid(domain)
+    if domain.order == 1:
+        first_order_extrapolator(domain)
+    elif domain.order == 2:
+        second_order_extrapolator(domain)
+    else:
+        raise 'Unknown order'
 def protect_against_infinitesimal_heights_centroid(domain):
 …
+def first_order_limiter(domain):
+    """First order limiter function, specific to the shallow water wave
+    equation.
+def first_order_extrapolator(domain):
+    """First order extrapolator function, specific
+    to the shallow water wave equation.
     It will ensure that h (w-z) is always non-negative even in the
 …
     """
-    #FIXME: Might go elsewhere
-    protect_against_infinitesimal_heights_centroid(domain)
     #Update conserved quantities using straight first order
     for name in domain.conserved_quantities:
         Q = domain.quantities[name]
+        Q.first_order_limiter()
+        Q.first_order_extrapolator()
+        Q.interpolate_from_vertices_to_edges()
     #Water levels at centroids
     wc = domain.quantities['level'].centroid_values
 …
             alpha = 1.0
         #Update stage
+        #Update water levels
         for i in range(3):
             #wv[k,i] = zv[k,i] + hc[k] + alpha*(zc[k]-zv[k,i])
 …
 def second_order_limiter(domain):
+def second_order_extrapolator(domain):
     """Second order limiter function, specific to the shallow water wave
     equation.
 …
     """
+    #FIXME: Might go elsewhere
+    protect_against_infinitesimal_heights_centroid(domain)
+    #FIXME: first and second order migh merge
+    from Numeric import minimum, maximum
     #Update conserved quantities using straight second order
     for name in domain.conserved_quantities:
         Q = domain.quantities[name]
+        Q.second_order_limiter()
+        Q.second_order_extrapolator()
+        Q.limiter()
+        Q.interpolate_from_vertices_to_edges()
     #Water levels at centroids
     wc = domain.quantities['level'].centroid_values
 …
     hc = wc - zc
-    #Momentums at centroids
-    xmomc = domain.quantities['xmomentum'].centroid_values
-    ymomc = domain.quantities['ymomentum'].centroid_values
     #Water levels at vertices
     #wv = domain.quantities['level'].vertex_values
+    #
+    wv = domain.quantities['level'].vertex_values
     #Bed elevations at vertices
+    #zv = domain.quantities['elevation'].vertex_values
+    #
+    #Momentums at vertices
+    #xmomv = domain.quantities['xmomentum'].vertex_values
+    #ymomv = domain.quantities['ymomentum'].vertex_values
+    zv = domain.quantities['elevation'].vertex_values
+    #Water depths at vertices
+    hv = wv-zv
+    #Computed linear combination between constant levels and and
+    #levels parallel to the bed elevation.
+    for k in range(domain.number_of_elements):
+        #Compute maximal variation in bed elevation
+        #  This quantitiy is
+        #    dz = max_i abs(z_i - z_c)
+        #  and it is independent of dimension
+        #  In the 1d case zc = (z0+z1)/2
+        #  In the 2d case zc = (z0+z1+z2)/3
+        #
+        z_range = max(abs(zv[k,0]-zc[k]),
+                      abs(zv[k,1]-zc[k]),
+                      abs(zv[k,2]-zc[k]))
+        hmin = minimum( hv[k, :] )
+        #Create alpha in [0,1], where alpha==0 means using shallow
+        #first order scheme and alpha==1 means using the limited
+        #second order scheme for the stage w
+        #If hmin < 0 then alpha=0 reverrting to first order.
+        if z_range > 0.0:
+            alpha = max( min( 2*hmin/z_range, 1.0), 0.0)
+        else:
+            alpha = 1.0
+        #Update water levels
+        #  (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) =
+        #   zvi + hc + alpha*(hvi - hc)
+        for i in range(3):
+            wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
+        #Momentums at centroids
+        xmomc = domain.quantities['xmomentum'].centroid_values
+        ymomc = domain.quantities['ymomentum'].centroid_values
+        #Momentums at vertices
+        xmomv = domain.quantities['xmomentum'].vertex_values
+        ymomv = domain.quantities['ymomentum'].vertex_values
+        # Update momentum as a linear combination of
+        # xmomc and ymomc (shallow) and momentum
+        # from extrapolator xmomv and ymomv (deep).
+        ##f##or i in range(3):
+        #####    xmomv[k,i] = (1-alpha)*xmomc[k] + alpha*xmomv[k,i];
+        xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:];
+    #Finally, protect against infinitesimal heights and high speeds
+    #Water depths at vertices
+    hv = wv-zv
+    hc = wc-zc
     for k in range(domain.number_of_elements):
+        #Protect against infinitesimal heights and high velocities
+        if hc[k] < domain.minimum_allowed_height:
+            #Control level and height
+        hmax = maximum(hv[k,:])
+        if hmax < minimum_allowed_height:
+            #Reset negative heights to bed elevation
             if hc[k] < 0.0:
+                wc[k] = zc[k]; hc[k] = 0.0
+            #Control momentum as well!
+            xmomc[k] = ymomc[k] = 0.0
+    #Compute second order limiter for each conserved quantity
+    for name in domain.conserved_quantities:
+        domain.quantities[name].second_order_limiter()
+#   //Adjusted water heights at vertices
+#   hv0 = qv0[0]-zv0;
+#   hv1 = qv1[0]-zv1;
+#   hv2 = qv2[0]-zv2;
+#   hmin = min( min(hv0, hv1), hv2 );
+#   // -----------------------------------
+#   //  Second run - Linear combination with
+#   //  first order scheme for shallow depths
+#   // -----------------------------------
+#   //Compute maximal variation in bed elevation
+#   //
+#   //This quantitiy is
+#   // dz = max_i abs(z_i - z_c)
+#   //and it is independent of dimension
+#   //In the 1d case zc = (z0+z1)/2
+#   //In the 2d case zc = (z0+z1+z2)/3
+#   //
+#   z_range = max( max(fabs(zv0-zc), fabs(zv1-zc)), fabs(zv2-zc) );
+#   //Create alpha in [0,1], where alpha==0 means using shallow
+#   //first order scheme and alpha==1 means using the limited
+#   //second order scheme for the stage w
+#   //If hmin < 0 then alpha=0 and first order method is
+#   if (z_range>0.0) {
+#     alpha = max( min(2*hmin/z_range, 1.0), 0.0);
+#   } else {
+#     alpha = 1.0;
+#   }
+#   if (alpha < 1) {
+#     //Update stage
+#     //
+#     //(1-alpha)*(zvi+hc) + alpha*(zvi+hvi) =
+#     //zvi + hc + alpha*(hvi - hc)
+#     //
+#     qv0[0] = zv0 + hc + alpha*(hv0-hc);
+#     qv1[0] = zv1 + hc + alpha*(hv1-hc);
+#     qv2[0] = zv2 + hc + alpha*(hv2-hc);
+#     //Update momentum as a linear combination of
+#     //qc[j] (shallow) and momentum from gradient limiter (deep).
+#     for (j=1; j<number_of_conserved_quantities; j++) {
+#       qv0[j] = (1-alpha)*qc[j] + alpha*qv0[j];
+#       qv1[j] = (1-alpha)*qc[j] + alpha*qv1[j];
+#       qv2[j] = (1-alpha)*qc[j] + alpha*qv2[j];
+#     }
+#   }
+#   //Finally, protect against infinitesimal heights and high speeds
+#   hv0 = qv0[0]-zv0; hv1 = qv1[0]-zv1; hv2 = qv2[0]-zv2;
+#   hmax = max(hv0, max(hv1, hv2));
+#   if (hmax < minimum_allowed_height) {
+#     //Set negative heights to bed elevation
+#     if (hc < 0.0)  {qc[0] = zc; hc = 0.0;}
+#     if (hv0 < 0.0) {qv0[0] = zv0; hv0 = 0.0;}
+#     if (hv1 < 0.0) {qv1[0] = zv1; hv1 = 0.0;}
+#     if (hv2 < 0.0) {qv2[0] = zv2; hv2 = 0.0;}
+#     //FIXME: Do we need this desperately??
+#     //FIXME: We definitely do encounter negative heights,
+#     //but what about momentum?
+#     //Control momentum
+#     //qc[1] = 0.0; qc[2] = 0.0;
+#     //qv0[1] = 0.0; qv0[2] = 0.0;
+#     //qv1[1] = 0.0; qv1[2] = 0.0;
+#     //qv2[1] = 0.0; qv2[2] = 0.0;
+#   }
+# }
+def limiter_org(k):
+    """limiter function, specific to the shallow water wave equation.
+       It will ensure that h is always non-negative.
+       This version works on the consecutive data structure
+    """
+    #This version should be used to reproduce the bad limiting when levels
+    #are close to bed elevation
+    import sys
+    from config import beta
+    from Numeric import ones, Float
+    #conserved quantities at centroid
+    qc = Volume.conserved_quantities_centroid[k, :]
+    N = len(qc)
+    zc = Volume.field_values_centroid[k,0]         #bed elevation at centroid
+    wc = qc[0]                                     #Stage at centroid
+    hc = wc - zc                                   #Water depth at centroid
+    #conserved quant. at vertices
+    qv0 = Volume.conserved_quantities_vertex0[k,:]
+    qv1 = Volume.conserved_quantities_vertex1[k,:]
+    qv2 = Volume.conserved_quantities_vertex2[k,:]
+    #Bed elevation at vertices
+    zv0 = Volume.field_values_vertex0[k,0]
+    zv1 = Volume.field_values_vertex1[k,0]
+    zv2 = Volume.field_values_vertex2[k,0]
+    #deltas
+    dq0 = qv0-qc
+    dq1 = qv1-qc
+    dq2 = qv2-qc
+    #-------------------------------------
+    hmin = hc
+    qmin = qmax = qc
+    for i in range(3):
+        n = Volume.neighbours[k,i]
+        if n >= 0:
+            qn = Volume.conserved_quantities_centroid[n,:]
+            zn = Volume.field_values_centroid[n,0]
+            hn = qn[0]-zn
+            hmin = min(hmin, hn)
+            qmin = minimum(qmin, qn)
+            qmax = maximum(qmax, qn)
+    hmin = max(hmin,0.0)
+    #-----------------------------------
+    # First run - limit on conserved quantities
+    # based on neighbouring heights
+    #-----------------------------------
+    dqmax = qmax - qc
+    dqmin = qmin - qc
+    phi = ones(len(qc), Float)
+    for j in range(N):
+        #First find the gradient limiter (phi)
+        #Vertex0:
+        r = 1.0
+        if dq0[j] > 0: r = dqmax[j]/dq0[j]
+        if dq0[j] < 0: r = dqmin[j]/dq0[j]
+        phi[j] = min( min(r*beta, 1), phi[j] )
+        #Vertex1:
+        r = 1.0
+        if dq1[j] > 0: r = dqmax[j]/dq1[j]
+        if dq1[j] < 0: r = dqmin[j]/dq1[j]
+        phi[j] = min( min(r*beta, 1), phi[j] )
+        #Vertex2:
+        r = 1.0
+        if dq2[j] > 0: r = dqmax[j]/dq2[j]
+        if dq2[j] < 0: r = dqmin[j]/dq2[j]
+        phi[j] = min( min(r*beta, 1), phi[j] )
+        #The update using phi limiter
+        qv0[j] = qc[j] + phi[j]*dq0[j] #Vertex0
+        qv1[j] = qc[j] + phi[j]*dq1[j] #Vertex1
+        qv2[j] = qc[j] + phi[j]*dq2[j] #Vertex2
+    #-----------------------------------
+    # Second run - lower limit on height
+    # based on neighbouring heights for
+    # surfaces constructed below bed
+    #-----------------------------------
+    #Adjusted water heights at vertices
+    hv0 = qv0[0]-zv0
+    hv1 = qv1[0]-zv1
+    hv2 = qv2[0]-zv2
+    #Deltas
+    dhmin = hmin - hc
+    dh0 = hv0 - hc
+    dh1 = hv1 - hc
+    dh2 = hv2 - hc
+    #Work out if update is needed and calculate phi limiter
+    phi = 1.0
+    update_needed = False
+    if hv0 < 0.0:
+        update_needed = True
+        r = 1.0
+        if dh0 < 0: r = dhmin/dh0
+        phi = min( min(beta*r, 1), phi )
+    if hv1 < 0.0:
+        update_needed = True
+        r = 1.0
+        if dh1 < 0: r = dhmin/dh1
+        phi = min( min(beta*r, 1), phi )
+    if hv2 < 0.0:
+        update_needed = True
+        r = 1.0
+        if dh2 < 0: r = dhmin/dh2
+        phi = min( min(beta*r, 1), phi )
+    #Update stage at vertices
+    if update_needed:
+        hv0 = hc + phi*dh0
+        if hv0 < 0.0: raise 'hv0 = %12f hc = %12f ' % (hv0 , hc)
+        qv0[0] = zv0 + hv0
+        hv1 = hc + phi*dh1
+        if hv1 < 0.0: raise 'hv1 = %12f hc = %12f ' % (hv1 , hc)
+        qv1[0] = zv1 + hv1
+        hv2 = hc + phi*dh2
+        if hv2 < 0.0: raise 'hv2 = %12f hc = %12f ' % (hv2 , hc)
+        qv2[0] = zv2 + hv2
+    #In either case protect against infinitesimal heights
+    dry(k)
+def dry(k):
+    """Protect agains infinitesimal heights (and high speeds)
+    There are no states for any conserved quantity when h is
+    between minimum allowed height and 0.
+    """
+    #volume = Volume.instances[k]
+    maxhv = 0.0
+    wv0 = Volume.conserved_quantities_vertex0[k,0]
+    zv0 = Volume.field_values_vertex0[k,0]
+    hv0 = wv0-zv0
+    wv1 = Volume.conserved_quantities_vertex1[k,0]
+    zv1 = Volume.field_values_vertex1[k,0]
+    hv1 = wv1-zv1
+    wv2 = Volume.conserved_quantities_vertex2[k,0]
+    zv2 = Volume.field_values_vertex2[k,0]
+    hv2 = wv2-zv2
+    maxhv = max(hv0, hv1, hv2)
+    if maxhv < minimum_allowed_height:
+        #Set water level to bed elevation and zero momentums.
+        Volume.conserved_quantities_vertex0[k,:] =\
+                                                 array([Volume.field_values_vertex0[k,0], 0, 0])
+        Volume.conserved_quantities_vertex1[k,:] =\
+                                                 array([Volume.field_values_vertex1[k,0], 0, 0])
+        Volume.conserved_quantities_vertex2[k,:] =\
+                                                 array([Volume.field_values_vertex2[k,0], 0, 0])
+        #Set only momentums to zero
+        #Volume.conserved_quantities_vertex0[k,1:] = array([0., 0.])
+        #Volume.conserved_quantities_vertex1[k,1:] = array([0., 0.])
+        #Volume.conserved_quantities_vertex2[k,1:] = array([0., 0.])
+def distribute_to_vertices_and_edges(domain):
+    """Extrapolate conserved quantities from centroid to
+    vertices and edge-midpoints for each volume
+    """
+    from config import minimum_allowed_height
+    from mesh import Point
+    order = domain.order
+    limiter = domain.limiter
+    if limiter is None:
+        from shallow_water import limit_on_w_then_h as limiter
+    for k in range(Volume.number_of_instances):
+        #conserved quantities at centroid
+        q = Volume.conserved_quantities_centroid[k,:]
+        if order == 1:
+            ########################
+            #NOTE: This is now specific to the shallow water wave equation
+            zc = Volume.field_values_centroid[k,0] #bed elevation at centroid
+            wc = q[0]                              #Stage at centroid
+            hc = wc - zc                           #Water depth at centroid
+            #Protect against infinitesimal heights and high velocities
+            if hc < minimum_allowed_height:
+                #Set momentums to zero
+                q[0] = zc
+                q[1] = 0.0
+                q[2] = 0.0
+            #Bed elevation at vertices
+            zv0 = Volume.field_values_vertex0[k,0]
+            zv1 = Volume.field_values_vertex1[k,0]
+            zv2 = Volume.field_values_vertex2[k,0]
+            #Update conserved quant. at vertices
+            Volume.conserved_quantities_vertex0[k,1:] = q[1:] #Momentum
+            Volume.conserved_quantities_vertex1[k,1:] = q[1:] #Momentum
+            Volume.conserved_quantities_vertex2[k,1:] = q[1:] #Momentum
+            #Stage
+            z_range = max( max(abs(zv0-zv1), abs(zv0-zv2)), abs(zv1-zv2) )
+            if z_range>0:
+                alpha = min(hc/z_range, 1.0)
+            else:
+                alpha = 1
+            zz = hc+alpha*zc
+            Volume.conserved_quantities_vertex0[k,0] = zz + (1-alpha)*zv0
+            Volume.conserved_quantities_vertex1[k,0] = zz + (1-alpha)*zv1
+            Volume.conserved_quantities_vertex2[k,0] = zz + (1-alpha)*zv2
+        elif order == 2:
+            a, b = compute_gradient(k)
+            #Compute extrapolated values at vertices
+            x,y = Point.coordinates[Volume.points[k,3]] #Centroid
+            #Get vertex coordinates
+            i = Volume.points[k,0]
+            x0 = Point.coordinates[i,0]
+            y0 = Point.coordinates[i,1]
+            i = Volume.points[k,1]
+            x1 = Point.coordinates[i,0]
+            y1 = Point.coordinates[i,1]
+            i = Volume.points[k,2]
+            x2 = Point.coordinates[i,0]
+            y2 = Point.coordinates[i,1]
+            #Extrapolate
+            Volume.conserved_quantities_vertex0[k, :] = q + a*(x0-x) + b*(y0-y)
+            Volume.conserved_quantities_vertex1[k, :] = q + a*(x1-x) + b*(y1-y)
+            Volume.conserved_quantities_vertex2[k, :] = q + a*(x2-x) + b*(y2-y)
+            if limiter is not None:
+                limiter(k)
+        else:
+            raise 'Unknown order'
+        #Compute conserved quantities as interpolated
+        #quantities at each edge midpoint
+        q0 = Volume.conserved_quantities_vertex0[k,:]
+        q1 = Volume.conserved_quantities_vertex1[k,:]
+        q2 = Volume.conserved_quantities_vertex2[k,:]
+        Volume.conserved_quantities_face0[k,:] = 0.5*(q1 + q2)
+        Volume.conserved_quantities_face1[k,:] = 0.5*(q0 + q2)
+        Volume.conserved_quantities_face2[k,:] = 0.5*(q0 + q1)
+                wc[k] = zc[k]
+                ###hc[k] = 0.0
+            for i in range(3):
+                if hv[k,i] < 0.0:
+                    wv[k,i] = zv[k,i]
+                    ##hv0 = 0.0;}
 …
+#########################
+#Standard forcing terms:
+def gravity(domain):
+    """Implement forcing function for bed slope working with
+    consecutive data structures of class Volume
+    """
+    from config import g
+    from util import gradient
+    from Numeric import zeros, Float, array, sum
+    Level = domain.quantities['level']
+    Xmom = domain.quantities['xmomentum']
+    Ymom = domain.quantities['ymomentum']
+    Elevation = domain.quantities['elevation']
+    h = Level.edge_values - Elevation.edge_values
+    V = domain.get_vertex_coordinates()
+    for k in range(domain.number_of_elements):
+        avg_h = sum( h[k,:] )/3
+        #Compute bed slope
+        x0, y0, x1, y1, x2, y2 = V[k,:]
+        z0, z1, z2 = Elevation.vertex_values[k,:]
+        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
+        #Update momentum
+        Xmom.explicit_update[k] += -g*zx*avg_h
+        Ymom.explicit_update[k] += -g*zy*avg_h
+def manning_friction(domain):
+    """Apply (Manning) friction to water momentum
+    """
+    import Numeric
+    from math import sqrt
+    from config import g, minimum_allowed_height
+    Level = domain.quantities['level']
+    Xmom = domain.quantities['xmomentum']
+    Ymom = domain.quantities['ymomentum']
+    Friction = domain.quantities['friction']
+    for k in range(domain.number_of_elements):
+        w = Level.centroid_values[k]
+        uh = Xmom.centroid_values[k]
+        vh = Ymom.centroid_values[k]
+        manning = Friction.centroid_values[k]
+        if w >= minimum_allowed_height:
+            S = -g * manning**2 * sqrt((uh**2 + vh**2))
+            S /= w**(7.0/3)
+            #Update momentum
+            Xmom.explicit_update[k] += S*uh
+            Ymom.explicit_update[k] += S*vh
 ##############################################

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 205 for inundation/ga/storm_surge/pyvolution/shallow_water.py

Legend:

inundation/ga/storm_surge/pyvolution/shallow_water.py

Download in other formats: