Changeset 5832


Ignore:
Timestamp:
Oct 11, 2008, 12:38:18 PM (16 years ago)
Author:
steve
Message:
 
Location:
anuga_work/development/anuga_1d
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/development/anuga_1d/comp_flux_ext.c

    r5831 r5832  
    188188                                // Update timestep based on edge i and possibly neighbour n
    189189                                if (max_speed > epsilon) {
    190                         // Original CFL calculation
    191                                        
    192                                         timestep = min(timestep, 0.5*cfl*areas[k]/max_speed); //Here, CFL=1.0 is assumed. ?????????????????????????????????????????????
    193                                         if (n>=0) {
    194                                                 timestep = min(timestep, 0.5*cfl*areas[n]/max_speed); //Here, CFL=1.0 is assumed. ?????????????????????????????????????????????
    195                                         }
     190                                    // Original CFL calculation
     191                                   
     192                                    timestep = min(timestep, 0.5*cfl*areas[k]/max_speed);
     193                                    if (n>=0) {
     194                                        timestep = min(timestep, 0.5*cfl*areas[n]/max_speed);
     195                                    }
    196196                                }
    197197            } // End edge i (and neighbour n)
    … …  
    263263  // Call underlying flux computation routine and update
    264264  // the explicit update arrays
    265   timestep = _compute_fluxes_ext(
    266                  cfl,
     265  timestep = _compute_fluxes_ext(cfl,
    267266                                 timestep,
    268267                                 epsilon,
    … …  
    329328    g                 = get_python_double(domain,"g");
    330329    h0                = get_python_double(domain,"h0");
    331         cfl                = get_python_double(domain,"cfl");
     330    cfl               = get_python_double(domain,"CFL");
    332331 
    333332   
    … …  
    358357    // the explicit update arrays
    359358    timestep = _compute_fluxes_ext(
    360                        cfl,
    361                                    timestep,
    362                                    epsilon,
    363                                    g,
    364                                    h0,
    365                                    (long*) neighbours -> data,
    366                                    (long*) neighbour_vertices -> data,
    367                                    (double*) normals -> data,
    368                                    (double*) areas -> data,
    369                                    (double*) stage_vertex_values -> data,
    370                                    (double*) xmom_vertex_values -> data,
    371                                    (double*) bed_vertex_values -> data,
    372                                    (double*) stage_boundary_values -> data,
    373                                    (double*) xmom_boundary_values -> data,
    374                                    (double*) stage_explicit_update -> data,
    375                                    (double*) xmom_explicit_update -> data,
    376                                    number_of_elements,
    377                                    (double*) max_speed_array -> data);
     359        cfl,
     360        timestep,
     361        epsilon,
     362        g,
     363        h0,
     364        (long*) neighbours -> data,
     365        (long*) neighbour_vertices -> data,
     366        (double*) normals -> data,
     367        (double*) areas -> data,
     368        (double*) stage_vertex_values -> data,
     369        (double*) xmom_vertex_values -> data,
     370        (double*) bed_vertex_values -> data,
     371        (double*) stage_boundary_values -> data,
     372        (double*) xmom_boundary_values -> data,
     373        (double*) stage_explicit_update -> data,
     374        (double*) xmom_explicit_update -> data,
     375        number_of_elements,
     376        (double*) max_speed_array -> data);
    378377
    379378

  • anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c

    r5827 r5832  
    3131                int i;
    3232        double z_left, z_right;
    33                 double ql[2], qr[2], flux_left[2], flux_right[2];
     33                double flux_left[2], flux_right[2];
    3434                double z_star, w_left, h_left, h_left_star, uh_left, soundspeed_left, u_left;
    3535                double w_right, h_right, h_right_star, uh_right, soundspeed_right, u_right;                          //h_right,
    … …  
    136136               
    137137                double flux[2], ql[3], qr[3], edgeflux[2];
    138                 double zl, zr, max_speed, normal;
     138                double max_speed, normal;
    139139                int k, i, ki, n, m, nm=0;
    140140               

  • anuga_work/development/anuga_1d/domain.py

    r5742 r5832  
    77   Geoscience Australia
    88"""
     9
    910from generic_boundary_conditions import *
    10 #from coordinate_transforms.geo_reference import Geo_reference
     11
    1112
    1213class Domain:
    1314
    14     def __init__(self, coordinates, boundary = None,
    15                  conserved_quantities = None, other_quantities = None,
     15    def __init__(self,
     16                 coordinates,
     17                 boundary = None,
     18                 conserved_quantities = None,
     19                 evolved_quantities = None,
     20                 other_quantities = None,
    1621                 tagged_elements = None):
    1722        """
    … …  
    2328        from config import timestepping_method
    2429        from config import CFL
    25 
    26 
    2730       
    2831        #Store Points
    2932        self.coordinates = array(coordinates)
    3033
    31 # #        if geo_reference is None:
    32 # #            self. = Geo_reference() #Use defaults
    33 # #        else:
    34 # #            self.geo_reference = geo_reference
    3534
    3635        #Register number of Elements
    … …  
    3837
    3938        self.beta = 1.0
    40         self.limiter = "minmod_kurganov"
    41         self.CFL = CFL
     39        self.set_limiter("minmod_kurganov")
     40        self.set_CFL(CFL)
    4241        self.set_timestepping_method(timestepping_method)
    4342
    … …  
    124123            self.conserved_quantities = conserved_quantities
    125124
     125        if evolved_quantities is None:
     126            self.evolved_quantities = self.conserved_quantities
     127        else:
     128            self.evolved_quantities = evolved_quantities
     129
    126130        if other_quantities is None:
    127131            self.other_quantities = []
    … …  
    133137        self.quantities = {}
    134138
     139        #print self.conserved_quantities
     140        #print self.evolved_quantities
     141       
     142
    135143        #FIXME: remove later - maybe OK, though....
    136         for name in self.conserved_quantities:
     144        for name in self.evolved_quantities:
    137145            self.quantities[name] = Quantity(self)
    138146        for name in self.other_quantities:
    … …  
    160168        self.number_of_steps = 0
    161169        self.number_of_first_order_steps = 0
    162         self.CFL = CFL
    163170
    164171        #Model time
    … …  
    169176        #Checkpointing and storage
    170177        from config import default_datadir
    171         self.datadir = default_datadir
     178        self.set_datadir(default_datadir)
    172179        self.filename = 'domain'
    173180        self.checkpoint = False
    … …  
    490497
    491498        return q
     499
     500
     501    def get_evolved_quantities(self, vol_id, vertex=None):#, edge=None):
     502        """Get evolved quantities at volume vol_id
     503
     504        If vertex is specified use it as index for vertex values
     505        If edge is specified use it as index for edge values
     506        If neither are specified use centroid values
     507        If both are specified an exeception is raised
     508
     509        Return value: Vector of length == number_of_evolved quantities
     510
     511        """
     512
     513        from Numeric import zeros, Float
     514
     515        #if not (vertex is None):# or edge is None):
     516        #    msg = 'Values for both vertex and edge was specified.'
     517        #    msg += 'Only one (or none) is allowed.'
     518       #     raise msg
     519
     520        q = zeros( len(self.evolved_quantities), Float)
     521
     522        for i, name in enumerate(self.evolved_quantities):
     523            Q = self.quantities[name]
     524            if vertex is not None:
     525                q[i] = Q.vertex_values[vol_id, vertex]
     526            #elif edge is not None:
     527            #    q[i] = Q.edge_values[vol_id, edge]
     528            else:
     529                q[i] = Q.centroid_values[vol_id]
     530
     531        return q
    492532       
    493533
    … …  
    677717    def check_integrity(self):
    678718        #Mesh.check_integrity(self)
     719
     720        #print self.quantities
     721        #print self.conserved_quantities
    679722       
    680723        for quantity in self.conserved_quantities:
    681724            msg = 'Conserved quantities must be a subset of all quantities'
    682725            assert quantity in self.quantities, msg
     726
     727        for quantity in self.evolved_quantities:
     728            msg = 'Evolved quantities must be a subset of all quantities'
     729            assert quantity in self.quantities, msg           
    683730
    684731        # #assert hasattr(self, 'boundary_objects')
    … …  
    720767        self.datadir = name
    721768
    722 
     769    def set_CFL(self, cfl):
     770        if cfl > 1.0:
     771            print 'WARNING: Setting CFL condition to %f which is greater than 1' % cfl
     772        self.CFL = cfl
     773
     774    def get_CFL(self):
     775        return self.CFL
     776   
     777    def set_filename(self, name):
     778        self.filename = name
     779
     780    def get_filename(self):
     781        return self.filename
     782
     783    def get_limiter(self):
     784        return self.limiter
     785
     786    def set_limiter(self,limiter):
     787
     788        possible_limiters = \
     789        ['pyvolution', 'steve_minmod', 'minmod', 'minmod_kurganov', 'superbee', 'vanleer', 'vanalbada']
     790
     791        if limiter in possible_limiters:
     792            self.limiter = limiter
     793            return
     794
     795        msg = '%s is an incorrect limiter type.\n'% limiter
     796        msg += 'Possible types are: '+ ", ".join(["%s" % el for el in possible_limiters])
     797        raise Exception, msg
     798
     799       
    723800    #--------------------------
    724801    # Main components of evolve
    … …  
    915992        """
    916993
    917         # Save initial initial conserved quantities values
     994        # Save initial conserved quantities values
    918995        self.backup_conserved_quantities()           
    919996
    … …  
    13821459            q = B.evaluate(vol_id, vertex_id)
    13831460
    1384             for j, name in enumerate(self.conserved_quantities):
     1461            for j, name in enumerate(self.evolved_quantities):
    13851462                Q = self.quantities[name]
    13861463                Q.boundary_values[i] = q[j]

  • anuga_work/development/anuga_1d/dry_dam_sudi.py

    r5827 r5832  
    8181import time
    8282
    83 finaltime = 20.0
    84 yieldstep = 20.0
     83finaltime = 5.0
     84yieldstep = finaltime
    8585L = 2000.0     # Length of channel (m)
    86 number_of_cells = [810]#,200,500,1000,2000,5000,10000,20000]
     86number_of_cells = [1610]#,200,500,1000,2000,5000,10000,20000]
    8787h_error = zeros(len(number_of_cells),Float)
    8888uh_error = zeros(len(number_of_cells),Float)
    … …  
    101101    domain.set_boundary({'exterior': Reflective_boundary(domain)})
    102102    domain.order = 2
    103     domain.set_timestepping_method('rk2')
     103    domain.set_timestepping_method('rk3')
    104104    domain.cfl = 1.0
    105105    domain.limiter = "vanleer"

  • anuga_work/development/anuga_1d/shallow_water_domain.py

    r5741 r5832  
    5454        conserved_quantities = ['stage', 'xmomentum']
    5555        other_quantities = ['elevation', 'friction']
    56         Generic_Domain.__init__(self, coordinates, boundary,
    57                                 conserved_quantities, other_quantities,
    58                                 tagged_elements)
     56        Generic_Domain.__init__(self,
     57                                coordinates = coordinates,
     58                                boundary = boundary,
     59                                conserved_quantities = conserved_quantities,
     60                                other_quantities = other_quantities,
     61                                tagged_elements = tagged_elements)
    5962       
    6063        from config import minimum_allowed_height, g, h0
    … …  
    6871        #print "\nI have Removed forcing terms line 64 1dsw"
    6972
    70         #Realtime visualisation
    71         self.visualiser = None
    72         self.visualise  = False
    73         self.visualise_color_stage = False
    74         self.visualise_stage_range = 1.0
    75         self.visualise_timer = True
    76         self.visualise_range_z = None
    7773       
    7874        #Stored output
    … …  
    8278
    8379        #Evolve parametrs
    84         self.cfl = 1.0
     80        self.set_CFL(1.0)
    8581       
    8682        #Reduction operation for get_vertex_values
    … …  
    8985        #self.reduction = min  #Looks better near steep slopes
    9086
    91         self.quantities_to_be_stored = ['stage','xmomentum']
     87        self.set_quantities_to_be_stored(['stage','xmomentum'])
    9288
    9389        self.__doc__ = 'shallow_water_domain'
    9490
     91        self.check_integrity()
    9592
    9693    def set_quantities_to_be_stored(self, q):
    … …  
    125122       
    126123
    127     def initialise_visualiser(self,scale_z=1.0,rect=None):
    128         #Realtime visualisation
    129         if self.visualiser is None:
    130             from realtime_visualisation_new import Visualiser
    131             self.visualiser = Visualiser(self,scale_z,rect)
    132             self.visualiser.setup['elevation']=True
    133             self.visualiser.updating['stage']=True
    134         self.visualise = True
    135         if self.visualise_color_stage == True:
    136             self.visualiser.coloring['stage'] = True
    137             self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
    138         print 'initialise visualiser'
    139         print self.visualiser.setup
    140         print self.visualiser.updating
    141 
    142124    def check_integrity(self):
    143         Generic_Domain.check_integrity(self)
     125
    144126        #Check that we are solving the shallow water wave equation
    145127
    … …  
    149131        assert self.conserved_quantities[1] == 'xmomentum', msg
    150132
     133        msg = 'First evolved quantity must be "stage"'
     134        assert self.evolved_quantities[0] == 'stage', msg
     135        msg = 'Second evolved quantity must be "xmomentum"'
     136        assert self.evolved_quantities[1] == 'xmomentum', msg       
     137
     138        Generic_Domain.check_integrity(self)
     139       
    151140    def extrapolate_second_order_sw(self):
    152141        #Call correct module function
    … …  
    584573        import sys
    585574
     575        cfl = domain.get_CFL()
    586576       
    587577        timestep = float(sys.maxint)
    … …  
    596586        #print 'g=',g
    597587        #print 'The type of g is',type(g)
     588
     589        h0 = domain.h0
    598590       
    599591        neighbours = domain.neighbours
    … …  
    652644        from comp_flux_ext import compute_fluxes_ext
    653645       
    654         domain.flux_timestep = compute_fluxes_ext(timestep,
    655                                   epsilon,
    656                                   g,
    657                                   neighbours,
    658                                   neighbour_vertices,
    659                                   normals,
    660                                   areas,
    661                                   stage_edge_values,
    662                                   xmom_edge_values,
    663                                   bed_edge_values,
    664                                   stage_boundary_values,
    665                                   xmom_boundary_values,
    666                                   stage_explicit_update,
    667                                   xmom_explicit_update,
    668                                   number_of_elements,
    669                                   max_speed_array)
     646        domain.flux_timestep = compute_fluxes_ext(cfl,
     647                                                  timestep,
     648                                                  epsilon,
     649                                                  g,
     650                                                  h0,
     651                                                  neighbours,
     652                                                  neighbour_vertices,
     653                                                  normals,
     654                                                  areas,
     655                                                  stage_edge_values,
     656                                                  xmom_edge_values,
     657                                                  bed_edge_values,
     658                                                  stage_boundary_values,
     659                                                  xmom_boundary_values,
     660                                                  stage_explicit_update,
     661                                                  xmom_explicit_update,
     662                                                  number_of_elements,
     663                                                  max_speed_array)
    670664
    671665

  • anuga_work/development/anuga_1d/shallow_water_domain_suggestion1.py

    r5827 r5832  
    5353
    5454        conserved_quantities = ['stage', 'xmomentum']
    55         other_quantities = ['elevation', 'friction', 'height', 'velocity']
    56         Generic_Domain.__init__(self, coordinates, boundary,
    57                                 conserved_quantities, other_quantities,
    58                                 tagged_elements)
     55        evolved_quantities = ['stage', 'xmomentum', 'elevation', 'height', 'velocity']
     56        other_quantities = ['friction']
     57        Generic_Domain.__init__(self,
     58                                coordinates = coordinates,
     59                                boundary = boundary,
     60                                conserved_quantities = conserved_quantities,
     61                                evolved_quantities = evolved_quantities,
     62                                other_quantities = other_quantities,
     63                                tagged_elements = tagged_elements)
    5964       
    6065        from config import minimum_allowed_height, g, h0
    … …  
    6873        #print "\nI have Removed forcing terms line 64 1dsw"
    6974
    70         #Realtime visualisation
    71         self.visualiser = None
    72         self.visualise  = False
    73         self.visualise_color_stage = False
    74         self.visualise_stage_range = 1.0
    75         self.visualise_timer = True
    76         self.visualise_range_z = None
    7775       
    7876        #Stored output
    … …  
    8179        self.smooth = True
    8280
    83         #Evolve parametrs
    84         self.cfl = 1.0
    8581       
    8682        #Reduction operation for get_vertex_values
    … …  
    8985        #self.reduction = min  #Looks better near steep slopes
    9086
    91         self.quantities_to_be_stored = ['stage','xmomentum']
     87        self.set_quantities_to_be_stored(['stage','xmomentum'])
    9288
    9389        self.__doc__ = 'shallow_water_domain'
     90
     91        self.check_integrity()
    9492
    9593
    … …  
    125123       
    126124
    127     def initialise_visualiser(self,scale_z=1.0,rect=None):
    128         #Realtime visualisation
    129         if self.visualiser is None:
    130             from realtime_visualisation_new import Visualiser
    131             self.visualiser = Visualiser(self,scale_z,rect)
    132             self.visualiser.setup['elevation']=True
    133             self.visualiser.updating['stage']=True
    134         self.visualise = True
    135         if self.visualise_color_stage == True:
    136             self.visualiser.coloring['stage'] = True
    137             self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
    138         print 'initialise visualiser'
    139         print self.visualiser.setup
    140         print self.visualiser.updating
    141 
    142125    def check_integrity(self):
    143         Generic_Domain.check_integrity(self)
     126
    144127        #Check that we are solving the shallow water wave equation
    145128
    … …  
    148131        msg = 'Second conserved quantity must be "xmomentum"'
    149132        assert self.conserved_quantities[1] == 'xmomentum', msg
     133
     134        msg = 'First evolved quantity must be "stage"'
     135        assert self.evolved_quantities[0] == 'stage', msg
     136        msg = 'Second evolved quantity must be "xmomentum"'
     137        assert self.evolved_quantities[1] == 'xmomentum', msg
     138        msg = 'Third evolved quantity must be "elevation"'
     139        assert self.evolved_quantities[2] == 'elevation', msg
     140        msg = 'Fourth evolved quantity must be "height"'
     141        assert self.evolved_quantities[3] == 'height', msg
     142        msg = 'Fifth evolved quantity must be "velocity"'
     143        assert self.evolved_quantities[4] == 'velocity', msg
     144
     145        Generic_Domain.check_integrity(self)
    150146
    151147    def extrapolate_second_order_sw(self):
    … …  
    10801076
    10811077        #Handy shorthands
    1082         #self.stage   = domain.quantities['stage'].edge_values
    1083         #self.xmom    = domain.quantities['xmomentum'].edge_values
    1084         #self.ymom    = domain.quantities['ymomentum'].edge_values
    1085         self.normals = domain.normals
    1086         self.stage   = domain.quantities['stage'].vertex_values
    1087         self.xmom    = domain.quantities['xmomentum'].vertex_values
     1078        self.normals  = domain.normals
     1079        self.stage    = domain.quantities['stage'].vertex_values
     1080        self.xmom     = domain.quantities['xmomentum'].vertex_values
     1081        self.bed      = domain.quantities['elevation'].vertex_values
     1082        self.height   = domain.quantities['height'].vertex_values
     1083        self.velocity = domain.quantities['velocity'].vertex_values
    10881084
    10891085        from Numeric import zeros, Float
    10901086        #self.conserved_quantities = zeros(3, Float)
    1091         self.conserved_quantities = zeros(2, Float)
     1087        self.evolved_quantities = zeros(5, Float)
    10921088
    10931089    def __repr__(self):
    … …  
    11001096        """
    11011097
    1102         q = self.conserved_quantities
     1098        q = self.evolved_quantities
    11031099        q[0] = self.stage[vol_id, edge_id]
    1104         q[1] = self.xmom[vol_id, edge_id]
    1105         #q[2] = self.ymom[vol_id, edge_id]
    1106         #normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
    1107         #normal = self.normals[vol_id, 2*edge_id:2*edge_id+1]
    1108         normal = self.normals[vol_id,edge_id]
    1109 
    1110         #r = rotate(q, normal, direction = 1)
    1111         #r[1] = -r[1]
    1112         #q = rotate(r, normal, direction = -1)
    1113         r = q
    1114         r[1] = normal*r[1]
    1115         r[1] = -r[1]
    1116         r[1] = normal*r[1]
    1117         q = r
    1118         #For start interval there is no outward momentum so do not need to
    1119         #reverse direction in this case
     1100        q[1] = -self.xmom[vol_id, edge_id]
     1101        q[2] = self.bed[vol_id, edge_id]
     1102        q[3] = self.height[vol_id, edge_id]
     1103        q[1] = -self.velocity[vol_id, edge_id]
     1104
    11201105
    11211106        return q
    … …  
    11271112
    11281113
    1129     def __init__(self, conserved_quantities=None):
     1114    def __init__(self, evolved_quantities=None):
    11301115        Boundary.__init__(self)
    11311116
    1132         if conserved_quantities is None:
     1117        if evolved_quantities is None:
    11331118            msg = 'Must specify one value for each conserved quantity'
    11341119            raise msg
    11351120
    11361121        from Numeric import array, Float
    1137         self.conserved_quantities=array(conserved_quantities).astype(Float)
     1122        self.evolved_quantities=array(evolved_quantities).astype(Float)
    11381123
    11391124    def __repr__(self):
    … …  
    11411126
    11421127    def evaluate(self, vol_id=None, edge_id=None):
    1143         return self.conserved_quantities
     1128        return self.evolved_quantities
    11441129
    11451130

  • anuga_work/development/anuga_1d/test_shallow_water.py

    r5742 r5832  
    55
    66
    7 from shallow_water_domain_new import *
    8 from shallow_water_domain_new import flux_function as domain_flux_function
     7#from shallow_water_domain import *
     8#from shallow_water_domain import flux_function as domain_flux_function
     9
     10from shallow_water_domain_suggestion1 import *
     11from shallow_water_domain_suggestion1 import flux_function as domain_flux_function
    912
    1013from Numeric import allclose, array, ones, Float, maximum, zeros
    … …  
    3639
    3740        stage_ud, xmom_ud = local_compute_fluxes(domain)
    38        
     41
     42        domain.distribute_to_vertices_and_edges()
    3943        domain.compute_fluxes()
    4044
    41         #print domain.quantities['stage'].explicit_update
    42         #print domain.quantities['xmomentum'].explicit_update
    43         #print stage_ud
     45        print domain.quantities['stage'].explicit_update
     46        print domain.quantities['xmomentum'].explicit_update
     47        print stage_ud
     48        print xmom_ud
    4449
    4550        assert allclose( domain.quantities['stage'].explicit_update, stage_ud )
    … …  
    316321            try:
    317322                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
    318                 timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
     323                timestep = min(timestep, domain.CFL*0.5*domain.areas[k]/max_speed)
    319324            except ZeroDivisionError:
    320325                pass
Note: See TracChangeset for help on using the changeset viewer.