Changeset 8450


Ignore:
Timestamp:
Jun 18, 2012, 9:43:05 PM (12 years ago)
Author:
davies
Message:

bal_dev: Changes to improve mass conservation

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/validation_tests/Tests/Analytical/trapezoidal_channel/results.tex

    r8430 r8450  
    2727\end{figure}
    2828
    29 Table~\ref{trapztab} shows the discharge computed at 3 cross-sections in the channel, at a number of time-steps on the way to near steady-state. By the end of the simulation they should all be essentially the same. Large variations may suggest mass conservation errors (small variations may be due to the way discharge is computed -- I don't think this will include the effect of velocity limiters during the flux computation?).
     29Table~\ref{trapztab} shows the discharge computed at 3 cross-sections in the channel, at a number of time-steps on the way to near steady-state. By the end of the simulation they should all be essentially the same. Large variations may suggest mass conservation errors (small variations are probably due to the interpolation that occurs in the 'compute\_flow\_through\_cross\_section' routine).
    3030
    3131\DTLloaddb{dischargeout}{Tests/Analytical/trapezoidal_channel/discharge_outputs.txt}
  • trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c

    r8446 r8450  
    443443
    444444            // Prevent outflow from 'seriously' dry cells
    445             if(stage_centroid_values[k]<min_bed_edgevalue[k]){
     445            if((stage_centroid_values[k]<=min_bed_edgevalue[k])|
     446               (ql[0]<=zl)){
    446447                if(edgeflux[0]>0.0){
    447448                    edgeflux[0]=0.0;
     
    451452            }
    452453            if(n>=0){
    453                 if(stage_centroid_values[n]<min_bed_edgevalue[n]){
     454                if((stage_centroid_values[n]<=min_bed_edgevalue[n])|
     455                   (qr[0]<=zr)){
    454456                    if(edgeflux[0]<0.0){
    455457                        edgeflux[0]=0.0;
     
    648650             //bmin =0.5*bmin + 0.5*min(zv[3*k],min(zv[3*k+1],zv[3*k+2]));
    649651             // Setting = minimum vertex value seems better, but might tend to be less smooth
    650              bmin =min(zv[3*k],min(zv[3*k+1],zv[3*k+2]));
     652             bmin =min(zv[3*k],min(zv[3*k+1],zv[3*k+2])) -minimum_allowed_height;
    651653             // Minimum allowed stage = bmin
    652654             // WARNING: ADDING MASS if wc[k]<bmin
    653655             if(wc[k]<bmin){
    654                  printf("Adding mass to dry cell\n");
     656                 printf("Adding mass to dry cell %d %f %f %f %f \n", k, zv[3*k], zv[3*k+1], zv[3*k+2], wc[k]);
    655657             }
    656              wc[k] = max(wc[k], bmin) ;
     658             wc[k] = max(wc[k], bmin);
    657659
    658660             // Set vertex values as well. Seems that this shouldn't be
  • trunk/anuga_work/development/gareth/tests/urban_flow/ideal_urban.py

    r8396 r8450  
    1313
    1414from balanced_dev import *
    15 #from balanced_dev import create_domain_from_regions as create_domain_from_regions
     15from balanced_dev import create_domain_from_regions as create_domain_from_regions
    1616#------------------------------------------------------------------------------
    1717# Useful parameters for controlling this case
     
    2626#bankwidth = 2. # Width of the bank regions -- note that these protrude into the channel
    2727man_n=0.01 # Manning's n
    28 l0 = 0.1 #0.2501 # Length scale associated with triangle side length in channel (min_triangle area = 0.5*l0^2)
     28l0 = 0.05 #0.2501 # Length scale associated with triangle side length in channel (min_triangle area = 0.5*l0^2)
    2929
    3030#assert chan_width < floodplain_width, \
  • trunk/anuga_work/development/gareth/tests/urban_flow/plotme.py

    r8384 r8450  
    77# Define gauge points, following Soares Frazao et al (2008)
    88offset=numpy.array([7.7,+1.8])
    9 G1=numpy.array([2.65,1.15]) + offset
    10 G2=numpy.array([2.65, -0.60]) + offset
    11 G3=numpy.array([4.00, 1.15]) + offset
    12 G4=numpy.array([4.00,-0.8]) + offset
    13 G5=numpy.array([5.20,0.3]) + offset
    14 G6=numpy.array([-1.87,1.10]) + offset
     9Gauge_xy=numpy.array([ [2.65,1.15] + offset,
     10                       [2.65, -0.60] + offset,
     11                       [4.00, 1.15] + offset,
     12                       [4.00,-0.8] + offset,
     13                       [5.20,0.3] + offset,
     14                       [-1.87,1.10] + offset ])
    1515
     16# Read model data
    1617p2=util.get_output('urban_flow0p05.sww',0.001)
    1718p=util.get_centroids(p2,velocity_extrapolation=True)
    1819
    19 G1ind = ((p.x-G1[0])**2 + (p.y-G1[1])**2).argmin()
    20 pyplot.plot(p.time, p.stage[:,G1ind], 'o-')
    21 pyplot.savefig('G1.png')
    2220
    23 pyplot.clf()
    24 G2ind = ((p.x-G2[0])**2 + (p.y-G2[1])**2).argmin()
    25 pyplot.plot(p.time, p.stage[:,G2ind], 'o-')
    26 pyplot.savefig('G2.png')
     21# Read experimental data
     22data_dir='/home/gareth/storage/anuga_clean/anuga_jan12/trunk/anuga_work/development/gareth/tests/urban_flow/experimental_data/Isolated_Obstacle/Isolated Obstacle/'
     23exp_data_h=data_dir+'building_gauges_h.txt'
     24exp_data_uv=data_dir+'/gareths_edits/data.txt'
     25
     26# Dir to save figures
     27fig_dir='./figures/'
     28h_gauges=numpy.genfromtxt(exp_data_h, skip_header=2)
     29uv_gauges=numpy.genfromtxt(exp_data_uv, skip_header=2)
     30
     31for i in range((Gauge_xy.shape[0])):
     32    Gauge=Gauge_xy[i,:]
     33    G_ind = ((p.x-Gauge[0])**2 + (p.y-Gauge[1])**2).argmin()
     34
     35    pyplot.plot(p.time, p.stage[:,G_ind], 'o-', label='model')
     36    pyplot.plot(h_gauges[:,0], h_gauges[:,(i+1)],'-', label='data')
     37    fig_filename=fig_dir+'G'+str(i+1)+'_stage.png'
     38    pyplot.savefig(fig_filename)
     39    pyplot.legend()
     40    pyplot.clf()
     41
     42    if(i<5):
     43        # Note -- Gauge 6 does not have velocity measurements
     44        pyplot.plot(p.time, p.xvel[:,G_ind], 'o-')
     45        pyplot.plot(uv_gauges[:,0], uv_gauges[:,2*i+1],'-')
     46        fig_filename=fig_dir+'G'+str(i+1)+'_xvel.png'
     47        pyplot.savefig(fig_filename)
     48        pyplot.clf()
    2749
    2850
    29 pyplot.clf()
    30 G3ind = ((p.x-G3[0])**2 + (p.y-G3[1])**2).argmin()
    31 pyplot.plot(p.time, p.stage[:,G3ind], 'o-')
    32 pyplot.savefig('G3.png')
    33 
    34 
    35 pyplot.clf()
    36 G4ind = ((p.x-G4[0])**2 + (p.y-G4[1])**2).argmin()
    37 pyplot.plot(p.time, p.stage[:,G4ind], 'o-')
    38 pyplot.savefig('G4.png')
    39 
    40 
    41 pyplot.clf()
    42 G5ind = ((p.x-G5[0])**2 + (p.y-G5[1])**2).argmin()
    43 pyplot.plot(p.time, p.stage[:,G5ind], 'o-')
    44 pyplot.savefig('G5.png')
    45 
    46 pyplot.clf()
    47 G6ind = ((p.x-G6[0])**2 + (p.y-G6[1])**2).argmin()
    48 pyplot.plot(p.time, p.stage[:,G6ind], 'o-')
    49 pyplot.savefig('G6.png')
     51        pyplot.plot(p.time, p.yvel[:,G_ind], 'o-')
     52        pyplot.plot(uv_gauges[:,0], uv_gauges[:,2*i+2],'-')
     53        fig_filename=fig_dir+'G'+str(i+1)+'_yvel.png'
     54        pyplot.savefig(fig_filename)
     55        pyplot.clf()
Note: See TracChangeset for help on using the changeset viewer.