Changeset 4218


Ignore:
Timestamp:
Feb 6, 2007, 11:47:41 AM (18 years ago)
Author:
ole
Message:

Implemented new Flux Limiter which I believe is more robust and will pave the way for new balanced slope limiters.
For backwards compatibility using fluxes as they were before today (5/2/7) set
domain.H0 = 0 instead of the new default which is minimum_allowable_height (0.001).

Location:
anuga_core/source/anuga/shallow_water
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r4200 r4218  
    147147
    148148        self.minimum_allowed_height = minimum_allowed_height
     149        self.H0 = minimum_allowed_height  # Minimal height for flux computation       
    149150        self.maximum_allowed_speed = maximum_allowed_speed
    150151        self.g = g
    … …  
    832833         compute_fluxes_ext_central as compute_fluxes_ext
    833834   
    834     domain.timestep = compute_fluxes_ext(timestep, domain.epsilon, domain.g,
     835
     836    domain.timestep = compute_fluxes_ext(timestep, domain.epsilon,
     837                                         domain.H0,                                         
     838                                         domain.g,
    835839                                         domain.neighbours,
    836840                                         domain.neighbour_edges,
    … …  
    20822086                    protect_against_infinitesimal_and_negative_heights_c
    20832087
    2084 
    20852088    #distribute_to_vertices_and_edges =\
    20862089    #              distribute_to_vertices_and_edges_c #(like MH's)

  • anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r4200 r4218  
    129129                  double z_left, double z_right,
    130130                  double n1, double n2,
    131                   double epsilon, double g,
     131                  double epsilon, double H0, double g,
    132132                  double *edgeflux, double *max_speed) {
    133133
    … …  
    152152  double flux_right[3], flux_left[3];
    153153
     154  double h0 = H0*H0; //This ensures a good balance when h approaches H0.
     155 
    154156  //Copy conserved quantities to protect from modification
    155157  for (i=0; i<3; i++) {
    … …  
    173175    u_left = 0.0;
    174176  } else {
    175     u_left = uh_left/h_left;
     177    u_left = uh_left/(h_left + h0/h_left);
    176178  }
    177179
    … …  
    184186    u_right = 0.0;
    185187  } else {
    186     u_right = uh_right/h_right;
     188    u_right = uh_right/(h_right + h0/h_right);
    187189  }
    188190
    … …  
    249251
    250252// Computational function for flux computation (using stage w=z+h)
     253// FIXME (Ole): Is this used anywhere??
    251254int flux_function_kinetic(double *q_left, double *q_right,
    252255                  double z_left, double z_right,
    253256                  double n1, double n2,
    254                   double epsilon, double g,
     257                  double epsilon, double H0, double g,
    255258                  double *edgeflux, double *max_speed) {
    256259
    … …  
    270273  double q_left_copy[3], q_right_copy[3];
    271274
     275  double h0 = H0*H0; //This ensures a good balance when h approaches H0.
    272276
    273277  //Copy conserved quantities to protect from modification
    … …  
    292296    u_left = 0.0;
    293297  } else {
    294     u_left = uh_left/h_left;
     298    u_left = uh_left/(h_left + h0/h_left);
    295299  }
    296300
    … …  
    303307    u_right = 0.0;
    304308  } else {
    305     u_right = uh_right/h_right;
     309    u_right = uh_right/(h_right + h0/h_right);
    306310  }
    307311
    … …  
    11751179    domain.timestep = compute_fluxes(timestep,
    11761180                                     domain.epsilon,
     1181                                     domain.H0,
    11771182                                     domain.g,
    11781183                                     domain.neighbours,
    … …  
    12211226
    12221227  //Local variables
    1223   double timestep, max_speed, epsilon, g;
     1228  double timestep, max_speed, epsilon, g, H0;
    12241229  double normal[2], ql[3], qr[3], zl, zr;
    12251230  double edgeflux[3]; //Work arrays for summing up fluxes
    … …  
    12311236
    12321237  // Convert Python arguments to C
    1233   if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOOOO",
     1238  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOO",
    12341239                        &timestep,
    12351240                        &epsilon,
     1241                        &H0,
    12361242                        &g,
    12371243                        &neighbours,
    … …  
    13001306      flux_function_central(ql, qr, zl, zr,
    13011307                    normal[0], normal[1],
    1302                     epsilon, g,
     1308                    epsilon, H0, g,
    13031309                    edgeflux, &max_speed);
    13041310      //update triangle k
    … …  
    13631369    domain.timestep = compute_fluxes(timestep,
    13641370                                     domain.epsilon,
     1371                                     domain.H0,
    13651372                                     domain.g,
    13661373                                     domain.neighbours,
    … …  
    14081415
    14091416  //Local variables
    1410   double timestep, max_speed, epsilon, g;
     1417  double timestep, max_speed, epsilon, g, H0;
    14111418  double normal[2], ql[3], qr[3], zl, zr;
    14121419  double edgeflux[3]; //Work arrays for summing up fluxes
    … …  
    14181425
    14191426  // Convert Python arguments to C
    1420   if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOOO",
     1427  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOO",
    14211428                        &timestep,
    14221429                        &epsilon,
     1430                        &H0,
    14231431                        &g,
    14241432                        &neighbours,
    … …  
    14851493      flux_function_kinetic(ql, qr, zl, zr,
    14861494                    normal[0], normal[1],
    1487                     epsilon, g,
     1495                    epsilon, H0, g,
    14881496                    edgeflux, &max_speed);
    14891497      //update triangle k

  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r4026 r4218  
    103103            raise 'Should have raised an exception'
    104104
     105
     106    #FIXME (Ole): Individual flux tests do NOT test C implementation directly.   
    105107    def test_flux_zero_case(self):
    106108        ql = zeros( 3, Float )
    … …  
    594596
    595597
     598    # FIXME (Ole): Need test like this for fluxes in very shallow water.   
    596599    def test_compute_fluxes3(self):
    597600        #Random values, incl momentum
    … …  
    969972        domain = Domain(points, vertices, boundary) # Create domain
    970973        domain.set_quantities_to_be_stored(None)
    971         domain.set_maximum_allowed_speed(100) #
     974        domain.set_maximum_allowed_speed(100) #
     975        domain.H0 = 0 # Backwards compatibility (6/2/7)
    972976       
    973977
    … …  
    10501054
    10511055
     1056       
    10521057        assert allclose(gauge_values[0], G0)
    10531058        assert allclose(gauge_values[1], G1)
    … …  
    23532358        domain.beta_vh     = 0.9
    23542359        domain.beta_vh_dry = 0.9
     2360        domain.H0 = 0 # Backwards compatibility (6/2/7)       
    23552361       
    23562362        # Boundary conditions
    … …  
    24532459        domain.smooth = False
    24542460        domain.default_order=1
     2461        domain.H0 = 0 # Backwards compatibility (6/2/7)       
    24552462
    24562463        # Boundary conditions
    … …  
    25032510        domain.beta_vh_dry = 0.9       
    25042511        #domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance
     2512        domain.H0 = 0 # Backwards compatibility (6/2/7)       
    25052513
    25062514        # Boundary conditions
    … …  
    25662574        domain.beta_vh_dry = 0.9       
    25672575        domain.maximum_allowed_speed = 0.0 #Makes it like the 'oldstyle'
     2576        domain.H0 = 0 # Backwards compatibility (6/2/7)       
    25682577
    25692578        # Boundary conditions
    … …  
    26152624        domain.beta_uh_dry = 0.9
    26162625        domain.beta_vh     = 0.9
    2617         domain.beta_vh_dry = 0.9       
     2626        domain.beta_vh_dry = 0.9
     2627        domain.H0 = 0 # Backwards compatibility (6/2/7)       
    26182628
    26192629        # Boundary conditions
    … …  
    27992809        domain.default_order=1
    28002810        domain.beta_h = 0.0 #Use first order in h-limiter
     2811        domain.H0 = 0 # Backwards compatibility (6/2/7)       
    28012812
    28022813        #Bed-slope and friction
    … …  
    29602971        domain.beta_vh_dry = 0.9
    29612972        domain.beta_h = 0.0 #Use first order in h-limiter
     2973        domain.H0 = 0 # Backwards compatibility (6/2/7)       
    29622974
    29632975        #Bed-slope and friction at vertices (and interpolated elsewhere)
    … …  
    30583070        domain.beta_vh_dry = 0.9
    30593071        domain.beta_h = 0.0 #Use first order in h-limiter
     3072        domain.H0 = 0 # Backwards compatibility (6/2/7)       
    30603073
    30613074        #Bed-slope and friction at vertices (and interpolated elsewhere)
    … …  
    31523165        domain.beta_vh_dry = 0.9
    31533166        domain.beta_h = 0.0 #Use first order in h-limiter
     3167        domain.H0 = 0 # Backwards compatibility (6/2/7)       
    31543168
    31553169        #Bed-slope and friction at vertices (and interpolated elsewhere)
    … …  
    32553269
    32563270
     3271
     3272    def test_bedslope_problem_second_order_more_steps_feb_2007(self):
     3273        """test_bedslope_problem_second_order_more_steps_feb_2007
     3274
     3275        Test shallow water finite volumes, using parameters from feb 2007 rather
     3276        than backward compatibility ad infinitum
     3277       
     3278        """
     3279        from mesh_factory import rectangular
     3280        from Numeric import array
     3281
     3282        #Create basic mesh
     3283        points, vertices, boundary = rectangular(6, 6)
     3284
     3285        #Create shallow water domain
     3286        domain = Domain(points, vertices, boundary)
     3287        domain.smooth = False
     3288        domain.default_order = 2
     3289        domain.beta_w      = 0.9
     3290        domain.beta_w_dry  = 0.9
     3291        domain.beta_uh     = 0.9
     3292        domain.beta_uh_dry = 0.9
     3293        domain.beta_vh     = 0.9
     3294        domain.beta_vh_dry = 0.9
     3295        domain.beta_h = 0.0 #Use first order in h-limiter
     3296        domain.H0 = 0.001
     3297
     3298        #Bed-slope and friction at vertices (and interpolated elsewhere)
     3299        def x_slope(x, y):
     3300            return -x/3
     3301
     3302        domain.set_quantity('elevation', x_slope)
     3303
     3304        # Boundary conditions
     3305        Br = Reflective_boundary(domain)
     3306        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
     3307
     3308        #Initial condition
     3309        domain.set_quantity('stage', expression = 'elevation + 0.05')
     3310        domain.check_integrity()
     3311
     3312        assert allclose(domain.quantities['stage'].centroid_values,
     3313                        [0.01296296, 0.03148148, 0.01296296,
     3314                        0.03148148, 0.01296296, 0.03148148,
     3315                        0.01296296, 0.03148148, 0.01296296,
     3316                        0.03148148, 0.01296296, 0.03148148,
     3317                        -0.04259259, -0.02407407, -0.04259259,
     3318                        -0.02407407, -0.04259259, -0.02407407,
     3319                        -0.04259259, -0.02407407, -0.04259259,
     3320                        -0.02407407, -0.04259259, -0.02407407,
     3321                        -0.09814815, -0.07962963, -0.09814815,
     3322                        -0.07962963, -0.09814815, -0.07962963,
     3323                        -0.09814815, -0.07962963, -0.09814815,
     3324                        -0.07962963, -0.09814815, -0.07962963,
     3325                        -0.1537037 , -0.13518519, -0.1537037,
     3326                        -0.13518519, -0.1537037, -0.13518519,
     3327                        -0.1537037 , -0.13518519, -0.1537037,
     3328                        -0.13518519, -0.1537037, -0.13518519,
     3329                        -0.20925926, -0.19074074, -0.20925926,
     3330                        -0.19074074, -0.20925926, -0.19074074,
     3331                        -0.20925926, -0.19074074, -0.20925926,
     3332                        -0.19074074, -0.20925926, -0.19074074,
     3333                        -0.26481481, -0.2462963, -0.26481481,
     3334                        -0.2462963, -0.26481481, -0.2462963,
     3335                        -0.26481481, -0.2462963, -0.26481481,
     3336                        -0.2462963, -0.26481481, -0.2462963])
     3337
     3338
     3339        #print domain.quantities['stage'].extrapolate_second_order()
     3340        #domain.distribute_to_vertices_and_edges()
     3341        #print domain.quantities['stage'].vertex_values[:,0]
     3342
     3343        #Evolution
     3344        for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5):
     3345            pass
     3346
     3347
     3348        assert allclose(domain.quantities['stage'].centroid_values,
     3349     [-0.02907614, -0.01462948, -0.02971293, -0.01435568, -0.02930822, -0.0141715,
     3350      -0.02900256, -0.01402774, -0.02897007, -0.01407298, -0.02958336, -0.01393198,
     3351      -0.07599907, -0.06253965, -0.07666738, -0.06312924, -0.07639753, -0.06265574,
     3352      -0.07572658, -0.06235626, -0.07569544, -0.0624551,  -0.07653298, -0.06289883,
     3353      -0.1236842,  -0.11090096, -0.12238737, -0.11116539, -0.12190705, -0.11072785,
     3354      -0.1208281,  -0.11001521, -0.12039674, -0.11011275, -0.1210331,  -0.11013222,
     3355      -0.16910086, -0.1583269,  -0.16731364, -0.15787261, -0.16655826, -0.15698876,
     3356      -0.16497539, -0.15560667, -0.16339389, -0.15509626, -0.16364638, -0.15425332,
     3357      -0.18773952, -0.19904678, -0.18906095, -0.1985912,  -0.18703781, -0.19698407,
     3358      -0.18337965, -0.19506463, -0.18190047, -0.19418413, -0.18587491, -0.19576587,
     3359      -0.1398922,  -0.14172812, -0.14134412, -0.14563187, -0.14097668, -0.14375562,
     3360      -0.13787839, -0.14035428, -0.13594691, -0.13938278, -0.13597595, -0.14218274])
     3361                     
     3362
     3363        assert allclose(domain.quantities['xmomentum'].centroid_values,
     3364     [ 0.00831991,  0.00327584,  0.0073476,   0.0034367,   0.00767331,  0.00356348,
     3365       0.00791201,  0.00364529,  0.00783242,  0.00349728,  0.0069742,   0.0031689,
     3366       0.02165378,  0.01421696,  0.02016895,  0.01318091,  0.02036575,  0.01369801,
     3367       0.02105492,  0.01400313,  0.02074176,  0.01356261,  0.01887467,  0.01232703,
     3368       0.03774493,  0.02854758,  0.03688063,  0.02759182,  0.03731697,  0.02811662,
     3369       0.03871488,  0.02912926,  0.03880077,  0.02803529,  0.0354565,   0.02599893,
     3370       0.06319926,  0.04729427,  0.05761625,  0.04591713,  0.05789906,  0.0468986,
     3371       0.05985286,  0.04870638,  0.06169141,  0.04811799,  0.05656696,  0.04415568,
     3372       0.08488718,  0.07187458,  0.07834844,  0.06842993,  0.07985979,  0.06981954,
     3373       0.08200942,  0.07216429,  0.08378261,  0.07273359,  0.08040488,  0.06646477,
     3374       0.01631518,  0.04691654,  0.02066439,  0.04441294,  0.02115705,  0.04560776,
     3375       0.02160783,  0.04664204,  0.02174952,  0.0479616,   0.02281756,  0.05667927])
     3376
     3377
     3378        assert allclose(domain.quantities['ymomentum'].centroid_values,
     3379     [  1.49908296e-04,  -3.32118806e-04,  -1.55139120e-04,  -2.97772609e-04,
     3380       -9.57477241e-05,  -3.11011790e-04,  -1.58896911e-04,  -3.76997605e-04,
     3381       -1.97659519e-04,  -3.34831296e-04,   6.54935308e-05,  -8.43493883e-06,
     3382        5.05063334e-04,  -1.43305846e-04,  -6.76061638e-05,  -5.01728590e-04,
     3383       -8.39003270e-05,  -4.64804117e-04,  -1.95135035e-04,  -5.88384357e-04,
     3384       -2.69800068e-04,  -5.35718006e-04,   2.59588334e-04,   3.00642498e-05,
     3385        5.15520349e-04,   1.05711270e-04,   9.26087960e-05,  -3.71855339e-04,
     3386        1.16386690e-04,  -3.82345749e-04,  -1.61741416e-04,  -6.31090292e-04,
     3387       -4.74530925e-04,  -6.95229436e-04,   6.08967544e-05,   2.20974020e-04,
     3388       -6.30000309e-04,   2.42320158e-04,  -5.89290588e-04,  -7.03283441e-05,
     3389       -4.18421888e-04,   6.62703090e-05,  -7.68647846e-04,  -3.40294284e-04,
     3390       -1.67585557e-03,  -7.40500723e-04,  -1.60020305e-03,   5.62070746e-05,
     3391       -1.48807206e-03,  -1.84455791e-03,  -2.27365819e-03,  -1.67381169e-03,
     3392       -1.95607481e-03,  -1.47497442e-03,  -1.73851968e-03,  -1.85314716e-03,
     3393       -2.01489344e-03,  -2.17608206e-03,  -1.66072261e-03,  -1.15856505e-03,
     3394       -1.18717624e-03,  -2.94595857e-03,  -3.59685615e-03,  -5.13811671e-03,
     3395       -6.17481232e-03,  -5.98871894e-03,  -6.00593324e-03,  -5.01282532e-03,
     3396       -4.51124363e-03,  -3.06579417e-03,   6.07680580e-04,  -4.80786234e-04])
     3397
     3398        os.remove(domain.get_name() + '.sww')
     3399
     3400
    32573401    def test_temp_play(self):
    32583402
    … …  
    32743418        domain.beta_vh_dry = 0.9
    32753419        domain.beta_h = 0.0 #Use first order in h-limiter
     3420        domain.H0 = 0 # Backwards compatibility (6/2/7)       
    32763421
    32773422        #Bed-slope and friction at vertices (and interpolated elsewhere)
Note: See TracChangeset for help on using the changeset viewer.