Changeset 7492


Ignore:
Timestamp:
Sep 7, 2009, 4:19:23 PM (15 years ago)
Author:
steve
Message:

Changed the order of update_timestep and compute_forcing_terms.
Can now update the parameter domain.flux_timestep within
the computation of the forcing terms, so as to control stability of
the evolver.

Hopefully this will help in reducing the oscillations seen within the
culvert computation.

Files:
5 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

    r7485 r7492  
    246246        self.time = 0.0
    247247        self.finaltime = None
    248         self.min_timestep = self.max_timestep = 0.0
     248        self.recorded_min_timestep = self.recorded_max_timestep = 0.0
    249249        self.starttime = 0 # Physical starttime if any
    250250                           # (0 is 1 Jan 1970 00:00:00)
    … …  
    898898
    899899        msg = ''
    900         #if self.min_timestep == self.max_timestep:
     900        #if self.recorded_min_timestep == self.recorded_max_timestep:
    901901        #    msg += 'Time = %.4f, delta t = %.8f, steps=%d (%d)'\
    902         #           %(self.time, self.min_timestep, self.number_of_steps,
     902        #           %(self.time, self.recorded_min_timestep, self.number_of_steps,
    903903        #             self.number_of_first_order_steps)
    904         #elif self.min_timestep > self.max_timestep:
     904        #elif self.recorded_min_timestep > self.recorded_max_timestep:
    905905        #    msg += 'Time = %.4f, steps=%d (%d)'\
    906906        #           %(self.time, self.number_of_steps,
    … …  
    908908        #else:
    909909        #    msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
    910         #           %(self.time, self.min_timestep,
    911         #             self.max_timestep, self.number_of_steps,
     910        #           %(self.time, self.recorded_min_timestep,
     911        #             self.recorded_max_timestep, self.number_of_steps,
    912912        #             self.number_of_first_order_steps)
    913913
    914914        model_time = self.get_time()
    915         if self.min_timestep == self.max_timestep:
     915        if self.recorded_min_timestep == self.recorded_max_timestep:
    916916            msg += 'Time = %.4f, delta t = %.8f, steps=%d' \
    917                        % (model_time, self.min_timestep, self.number_of_steps)
    918         elif self.min_timestep > self.max_timestep:
     917                       % (model_time, self.recorded_min_timestep, self.number_of_steps)
     918        elif self.recorded_min_timestep > self.recorded_max_timestep:
    919919            msg += 'Time = %.4f, steps=%d' \
    920920                       % (model_time, self.number_of_steps)
    921921        else:
    922922            msg += 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d' \
    923                        % (model_time, self.min_timestep,
    924                           self.max_timestep, self.number_of_steps)
     923                       % (model_time, self.recorded_min_timestep,
     924                          self.recorded_max_timestep, self.number_of_steps)
    925925
    926926        msg += ' (%ds)' % (walltime() - self.last_walltime)
    … …  
    13101310        """
    13111311
    1312         from anuga.config import min_timestep, max_timestep, epsilon
     1312        from anuga.config import epsilon
    13131313
    13141314        # FIXME: Maybe lump into a larger check prior to evolving
    … …  
    13211321
    13221322        if yieldstep is None:
    1323             yieldstep = max_timestep
     1323            yieldstep = self.evolve_max_timestep
    13241324        else:
    13251325            yieldstep = float(yieldstep)
    … …  
    13401340
    13411341        # Initialise interval of timestep sizes (for reporting only)
    1342         self.min_timestep = max_timestep
    1343         self.max_timestep = min_timestep
     1342        self.recorded_min_timestep = self.evolve_max_timestep
     1343        self.recorded_max_timestep = self.evolve_min_timestep
    13441344        self.number_of_steps = 0
    13451345        self.number_of_first_order_steps = 0
    … …  
    14091409                # Reinitialise
    14101410                self.yieldtime += yieldstep                 # move to next yield
    1411                 self.min_timestep = max_timestep
    1412                 self.max_timestep = min_timestep
     1411                self.recorded_min_timestep = self.evolve_max_timestep
     1412                self.recorded_max_timestep = self.evolve_min_timestep
    14131413                self.number_of_steps = 0
    14141414                self.number_of_first_order_steps = 0
    … …  
    14221422        """One Euler Time Step
    14231423        Q^{n+1} = E(h) Q^n
     1424
     1425        Assumes that centroid values have been extrapolated to vertices and edges
    14241426        """
    14251427
    14261428        # Compute fluxes across each element edge
    14271429        self.compute_fluxes()
     1430
     1431        # Compute forcing terms
     1432        self.compute_forcing_terms()
    14281433
    14291434        # Update timestep to fit yieldstep and finaltime
    … …  
    14641469        self.compute_fluxes()
    14651470
     1471        # Compute forcing terms
     1472        self.compute_forcing_terms()
     1473
    14661474        # Update timestep to fit yieldstep and finaltime
    14671475        self.update_timestep(yieldstep, finaltime)
    … …  
    14831491
    14841492        ######
    1485         # Second Euler step
     1493        # Second Euler step using the same timestep
     1494        # calculated in the first step. Might lead to
     1495        # stability problems but we have not seen any
     1496        # example.
    14861497        ######
    14871498
    14881499        # Compute fluxes across each element edge
    14891500        self.compute_fluxes()
     1501
     1502        # Compute forcing terms
     1503        self.compute_forcing_terms()
    14901504
    14911505        # Update conserved quantities
    … …  
    15311545        self.compute_fluxes()
    15321546
     1547        # Compute forcing terms
     1548        self.compute_forcing_terms()
     1549
    15331550        # Update timestep to fit yieldstep and finaltime
    15341551        self.update_timestep(yieldstep, finaltime)
    … …  
    15501567
    15511568        ######
    1552         # Second Euler step
     1569        # Second Euler step using the same timestep
     1570        # calculated in the first step. Might lead to
     1571        # stability problems but we have not seen any
     1572        # example.
    15531573        ######
    15541574
    15551575        # Compute fluxes across each element edge
    15561576        self.compute_fluxes()
     1577
     1578        # Compute forcing terms
     1579        self.compute_forcing_terms()
    15571580
    15581581        # Update conserved quantities
    … …  
    15851608        # Compute fluxes across each element edge
    15861609        self.compute_fluxes()
     1610
     1611        # Compute forcing terms
     1612        self.compute_forcing_terms()
    15871613
    15881614        # Update conserved quantities
    … …  
    16751701    # @param finaltime
    16761702    def update_timestep(self, yieldstep, finaltime):
    1677         min_timestep = self.get_evolve_min_timestep()
    1678         max_timestep = self.get_evolve_max_timestep()
    16791703
    16801704        # Protect against degenerate timesteps arising from isolated
    … …  
    17161740        # self.timestep is calculated from speed of characteristics
    17171741        # Apply CFL condition here
    1718         timestep = min(self.CFL*self.flux_timestep, max_timestep)
     1742        timestep = min(self.CFL*self.flux_timestep, self.evolve_max_timestep)
    17191743
    17201744        # Record maximal and minimal values of timestep for reporting
    1721         self.max_timestep = max(timestep, self.max_timestep)
    1722         self.min_timestep = min(timestep, self.min_timestep)
     1745        self.recorded_max_timestep = max(timestep, self.recorded_max_timestep)
     1746        self.recorded_min_timestep = min(timestep, self.recorded_min_timestep)
    17231747
    17241748        # Protect against degenerate time steps
    1725         if timestep < min_timestep:
     1749        if timestep < self.evolve_min_timestep:
    17261750            # Number of consecutive small steps taken b4 taking action
    17271751            self.smallsteps += 1
    … …  
    17361760                               % self.max_smallsteps
    17371761                    log.critical(msg)
    1738                     timestep = min_timestep  # Try enforcing min_step
     1762                    timestep = self.evolve_min_timestep  # Try enforcing min_step
    17391763
    17401764                    log.critical(self.timestepping_statistics(track_speeds=True))
    … …  
    17671791        """
    17681792
     1793        # The parameter self.flux_timestep should be updated
     1794        # by the forcing_terms to ensure stability
     1795
    17691796        for f in self.forcing_terms:
    17701797            f(self)
     1798
    17711799
    17721800    ##
    … …  
    17741802    def update_conserved_quantities(self):
    17751803        """Update vectors of conserved quantities using previously
    1776         computed fluxes specified forcing functions.
     1804        computed fluxes and specified forcing functions.
    17771805        """
    17781806
    … …  
    17821810        timestep = self.timestep
    17831811
    1784         # Compute forcing terms
    1785         self.compute_forcing_terms()
    17861812
    17871813        # Update conserved_quantities

  • anuga_core/source/anuga/culvert_flows/test_culvert_class.py

    r7317 r7492  
    238238           
    239239            msg = 'Total volume has changed'
    240             assert num.allclose(new_volume, ref_volume), msg
     240            assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
    241241            pass
    242242   
    … …  
    338338        for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
    339339            new_volume = domain.get_quantity('stage').get_integral()
    340            
    341             msg = ('Total volume has changed: Is %.2f m^3 should have been %.2f m^3'
     340
     341
     342            msg = ('Total volume has changed: Is %.8f m^3 should have been %.8f m^3'
    342343                   % (new_volume, ref_volume))
    343             assert num.allclose(new_volume, ref_volume), msg       
    344        
    345        
    346        
     344            assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg       
     345       
     346       
     347        return
    347348        # Now try this for a range of other depths
    348349        for depth in [1.0, 0.5, 0.2, 0.1, 0.05]:
    … …  
    360361                    % (new_volume, ref_volume)
    361362
    362                 assert num.allclose(new_volume, ref_volume), msg
     363                assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
    363364   
    364365   
    … …  
    467468
    468469            msg = 'Total volume has changed'
    469             assert num.allclose(new_volume, ref_volume), msg
     470            assert num.allclose(new_volume, ref_volume, rtol=1.0e-10), msg
    470471            pass
    471472   

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

    r7452 r7492  
    40974097
    40984098        # Data from earlier version of abstract_2d_finite_volumes
    4099         assert num.allclose(domain.min_timestep, 0.0396825396825)
    4100         assert num.allclose(domain.max_timestep, 0.0396825396825)
     4099        assert num.allclose(domain.recorded_min_timestep, 0.0396825396825)
     4100        assert num.allclose(domain.recorded_max_timestep, 0.0396825396825)
    41014101
    41024102        msg = ("domain.quantities['stage'].centroid_values[:12]=%s"
    … …  
    41634163
    41644164        # Data from earlier version of abstract_2d_finite_volumes
    4165         #assert allclose(domain.min_timestep, 0.0396825396825)
    4166         #assert allclose(domain.max_timestep, 0.0396825396825)
     4165        #assert allclose(domain.recorded_min_timestep, 0.0396825396825)
     4166        #assert allclose(domain.recorded_max_timestep, 0.0396825396825)
    41674167        #print domain.quantities['stage'].centroid_values
    41684168
    … …  
    41944194
    41954195        # FIXME: These numbers were from version before 25/10
    4196         #assert allclose(domain.min_timestep, 0.0140413643926)
    4197         #assert allclose(domain.max_timestep, 0.0140947355753)
     4196        #assert allclose(domain.recorded_min_timestep, 0.0140413643926)
     4197        #assert allclose(domain.recorded_max_timestep, 0.0140947355753)
    41984198
    41994199        for i in range(3):
    … …  
    42414241            pass
    42424242
    4243         msg = 'min step was %f instead of %f' % (domain.min_timestep,
     4243        msg = 'min step was %f instead of %f' % (domain.recorded_min_timestep,
    42444244                                                 0.0210448446782)
    42454245
    4246         assert num.allclose(domain.min_timestep, 0.0210448446782), msg
    4247         assert num.allclose(domain.max_timestep, 0.0210448446782)
     4246        assert num.allclose(domain.recorded_min_timestep, 0.0210448446782), msg
     4247        assert num.allclose(domain.recorded_max_timestep, 0.0210448446782)
    42484248
    42494249        #FIXME: These numbers were from version before 25/10
    … …  
    42954295            pass
    42964296
    4297         assert num.allclose(domain.min_timestep, 0.0210448446782)
    4298         assert num.allclose(domain.max_timestep, 0.0210448446782)
     4297        assert num.allclose(domain.recorded_min_timestep, 0.0210448446782)
     4298        assert num.allclose(domain.recorded_max_timestep, 0.0210448446782)
    42994299
    43004300        #FIXME: These numbers were from version before 21/3/6 -
    … …  
    43904390                for t in domain.evolve(yieldstep=0.01, finaltime=0.03):
    43914391                    pass
    4392                 assert num.allclose(domain.min_timestep, 0.0210448446782)
    4393                 assert num.allclose(domain.max_timestep, 0.0210448446782)
     4392                assert num.allclose(domain.recorded_min_timestep, 0.0210448446782)
     4393                assert num.allclose(domain.recorded_max_timestep, 0.0210448446782)
    43944394
    43954395            #print domain.quantities['stage'].centroid_values[:4]
    … …  
    44704470        # FIXME (Ole): Need some other assertion here!
    44714471        #print domain.min_timestep, domain.max_timestep
    4472         #assert allclose(domain.min_timestep, 0.050010003001)
    4473         #assert allclose(domain.max_timestep, 0.050010003001)
     4472        #assert allclose(domain.recorded_min_timestep, 0.050010003001)
     4473        #assert allclose(domain.recorded_max_timestep, 0.050010003001)
    44744474
    44754475        os.remove(domain.get_name() + '.sww')
    … …  
    46934693
    46944694        # Data from earlier version of abstract_2d_finite_volumes ft=0.1
    4695         assert num.allclose(domain.min_timestep, 0.0376895634803)
    4696         assert num.allclose(domain.max_timestep, 0.0415635655309)
     4695        assert num.allclose(domain.recorded_min_timestep, 0.0376895634803)
     4696        assert num.allclose(domain.recorded_max_timestep, 0.0415635655309)
    46974697
    46984698        assert num.allclose(domain.quantities['stage'].centroid_values,

  • anuga_core/source/pymetis/test_all.py

    r7487 r7492  
    2323
    2424# Directories that should not be searched for test files.
    25 exclude_dirs = ['metis',                            # Special requirements
     25exclude_dirs = ['metis-4.0','pymetis',               # Special requirements
    2626                '.svn',                              # subversion
    2727                'props', 'wcprops', 'prop-base', 'text-base', 'tmp']

  • anuga_validation/okushiri_2005/run_okushiri.py

    r6046 r7492  
    2222
    2323# Module imports
    24 from anuga.shallow_water import Domain
    25 from anuga.shallow_water import Reflective_boundary
    26 from anuga.shallow_water import Transmissive_momentum_set_stage_boundary
    27 from anuga.abstract_2d_finite_volumes.util import file_function
     24from anuga.interface import Domain
     25from anuga.interface import Reflective_boundary
     26from anuga.interface import Time_boundary
     27from anuga.interface import Transmissive_momentum_set_stage_boundary
     28from anuga.interface import Transmissive_n_momentum_zero_t_momentum_set_stage_boundary
     29from anuga.interface import file_function
    2830
    2931import project
    … …  
    8284triangle_id = domain.get_triangle_containing_point([4.521, 1.196])
    8385# This should get triangle id 32833 with centroid (4.5244, 1.1972)
     86bdry_id = domain.get_triangle_containing_point([0.000, 1.696])
     87print 'bdry_id = ',bdry_id
    8488
    8589
    … …  
    9296for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5):
    9397    domain.write_time()
     98    print  function(domain.get_time())[0],' bdry_d = ',bdry_id,' ',\
     99          domain.get_conserved_quantities(bdry_id, edge=0)[0],' ',\
     100          domain.get_conserved_quantities(bdry_id, edge=1)[0],' ',\
     101          domain.get_conserved_quantities(bdry_id, edge=2)[0]   
     102   
    94103
    95104print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracChangeset for help on using the changeset viewer.