Changeset 4815


Ignore:
Timestamp:
Nov 14, 2007, 5:43:02 PM (16 years ago)
Author:
ole
Message:

Work towards Froude_number adjusted, balanced limiters. This is currently
only in place for tight_slope_limiters == 1.

Location:
anuga_core/source/anuga
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

    r4808 r4815  
    2424from anuga.fit_interpolate.fit import fit_to_mesh
    2525from anuga.config import points_file_block_line_size as default_block_line_size
     26from anuga.config import epsilon
    2627
    2728class Quantity:
     
    109110        (except for a filename or points, attributes (for now))
    110111        - see set_values for details
    111 
    112         Note that if two quantitites q1 and q2 are multiplied,
    113         vertex values are multiplied entry by entry
    114         while centroid and edge values are re-interpolated.
    115         Hence they won't be the product of centroid or edge values
    116         from q1 and q2.
    117         """
    118 
    119         Q = Quantity(self.domain)
    120         Q.set_values(other)
     112        """
     113
     114        if isinstance(other, Quantity):
     115            Q = other
     116        else:   
     117            Q = Quantity(self.domain)
     118            Q.set_values(other)
    121119
    122120        result = Quantity(self.domain)
    123         result.set_values(self.vertex_values * Q.vertex_values)
     121
     122        # The product of vertex_values, edge_values and centroid_values
     123        # are calculated and assigned directly without using
     124        # set_values (which calls interpolate). Otherwise
     125        # edge and centroid values wouldn't be products from q1 and q2
     126        result.vertex_values = self.vertex_values * Q.vertex_values
     127        result.edge_values = self.edge_values * Q.edge_values
     128        result.centroid_values = self.centroid_values * Q.centroid_values
     129       
    124130        return result
    125131
     
    129135        return self * other
    130136
     137    def __div__(self, other):
     138        """Divide self with anything that could populate a quantity
     139
     140        E.g other can be a constant, an array, a function, another quantity
     141        (except for a filename or points, attributes (for now))
     142        - see set_values for details
     143
     144        Zero division is dealt with by adding an epsilon to the divisore
     145        FIXME (Ole): Replace this with native INF once we migrate to NumPy
     146        """
     147
     148        if isinstance(other, Quantity):
     149            Q = other
     150        else:   
     151            Q = Quantity(self.domain)
     152            Q.set_values(other)
     153
     154        result = Quantity(self.domain)
     155
     156        # The quotient of vertex_values, edge_values and centroid_values
     157        # are calculated and assigned directly without using
     158        # set_values (which calls interpolate). Otherwise
     159        # edge and centroid values wouldn't be quotient of q1 and q2
     160        result.vertex_values = self.vertex_values/(Q.vertex_values + epsilon)
     161        result.edge_values = self.edge_values/(Q.edge_values + epsilon)
     162        result.centroid_values = self.centroid_values/(Q.centroid_values + epsilon)
     163
     164        return result
     165
     166    def __rdiv__(self, other):
     167        """Handle cases like 3/Q, where Q is an instance of class Quantity
     168        """
     169        return self / other
     170
    131171    def __pow__(self, other):
    132172        """Raise quantity to (numerical) power
     
    140180        """
    141181
     182        if isinstance(other, Quantity):
     183            Q = other
     184        else:   
     185            Q = Quantity(self.domain)
     186            Q.set_values(other)
     187
    142188        result = Quantity(self.domain)
    143         result.set_values(self.vertex_values**other)
     189
     190        # The power of vertex_values, edge_values and centroid_values
     191        # are calculated and assigned directly without using
     192        # set_values (which calls interpolate). Otherwise
     193        # edge and centroid values wouldn't be correct
     194        result.vertex_values = self.vertex_values ** other
     195        result.edge_values = self.edge_values ** other
     196        result.centroid_values = self.centroid_values ** other
     197
    144198        return result
    145199
    146 
     200    #def __sqrt__(self, other):
     201    #    """Define in terms of x**0.5
     202    #    """
     203    #    pass
     204       
    147205
    148206    def interpolate(self):
  • anuga_core/source/anuga/config.py

    r4805 r4815  
    153153
    154154
     155maximum_froude_number = 100.0 # To be used in limiters.
     156
    155157minimum_storable_height = 1.0e-5 # Water depth below which it is *stored* as 0
    156158
  • anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r4769 r4815  
    141141
    142142
    143 void adjust_froude_number(double *uh,
    144                           double h,
    145                           double g) {
     143double compute_froude_number(double uh,
     144                             double h,
     145                             double g,
     146                             double epsilon) {
    146147                         
    147   // Adjust momentum if Froude number is excessive
    148   double max_froude_number = 20.0;                       
     148  // Compute Froude number; v/sqrt(gh)
     149 
    149150  double froude_number;
    150151 
    151152  //Compute Froude number (stability diagnostics)
    152   froude_number = *uh/sqrt(g*h)/h;
    153 
    154   if (froude_number > max_froude_number) {
    155     printf("---------------------------------------------\n");
    156     printf("froude_number=%f (uh=%f, h=%f)\n", froude_number, *uh, h);
    157    
    158     *uh = *uh/fabs(*uh) * max_froude_number * sqrt(g*h)*h;
    159    
    160     froude_number = *uh/sqrt(g*h)/h;   
    161     printf("Adjusted froude_number=%f (uh=%f, h=%f)\n", froude_number, *uh, h);
    162     printf("---------------------------------------------\n");   
    163   }
     153  if (h > epsilon) {
     154    froude_number = uh/sqrt(g*h)/h;
     155  } else {
     156    froude_number = 0.0;
     157    // FIXME (Ole): What should it be when dry??
     158  }
     159 
     160  return froude_number;
    164161}
    165162
     
    176173  double u;
    177174
    178   //adjust_froude_number(uh, *h, 9.81); // Highly experimental and
    179                                         // probably unneccessary
    180175 
    181176  if (*h < epsilon) {
     
    484479                              double alpha_balance) {
    485480
    486   int k, k3, i;
     481  int k, k3, i, excessive_froude_number=0;
    487482  double dz, hmin, alpha, h_diff, hc_k;
    488   double epsilon = 1.0e-6; // Temporary measure
    489   double hv[3]; // Depths at vertices
     483  double epsilon = 1.0e-6; // FIXME: Temporary measure
     484  double g = 9.81; // FIXME: Temporary measure
     485  double hv[3], h; // Depths at vertices
     486  double Fx, Fy; // Froude numbers
    490487
    491488  // Compute linear combination between w-limited stages and
     
    511508    }
    512509
    513     // Calculate depth at vertices
     510    // Calculate depth at vertices (possibly negative here!)
    514511    hv[0] = wv[k3] -   zv[k3];
    515512    hv[1] = wv[k3+1] - zv[k3+1];
     
    606603    //   Momentum is balanced between constant and limited
    607604
    608     if (alpha < 1) {
     605
     606    if (alpha < 1) {     
    609607      for (i=0; i<3; i++) {
    610      
    611608        // FIXME (Ole): Simplify when (if) hvbar gets retired       
    612609        if (beta_h > epsilon) {   
     
    622619        xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
    623620        ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
     621      }
     622    }
     623       
     624
     625       
     626    if (tight_slope_limiters == 1) {                   
     627   
     628      // Ensure that the Froude number is kept realistic at vertices
     629      // FIXME (Ole): I think it could be used to adjust alpha down
     630      // whenever Fr is too large. Possible make sure it doesn't deviate
     631      // too much from the value at the centroid. I like this idea!
     632     
     633      // FIXME (Ole): currently only used with tights_SL
     634
     635      excessive_froude_number=0;   
     636      for (i=0; i<3; i++) {   
     637     
     638        // Recalculate depth at vertex i
     639        h = wv[k3+i] - zv[k3+i];
     640       
     641        Fx = compute_froude_number(xmomv[k3+i], h, g, epsilon);
     642        Fy = compute_froude_number(ymomv[k3+i], h, g, epsilon);
     643       
     644        if ( (fabs(Fx) > 100.0) || (fabs(Fy) > 100.0)) {
     645          // FIXME: use max_froude - or base it on centroid value of F
     646          // printf("Excessive Froude number detected: %f or %f\n", Fx, Fy);
     647          excessive_froude_number=1;         
     648        }
     649      }
     650     
     651      if (excessive_froude_number) {
     652   
     653        // printf("Adjusting momentum to first order.\n");
     654        // Go back to first order (alpha = 0) for this triangle
     655        for (i=0; i<3; i++) {         
     656          xmomv[k3+i] = xmomc[k];
     657          ymomv[k3+i] = ymomc[k];
     658         
     659          if (beta_h > epsilon) {         
     660            wv[k3+i] = zv[k3+i] + hvbar[k3+i];
     661          } else {
     662            wv[k3+i] = zv[k3+i] + hc_k;
     663          }
     664        }       
    624665      }
    625666    }
  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r4805 r4815  
    3333    def __call__(self, x, y):
    3434        from Numeric import zeros, Float
    35         from math import sqrt
    3635
    3736        N = len(x)
     
    176175    """
    177176
    178     from math import sqrt, exp, cos, pi
     177    from math import exp, cos, pi
    179178
    180179    x = array(x)
     
    206205    """Rotating field
    207206    """
    208     from math import sqrt, atan, pi
     207    from math import atan, pi
    209208
    210209    x = array(x)
     
    366365        assert allclose(max_speed, 1.31414103233)
    367366
    368     def test_flux_computation(self):   
     367    def test_flux_computation(self):   
    369368        """test_flux_computation - test flux calculation (actual C implementation)
    370         This one tests the constant case where only the pressure term contributes to each edge and cancels out
    371         once the total flux has been summed up.
    372         """
    373                
     369        This one tests the constant case where only the pressure term contributes to each edge and cancels out
     370        once the total flux has been summed up.
     371        """
     372               
    374373        a = [0.0, 0.0]
    375374        b = [0.0, 2.0]
     
    386385        domain.check_integrity()
    387386
    388         # The constant case             
    389         domain.set_quantity('elevation', -1)
    390         domain.set_quantity('stage', 1)
    391        
    392         domain.compute_fluxes()
    393         assert allclose(domain.get_quantity('stage').explicit_update[1], 0) # Central triangle
    394        
    395 
    396         # The more general case                 
     387        # The constant case             
     388        domain.set_quantity('elevation', -1)
     389        domain.set_quantity('stage', 1)
     390       
     391        domain.compute_fluxes()
     392        assert allclose(domain.get_quantity('stage').explicit_update[1], 0) # Central triangle
     393       
     394
     395        # The more general case                 
    397396        def surface(x,y):
    398397            return -x/2                   
    399        
    400         domain.set_quantity('elevation', -10)
    401         domain.set_quantity('stage', surface)   
    402         domain.set_quantity('xmomentum', 1)             
    403        
    404         domain.compute_fluxes()
    405        
    406         #print domain.get_quantity('stage').explicit_update
    407         # FIXME (Ole): TODO the general case
    408         #assert allclose(domain.get_quantity('stage').explicit_update[1], ........??)
    409                
    410        
    411                
     398       
     399        domain.set_quantity('elevation', -10)
     400        domain.set_quantity('stage', surface)   
     401        domain.set_quantity('xmomentum', 1)             
     402       
     403        domain.compute_fluxes()
     404       
     405        #print domain.get_quantity('stage').explicit_update
     406        # FIXME (Ole): TODO the general case
     407        #assert allclose(domain.get_quantity('stage').explicit_update[1], ........??)
     408               
     409       
     410               
    412411    def test_sw_domain_simple(self):
    413412        a = [0.0, 0.0]
     
    17511750    def test_constant_wind_stress(self):
    17521751        from anuga.config import rho_a, rho_w, eta_w
    1753         from math import pi, cos, sin, sqrt
     1752        from math import pi, cos, sin
    17541753
    17551754        a = [0.0, 0.0]
     
    18031802    def test_variable_wind_stress(self):
    18041803        from anuga.config import rho_a, rho_w, eta_w
    1805         from math import pi, cos, sin, sqrt
     1804        from math import pi, cos, sin
    18061805
    18071806        a = [0.0, 0.0]
     
    18721871    def test_windfield_from_file(self):
    18731872        from anuga.config import rho_a, rho_w, eta_w
    1874         from math import pi, cos, sin, sqrt
     1873        from math import pi, cos, sin
    18751874        from anuga.config import time_format
    18761875        from anuga.abstract_2d_finite_volumes.util import file_function
     
    19821981    def test_windfield_from_file_seconds(self):
    19831982        from anuga.config import rho_a, rho_w, eta_w
    1984         from math import pi, cos, sin, sqrt
     1983        from math import pi, cos, sin
    19851984        from anuga.config import time_format
    19861985        from anuga.abstract_2d_finite_volumes.util import file_function
     
    20962095
    20972096        from anuga.config import rho_a, rho_w, eta_w
    2098         from math import pi, cos, sin, sqrt
     2097        from math import pi, cos, sin
    20992098
    21002099        a = [0.0, 0.0]
     
    25862585    def test_balance_deep_and_shallow(self):
    25872586        """Test that balanced limiters preserve conserved quantites.
     2587        This test is using old depth based balanced limiters
    25882588        """
    25892589        import copy
     
    26012601        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
    26022602
    2603         mesh = Domain(points, elements)
    2604         mesh.check_integrity()
     2603        domain = Domain(points, elements)
     2604        domain.check_integrity()
    26052605
    26062606        #Create a deliberate overshoot
    2607         mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
    2608         mesh.set_quantity('elevation', 0) #Flat bed
    2609         stage = mesh.quantities['stage']
     2607        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
     2608        domain.set_quantity('elevation', 0) #Flat bed
     2609        stage = domain.quantities['stage']
    26102610
    26112611        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
    26122612
    26132613        #Limit
    2614         mesh.distribute_to_vertices_and_edges()
     2614        domain.tight_slope_limiters = 0               
     2615        domain.distribute_to_vertices_and_edges()
    26152616
    26162617        #Assert that quantities are conserved
    26172618        from Numeric import sum
    2618         for k in range(len(mesh)):
     2619        for k in range(len(domain)):
    26192620            assert allclose (ref_centroid_values[k],
    26202621                             sum(stage.vertex_values[k,:])/3)
     
    26232624        #Now try with a non-flat bed - closely hugging initial stage in places
    26242625        #This will create alphas in the range [0, 0.478260, 1]
    2625         mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
    2626         mesh.set_quantity('elevation', [[0,0,0],
     2626        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
     2627        domain.set_quantity('elevation', [[0,0,0],
    26272628                                        [1.8,1.9,5.9],
    26282629                                        [4.6,0,0],
    26292630                                        [0,2,4]])
    2630         stage = mesh.quantities['stage']
     2631        stage = domain.quantities['stage']
    26312632
    26322633        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
     
    26342635
    26352636        #Limit
    2636         mesh.distribute_to_vertices_and_edges()
     2637        domain.tight_slope_limiters = 0       
     2638        domain.distribute_to_vertices_and_edges()
    26372639
    26382640
    26392641        #Assert that all vertex quantities have changed
    2640         for k in range(len(mesh)):
     2642        for k in range(len(domain)):
    26412643            #print ref_vertex_values[k,:], stage.vertex_values[k,:]
    26422644            assert not allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])
    26432645        #and assert that quantities are still conserved
    26442646        from Numeric import sum
    2645         for k in range(len(mesh)):
     2647        for k in range(len(domain)):
     2648            assert allclose (ref_centroid_values[k],
     2649                             sum(stage.vertex_values[k,:])/3)
     2650
     2651
     2652        # Check actual results
     2653        assert allclose (stage.vertex_values,
     2654                         [[2,2,2],
     2655                          [1.93333333, 2.03333333, 6.03333333],
     2656                          [6.93333333, 4.53333333, 4.53333333],
     2657                          [5.33333333, 5.33333333, 5.33333333]])
     2658
     2659
     2660    def test_balance_deep_and_shallow_tight_SL(self):
     2661        """Test that balanced limiters preserve conserved quantites.
     2662        This test is using Tight Slope Limiters
     2663        """
     2664        import copy
     2665
     2666        a = [0.0, 0.0]
     2667        b = [0.0, 2.0]
     2668        c = [2.0, 0.0]
     2669        d = [0.0, 4.0]
     2670        e = [2.0, 2.0]
     2671        f = [4.0, 0.0]
     2672
     2673        points = [a, b, c, d, e, f]
     2674
     2675        #bac, bce, ecf, dbe
     2676        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
     2677
     2678        domain = Domain(points, elements)
     2679        domain.check_integrity()
     2680
     2681        #Create a deliberate overshoot
     2682        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
     2683        domain.set_quantity('elevation', 0) #Flat bed
     2684        stage = domain.quantities['stage']
     2685
     2686        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
     2687
     2688        #Limit
     2689        domain.tight_slope_limiters = 1               
     2690        domain.distribute_to_vertices_and_edges()
     2691
     2692        #Assert that quantities are conserved
     2693        from Numeric import sum
     2694        for k in range(len(domain)):
     2695            assert allclose (ref_centroid_values[k],
     2696                             sum(stage.vertex_values[k,:])/3)
     2697
     2698
     2699        #Now try with a non-flat bed - closely hugging initial stage in places
     2700        #This will create alphas in the range [0, 0.478260, 1]
     2701        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
     2702        domain.set_quantity('elevation', [[0,0,0],
     2703                                        [1.8,1.9,5.9],
     2704                                        [4.6,0,0],
     2705                                        [0,2,4]])
     2706        stage = domain.quantities['stage']
     2707
     2708        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
     2709        ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
     2710
     2711        #Limit
     2712        domain.tight_slope_limiters = 1       
     2713        domain.distribute_to_vertices_and_edges()
     2714
     2715
     2716        #Assert that all vertex quantities have changed
     2717        for k in range(len(domain)):
     2718            #print ref_vertex_values[k,:], stage.vertex_values[k,:]
     2719            assert not allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])
     2720        #and assert that quantities are still conserved
     2721        from Numeric import sum
     2722        for k in range(len(domain)):
    26462723            assert allclose (ref_centroid_values[k],
    26472724                             sum(stage.vertex_values[k,:])/3)
     
    26572734        #                  [5.33333333, 5.33333333, 5.33333333]])
    26582735
     2736
     2737
     2738    def test_balance_deep_and_shallow_Froude(self):
     2739        """Test that balanced limiters preserve conserved quantites -
     2740        and also that excessive Froude numbers are dealt with.
     2741        This test is using tight slope limiters.
     2742        """
     2743        import copy
     2744        from Numeric import sqrt, absolute
     2745
     2746        a = [0.0, 0.0]
     2747        b = [0.0, 2.0]
     2748        c = [2.0, 0.0]
     2749        d = [0.0, 4.0]
     2750        e = [2.0, 2.0]
     2751        f = [4.0, 0.0]
     2752
     2753        points = [a, b, c, d, e, f]
     2754
     2755        # bac, bce, ecf, dbe
     2756        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
     2757
     2758        domain = Domain(points, elements)
     2759        domain.check_integrity()
     2760
     2761        # Create non-flat bed - closely hugging initial stage in places
     2762        # This will create alphas in the range [0, 0.478260, 1]
     2763        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
     2764        domain.set_quantity('elevation', [[0,0,0],
     2765                                        [1.8,1.999,5.999],
     2766                                        [4.6,0,0],
     2767                                        [0,2,4]])
     2768
     2769        # Create small momenta, that nonetheless will generate large speeds
     2770        # due to shallow depth at isolated vertices
     2771        domain.set_quantity('xmomentum', -0.0058)
     2772        domain.set_quantity('ymomentum', 0.0890)
     2773
     2774
     2775
     2776       
     2777        stage = domain.quantities['stage']
     2778        elevation = domain.quantities['elevation']
     2779        xmomentum = domain.quantities['xmomentum']
     2780        ymomentum = domain.quantities['ymomentum']       
     2781
     2782        # Setup triangle #1 to mimick real Froude explosion observed
     2783        # in the Onslow example 13 Nov 2007.
     2784
     2785        stage.vertex_values[1,:] = [1.6385, 1.6361, 1.2953]
     2786        elevation.vertex_values[1,:] = [1.6375, 1.6336, 0.4647]       
     2787        xmomentum.vertex_values[1,:] = [-0.0058, -0.0050, -0.0066]
     2788        ymomentum.vertex_values[1,:] = [0.0890, 0.0890, 0.0890]
     2789
     2790        xmomentum.interpolate()
     2791        ymomentum.interpolate()       
     2792        stage.interpolate()       
     2793        elevation.interpolate()
     2794
     2795        # Verify interpolation
     2796        assert allclose(stage.centroid_values[1], 1.5233)
     2797        assert allclose(elevation.centroid_values[1], 1.2452667)
     2798        assert allclose(xmomentum.centroid_values[1], -0.0058)
     2799        assert allclose(ymomentum.centroid_values[1], 0.089)
     2800
     2801        # Derived quantities
     2802        depth = stage-elevation
     2803        u = xmomentum/depth
     2804        v = ymomentum/depth
     2805
     2806        denom = (depth*g)**0.5
     2807        Fx = u/denom
     2808        Fy = v/denom
     2809       
     2810   
     2811        # Verify against Onslow example (14 Nov 2007)
     2812        assert allclose(depth.centroid_values[1], 0.278033)
     2813        assert allclose(u.centroid_values[1], -0.0208608)
     2814        assert allclose(v.centroid_values[1], 0.3201055)
     2815
     2816        assert allclose(denom.centroid_values[1],
     2817                        sqrt(depth.centroid_values[1]*g))
     2818
     2819        assert allclose(u.centroid_values[1]/denom.centroid_values[1],
     2820                        -0.012637746977)
     2821        assert allclose(Fx.centroid_values[1],
     2822                        u.centroid_values[1]/denom.centroid_values[1])
     2823
     2824        # Check that Froude numbers are small at centroids.
     2825        assert allclose(Fx.centroid_values[1], -0.012637746977)
     2826        assert allclose(Fy.centroid_values[1], 0.193924048435)
     2827
     2828
     2829        # But Froude numbers are huge at some vertices and edges
     2830        assert allclose(Fx.vertex_values[1,:], [-5.85888475e+01,
     2831                                                -1.27775313e+01,
     2832                                                -2.78511420e-03])
     2833
     2834        assert allclose(Fx.edge_values[1,:], [-6.89150773e-03,
     2835                                              -7.38672488e-03,
     2836                                              -2.35626238e+01])
     2837
     2838        assert allclose(Fy.vertex_values[1,:], [8.99035764e+02,
     2839                                                2.27440057e+02,
     2840                                                3.75568430e-02])
     2841
     2842        assert allclose(Fy.edge_values[1,:], [1.05748998e-01,
     2843                                              1.06035244e-01,
     2844                                              3.88346947e+02])
     2845       
     2846       
     2847        # The task is now to arrange the limiters such that Froude numbers
     2848        # remain under control whil at the same time obeying the conservation
     2849        # laws.
     2850
     2851       
     2852        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
     2853        ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
     2854
     2855        # Limit (and invoke balance_deep_and_shallow)
     2856        domain.tight_slope_limiters = 1
     2857        domain.distribute_to_vertices_and_edges()
     2858
     2859        # Redo derived quantities
     2860        depth = stage-elevation
     2861        u = xmomentum/depth
     2862        v = ymomentum/depth
     2863
     2864        denom = (depth*g)**0.5
     2865        Fx = u/denom
     2866        Fy = v/denom
     2867
     2868
     2869        # Assert that Froude numbers are less than max value (TBA)
     2870        # at vertices, edges and centroids.
     2871        from anuga.config import maximum_froude_number
     2872        assert alltrue(absolute(Fx.vertex_values[1,:]) < maximum_froude_number)
     2873        assert alltrue(absolute(Fy.vertex_values[1,:]) < maximum_froude_number)
     2874
     2875
     2876        # Assert that all vertex quantities have changed
     2877        for k in range(len(domain)):
     2878            #print ref_vertex_values[k,:], stage.vertex_values[k,:]
     2879            assert not allclose (ref_vertex_values[k,:],
     2880                                 stage.vertex_values[k,:])
     2881           
     2882        # Assert that quantities are still conserved
     2883        from Numeric import sum
     2884        for k in range(len(domain)):
     2885            assert allclose (ref_centroid_values[k],
     2886                             sum(stage.vertex_values[k,:])/3)
     2887
     2888
     2889       
     2890        return
     2891   
     2892        qwidth = 12
     2893        for k in [1]: #range(len(domain)):
     2894            print 'Triangle %d (C, V, E)' %k
     2895           
     2896            print 'stage'.ljust(qwidth), stage.centroid_values[k],\
     2897                  stage.vertex_values[k,:], stage.edge_values[k,:]
     2898            print 'elevation'.ljust(qwidth), elevation.centroid_values[k],\
     2899                  elevation.vertex_values[k,:], elevation.edge_values[k,:]
     2900            print 'depth'.ljust(qwidth), depth.centroid_values[k],\
     2901                  depth.vertex_values[k,:], depth.edge_values[k,:]
     2902            print 'xmomentum'.ljust(qwidth), xmomentum.centroid_values[k],\
     2903                  xmomentum.vertex_values[k,:], xmomentum.edge_values[k,:]
     2904            print 'ymomentum'.ljust(qwidth), ymomentum.centroid_values[k],\
     2905                  ymomentum.vertex_values[k,:], ymomentum.edge_values[k,:]
     2906            print 'u'.ljust(qwidth),u.centroid_values[k],\
     2907                  u.vertex_values[k,:], u.edge_values[k,:]
     2908            print 'v'.ljust(qwidth), v.centroid_values[k],\
     2909                  v.vertex_values[k,:], v.edge_values[k,:]
     2910            print 'Fx'.ljust(qwidth), Fx.centroid_values[k],\
     2911                  Fx.vertex_values[k,:], Fx.edge_values[k,:]
     2912            print 'Fy'.ljust(qwidth), Fy.centroid_values[k],\
     2913                  Fy.vertex_values[k,:], Fy.edge_values[k,:]
     2914           
     2915           
     2916       
    26592917
    26602918
     
    53065564
    53075565    suite = unittest.makeSuite(Test_Shallow_Water,'test')
     5566
     5567    #suite = unittest.makeSuite(Test_Shallow_Water,'test_balance_deep_and_shallow_Froude')
     5568   
    53085569    #suite = unittest.makeSuite(Test_Shallow_Water,'test_fitting_using_shallow_water_domain')   
    53095570    #suite = unittest.makeSuite(Test_Shallow_Water,'test_tight_slope_limiters')
Note: See TracChangeset for help on using the changeset viewer.