Python Domain.evolve方法代码示例

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import evolve [as 别名]
    def concept_ungenerateIII(self):
        from anuga import Domain, Reflective_boundary, \
                            Dirichlet_boundary
        from anuga.pmesh.mesh_interface import create_mesh_from_regions

        # These are the absolute values
        polygon = [[0,0], [100,0], [100,100], [0,100]]

        boundary_tags = {'wall': [0,1,3], 'wave': [2]}
        inner1_polygon = [[10,10], [20,10], [20,20], [10,20]]
        inner2_polygon = [[30,30], [40,30], [40,40], [30,40]]

        max_area = 1
        interior_regions = [(inner1_polygon, 5), (inner2_polygon, 10)]
        m = create_mesh_from_regions(polygon,
                                     boundary_tags,
                                     max_area,
                                     interior_regions=interior_regions)

        fileName = tempfile.mktemp('.txt')
        file = open(fileName, 'w')
        file.write('         1       ??      ??\n\
       90.0       90.0\n\
       81.0       90.0\n\
       81.0       81.0\n\
       90.0       81.0\n\
       90.0       90.0\n\
END\n\
         2      ?? ??\n\
       10.0       80.0\n\
       10.0       90.0\n\
       20.0       90.0\n\
       10.0       80.0\n\
END\n\
END\n')
        file.close()

        m.import_ungenerate_file(fileName)
        os.remove(fileName)
        m.generate_mesh(maximum_triangle_area=max_area, verbose=False)
        mesh_filename = 'mesh.tsh'
        m.export_mesh_file(mesh_filename)

        domain = Domain(mesh_filename, use_cache=False)

        Br = Reflective_boundary(domain)
        Bd = Dirichlet_boundary([3, 0, 0])
        domain.set_boundary({'wall': Br, 'wave': Bd})
        yieldstep = 0.1
        finaltime = 10
        for t in domain.evolve(yieldstep, finaltime):
            domain.write_time()

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import evolve [as 别名]
    def test_runup_sinusoid(self):
        """ Run a version of the validation test runup_sinusoid
        to ensure limiting solution has small velocity
        """

        points, vertices, boundary = anuga.rectangular_cross(20,20, len1=1., len2=1.)


        domain=Domain(points,vertices,boundary)    # Create Domain
        domain.set_flow_algorithm('DE0')
        
        domain.set_name('runup_sinusoid_v2')                         # Output to file runup.sww
        domain.set_datadir('.')                          # Use current folder
        domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
        #domain.set_store_vertices_uniquely(True)
        #------------------
        # Define topography
        #------------------
        scale_me=1.0

        def topography(x,y):
            return (-x/2.0 +0.05*num.sin((x+y)*50.0))*scale_me

        def stagefun(x,y):
            stge=-0.2*scale_me #+0.01*(x>0.9)
            return stge

        domain.set_quantity('elevation',topography)     # Use function for elevation
        domain.get_quantity('elevation').smooth_vertex_values()
        domain.set_quantity('friction',0.03)            # Constant friction


        domain.set_quantity('stage', stagefun)             # Constant negative initial stage
        domain.get_quantity('stage').smooth_vertex_values()


        #--------------------------
        # Setup boundary conditions
        #--------------------------
        Br=anuga.Reflective_boundary(domain)                 # Solid reflective wall
        Bd=anuga.Dirichlet_boundary([-0.1*scale_me,0.,0.])   # Constant boundary values -- not used in this example

        #----------------------------------------------
        # Associate boundary tags with boundary objects
        #----------------------------------------------
        domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom':Br})

        #------------------------------
        #Evolve the system through time
        #------------------------------

        for t in domain.evolve(yieldstep=7.0,finaltime=7.0):
            #print domain.timestepping_statistics()
            xx = domain.quantities['xmomentum'].centroid_values
            yy = domain.quantities['ymomentum'].centroid_values
            dd = domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values
            #dd_raw=1.0*dd
            dd = (dd)*(dd>1.0e-03)+1.0e-03
            vv = ( (xx/dd)**2 + (yy/dd)**2)**0.5
            vv = vv*(dd>1.0e-03)
            #print 'Peak velocity is: ', vv.max(), vv.argmax()
            #print 'Volume is', sum(dd_raw*domain.areas)


        #print vv.max()

        assert num.all(vv<1.01e-01)

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import evolve [as 别名]
# Associate boundary tags with boundary objects
domain.set_boundary({'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})




#------------------------------------------------------------------------------
# Produce a documentation of parameters
#------------------------------------------------------------------------------
if myid == 0:
    parameter_file=open('parameters.tex', 'w')
    parameter_file.write('\\begin{verbatim}\n')
    from pprint import pprint
    pprint(domain.get_algorithm_parameters(),parameter_file,indent=4)
    parameter_file.write('\\end{verbatim}\n')
    parameter_file.close()

#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------
for t in domain.evolve(yieldstep = 1.0, finaltime = 20.):
    #print domain.timestepping_statistics(track_speeds=True)
    if myid == 0 and verbose: print domain.timestepping_statistics()


domain.sww_merge(delete_old=True)


finalize()

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import evolve [as 别名]
domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})

#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------
polygon1 = [ [10.0, 0.0], [11.0, 0.0], [11.0, 5.0], [10.0, 5.0] ]
polygon2 = [ [12.0, 2.0], [13.0, 2.0], [13.0, 3.0], [12.0, 3.0] ]

from anuga.operators.rate_operators import Rate_operator

op1 = Rate_operator(domain, rate=lambda t: 10.0 if (t>=0.0) else 0.0, polygon=polygon2)
op2 = Rate_operator(domain, rate=lambda t: 10.0 if (t>=0.0) else 0.0, radius=0.5, center=(10.0, 3.0))


domain.set_starttime(-0.1)
for t in domain.evolve(yieldstep=0.01, finaltime=0.0):
    domain.print_timestepping_statistics()
    domain.print_operator_timestepping_statistics()

    stage = domain.get_quantity('stage')
    elev  = domain.get_quantity('elevation')
    height = stage - elev

    print 'integral = ', height.get_integral()


for t in domain.evolve(yieldstep=0.1, duration=5.0):

    domain.print_timestepping_statistics()
    domain.print_operator_timestepping_statistics()

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import evolve [as 别名]
Bw5 = anuga.Flather_external_stage_zero_velocity_boundary(domain, zero_fun)

# Associate boundary tags with boundary objects
domain.set_boundary({'left': Bw2, 'right': Bw5, 'top': Br, 'bottom': Br})

#------------------------------------------------------------------------------
# Produce a documentation of parameters
#------------------------------------------------------------------------------
if myid == 0:
    parameter_file=open('parameters.tex', 'w')
    parameter_file.write('\\begin{verbatim}\n')
    from pprint import pprint
    pprint(domain.get_algorithm_parameters(),parameter_file,indent=4)
    parameter_file.write('\\end{verbatim}\n')
    parameter_file.close()

#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------

for t in domain.evolve(yieldstep = 10, finaltime = 2e4):
    domain.write_time()
    if interactive_visualisation:
        vis.update()

domain.sww_merge(delete_old=True)

finalize()

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import evolve [as 别名]
    return [wO,  pO, 0.0]        # dimensional

Br = anuga.Reflective_boundary(domain)      # Solid reflective wall
Bt = anuga.Transmissive_boundary(domain)    # Continue all values on boundary 
#Bd = anuga.Dirichlet_boundary([1,0.,0.])    # Constant boundary values
BTime = anuga.Time_boundary(domain,f_CG)
# Associate boundary tags with boundary objects
domain.set_boundary({'left': BTime, 'right': Bt, 'top': Br, 'bottom': Br})

#------------------------------------------------------------------------------
# Produce a documentation of parameters
#------------------------------------------------------------------------------
if myid == 0:
    parameter_file=open('parameters.tex', 'w')
    parameter_file.write('\\begin{verbatim}\n')
    from pprint import pprint
    pprint(domain.get_algorithm_parameters(),parameter_file,indent=4)
    parameter_file.write('\\end{verbatim}\n')
    parameter_file.close()

#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------
for t in domain.evolve(yieldstep = Tp/48., finaltime = 7000.):
    if myid == 0 and verbose: print domain.timestepping_statistics()


domain.sww_merge(delete_old=True)

finalize()

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import evolve [as 别名]
factor = 1e-3
op3 = Rate_operator(domain, rate = rain, factor=factor)
Q3 = numpy.sum(op3.get_spatial_rate()*domain.areas)*factor


#op3()
#domain.fractional_step_operators.remove(op3)


#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------
accum = 0.0
yieldstep = 0.1
finaltime = 5.0
for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
    domain.print_timestepping_statistics()
    domain.print_operator_timestepping_statistics()

    stage = domain.get_quantity('stage')
    elev  = domain.get_quantity('elevation')
    height = stage - elev

    print indent + 'Integral = ', height.get_integral()
    print indent + 'Exact accumultion = ', accum
    
    dd = max(min(yieldstep,4.0-t),0.0)
    accum += (Q1+Q2)*yieldstep + dd*Q3

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import evolve [as 别名]
##vis.start()
#===============================================================================


if myid == 0:
    #------------------------------------------------------------------------------
    # Produce a documentation of parameters
    #------------------------------------------------------------------------------
    parameter_file=open('parameters.tex', 'w')
    parameter_file.write('\\begin{verbatim}\n')
    from pprint import pprint
    pprint(domain.get_algorithm_parameters(),parameter_file,indent=4)
    parameter_file.write('\\end{verbatim}\n')
    parameter_file.close()

#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------
for t in domain.evolve(yieldstep = 1000., finaltime = 30000.):
    #print domain.timestepping_statistics(track_speeds=True)
    if myid == 0 and verbose:
        print domain.timestepping_statistics()
    #vis.update()


domain.sww_merge(delete_old=True)

finalize()

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import evolve [as 别名]
    def concept_ungenerateII(self):
        from anuga import Domain, Reflective_boundary, Dirichlet_boundary

        x=0
        y=0
        mesh_geo = geo_reference=Geo_reference(56, x, y)

        # These are the absolute values
        polygon_absolute = [[0,0], [100,0], [100,100], [0,100]]
        x_p = -10
        y_p = -40
        geo_ref_poly = Geo_reference(56, x_p, y_p)
        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)

        boundary_tags = {'wall': [0,1,3], 'wave': [2]}

        inner1_polygon_absolute = [[10,10], [20,10], [20,20], [10,20]]
        inner1_polygon = geo_ref_poly.\
                            change_points_geo_ref(inner1_polygon_absolute)

        inner2_polygon_absolute = [[30,30], [40,30], [40,40], [30,40]]
        inner2_polygon = geo_ref_poly.\
                            change_points_geo_ref(inner2_polygon_absolute)

        max_area = 1
        interior_regions = [(inner1_polygon, 5), (inner2_polygon, 10)]
        m = create_mesh_from_regions(polygon,
                                     boundary_tags,
                                     max_area,
                                     interior_regions=interior_regions,
                                     poly_geo_reference=geo_ref_poly,
                                     mesh_geo_reference=mesh_geo)

        m.export_mesh_file('a_test_mesh_iknterface.tsh')

        fileName = tempfile.mktemp('.txt')
        file = open(fileName, 'w')
        file.write('         1       ??      ??\n\
       90.0       90.0\n\
       81.0       90.0\n\
       81.0       81.0\n\
       90.0       81.0\n\
       90.0       90.0\n\
END\n\
         2      ?? ??\n\
       10.0       80.0\n\
       10.0       90.0\n\
       20.0       90.0\n\
       10.0       80.0\n\
END\n\
END\n')
        file.close()

        m.import_ungenerate_file(fileName)      #, tag='wall')
        os.remove(fileName)
        m.generate_mesh(maximum_triangle_area=max_area, verbose=False)
        mesh_filename = 'bento_b.tsh'
        m.export_mesh_file(mesh_filename)

        domain = Domain(mesh_filename, use_cache = False)

        Br = Reflective_boundary(domain)
        Bd = Dirichlet_boundary([3, 0, 0])
        domain.set_boundary({'wall': Br, 'wave': Bd})
        yieldstep = 0.1
        finaltime = 10
        for t in domain.evolve(yieldstep, finaltime):
            domain.write_time()

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import evolve [as 别名]
# Setup boundary conditions
#------------------------------------------------------------------------------
Bt = anuga.Transmissive_boundary(domain)
BdIN = anuga.Dirichlet_boundary([dana, fluxin, 0.0])
#BdOUT = anuga.Dirichlet_boundary([dana-10., fluxin, 0.0])
Br = anuga.Reflective_boundary(domain) # Solid reflective wall
domain.set_boundary({'left': BdIN, 'right': Bt, 'top': Br, 'bottom': Br})


#------------------------------------------------------------------------------
# Produce a documentation of parameters
#------------------------------------------------------------------------------
if myid == 0:
	parameter_file=open('parameters.tex', 'w')
	parameter_file.write('\\begin{verbatim}\n')
	from pprint import pprint
	pprint(domain.get_algorithm_parameters(),parameter_file,indent=4)
	parameter_file.write('\\end{verbatim}\n')
	parameter_file.close()

#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------
for t in domain.evolve(yieldstep=2.0, finaltime=200.0):
	if myid ==0 and verbose: print domain.timestepping_statistics()


domain.sww_merge(delete_old=True)

finalize()

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import evolve [as 别名]
#------------------------------------------------------------------------------
# Setup operators
#------------------------------------------------------------------------------

from anuga.operators.vegetation_operator import Vegetation_operator

veg_op = Vegetation_operator(domain)

Quantity(domain, name='veg_diameter', register=True)
domain.set_quantity('veg_diameter', 0.00064)

Quantity(domain, name='veg_spacing', register=True)
domain.set_quantity('veg_spacing', 0.15)

#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------
for t in domain.evolve(yieldstep = 1, finaltime = 30.0):
    domain.print_timestepping_statistics()

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import evolve [as 别名]
# Setup boundary conditions
#------------------------------------------------------------------------------
from math import sin, pi, exp
Br = anuga.Reflective_boundary(domain)      # Solid reflective wall
Bt = anuga.Transmissive_boundary(domain)    # Continue all values on boundary 
Bd = anuga.Dirichlet_boundary([1,0.,0.]) # Constant boundary values

# Associate boundary tags with boundary objects
domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})


#-------------------------------------------------------------------------
# Produce a documentation of parameters
#-------------------------------------------------------------------------
from anuga.validation_utilities import save_parameters_tex
save_parameters_tex(domain)


#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------
for t in domain.evolve(yieldstep = 0.05, finaltime = 3.0):
    if myid == 0 and verbose:
        print domain.timestepping_statistics()


domain.sww_merge(delete_old=True)

finalize()

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import evolve [as 别名]

#.........这里部分代码省略.........
            urs2sts(base_name,
                    basename_out=sts_file,
                    ordering_filename=order_file,
                    mean_stage=tide,
                    verbose=verbose)
            self.delete_mux(files)

            assert(os.access(sts_file+'.sts', os.F_OK))

            os.remove(order_file)

        barrier()
        boundary_polygon = create_sts_boundary(sts_file)

        # Append the remaining part of the boundary polygon to be defined by
        # the user
        bounding_polygon_utm=[]
        for point in bounding_polygon:
            zone,easting,northing=redfearn(point[0],point[1])
            bounding_polygon_utm.append([easting,northing])

        boundary_polygon.append(bounding_polygon_utm[3])
        boundary_polygon.append(bounding_polygon_utm[4])

        assert num.allclose(bounding_polygon_utm,boundary_polygon)


        extent_res=1000000
        meshname = 'urs_test_mesh' + '.tsh'
        interior_regions=None
        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
        
        # have to change boundary tags from last example because now bounding
        # polygon starts in different place.
        if myid==0:
            create_mesh_from_regions(boundary_polygon,
                                     boundary_tags=boundary_tags,
                                     maximum_triangle_area=extent_res,
                                     filename=meshname,
                                     interior_regions=interior_regions,
                                     verbose=verbose)

        barrier()
        
        domain_fbound = Domain(meshname)
        domain_fbound.set_quantities_to_be_stored(None)
        domain_fbound.set_quantity('stage', tide)
        if verbose: print "Creating file boundary condition"
        Bf = File_boundary(sts_file+'.sts',
                           domain_fbound,
                           boundary_polygon=boundary_polygon)
        Br = Reflective_boundary(domain_fbound)

        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})

        temp_fbound=num.zeros(int(finaltime/yieldstep)+1,num.float)
        if verbose: print "Evolving domain with file boundary condition"
        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,
                                                   finaltime=finaltime, 
                                                   skip_initial_step = False)):
            temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2]
            if verbose: domain_fbound.write_time()
            
        
        domain_drchlt = Domain(meshname)
        domain_drchlt.set_quantities_to_be_stored(None)
        domain_drchlt.set_starttime(time_step)
        domain_drchlt.set_quantity('stage', tide)
        Br = Reflective_boundary(domain_drchlt)
        #Bd = Dirichlet_boundary([2.0+tide,220+10*tide,-220-10*tide])
        Bd = Time_boundary(domain=domain_drchlt, f=lambda t: [2.0+t/finaltime+tide,220.+10.*tide+10.*t/finaltime,-220.-10.*tide-10.*t/finaltime])
        #Bd = Time_boundary(domain=domain_drchlt,f=lambda t: [2.0+num.sin(t)+tide,10.*(2+20.+num.sin(t)+tide),-10.*(2+20.+num.sin(t)+tide)])
        domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br})
        temp_drchlt=num.zeros(int(finaltime/yieldstep)+1,num.float)
        
        for i, t in enumerate(domain_drchlt.evolve(yieldstep=yieldstep,
                                                   finaltime=finaltime, 
                                                   skip_initial_step = False)):
            temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2]
            #domain_drchlt.write_time()
        
        #print domain_fbound.quantities['stage'].vertex_values
        #print domain_drchlt.quantities['stage'].vertex_values
                    
        assert num.allclose(temp_fbound,temp_drchlt),temp_fbound-temp_drchlt

        
        assert num.allclose(domain_fbound.quantities['stage'].vertex_values,
                            domain_drchlt.quantities['stage'].vertex_values)
                        
        assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values,
                            domain_drchlt.quantities['xmomentum'].vertex_values)                        
                        
        assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values,
                            domain_drchlt.quantities['ymomentum'].vertex_values)
        
        if not sys.platform == 'win32':
            if myid==0: os.remove(sts_file+'.sts')
        
        if myid==0: os.remove(meshname)

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import evolve [as 别名]

#.........这里部分代码省略.........
        #    P.plot(interpolation_points[:,0],interpolation_points[:,1],'ko')
        #    P.show()

        fbound_gauge_values = []
        fbound_proc_tri_ids = []
        for i, point in enumerate(interpolation_points):
            fbound_gauge_values.append([]) # Empty list for timeseries

            try:
                k = domain_fbound.get_triangle_containing_point(point)
                if domain_fbound.tri_full_flag[k] == 1:
                    fbound_proc_tri_ids.append(k)
                else:
                    fbound_proc_tri_ids.append(-1)            
            except:
                fbound_proc_tri_ids.append(-2)


        if verbose: print 'P%d has points = %s' %(myid, fbound_proc_tri_ids)

        #------------------------------------------------------------
        # Set boundary conditions
        #------------------------------------------------------------
        Bf = File_boundary(sts_file+'.sts',
                           domain_fbound,
                           boundary_polygon=boundary_polygon)
        Br = Reflective_boundary(domain_fbound)
    
        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})

        #------------------------------------------------------------
        # Evolve the domain on each processor
        #------------------------------------------------------------  
        for i, t in enumerate(domain_fbound.evolve(yieldstep=yieldstep,
                                                   finaltime=finaltime, 
                                                   skip_initial_step = False)):

            stage = domain_fbound.get_quantity('stage')
            for i in range(4):
                if fbound_proc_tri_ids[i] > -1:
                    fbound_gauge_values[i].append(stage.centroid_values[fbound_proc_tri_ids[i]])
        
        #------------------------------------------------------------
        # Create domain to be run sequntially on each processor
        #------------------------------------------------------------
        domain_drchlt = Domain(meshname)
        domain_drchlt.set_quantities_to_be_stored(None)
        domain_drchlt.set_starttime(time_step)
        domain_drchlt.set_quantity('stage', tide)
        Br = Reflective_boundary(domain_drchlt)
        #Bd = Dirichlet_boundary([2.0+tide,220+10*tide,-220-10*tide])
        Bd = Time_boundary(domain=domain_drchlt, function=lambda t: [2.0+t/finaltime+tide,220.+10.*tide+10.*t/finaltime,-220.-10.*tide-10.*t/finaltime])
        #Bd = Time_boundary(domain=domain_drchlt,function=lambda t: [2.0+num.sin(t)+tide,10.*(2+20.+num.sin(t)+tide),-10.*(2+20.+num.sin(t)+tide)])
        domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br})
       
        drchlt_gauge_values = []
        drchlt_proc_tri_ids = []
        for i, point in enumerate(interpolation_points):
            drchlt_gauge_values.append([]) # Empty list for timeseries

            try:
                k = domain_drchlt.get_triangle_containing_point(point)
                if domain_drchlt.tri_full_flag[k] == 1:
                    drchlt_proc_tri_ids.append(k)
                else:
                    drchlt_proc_tri_ids.append(-1)            

# 需要导入模块: from anuga import Domain [as 别名]
# 或者: from anuga.Domain import evolve [as 别名]
vtk_visualiser = True
if vtk_visualiser:
    from anuga.visualiser import RealtimeVisualiser
    vis = RealtimeVisualiser(domain)
    vis.render_quantity_height("height",dynamic=True)
    #vis.render_quantity_height("stage", zScale =1.0, dynamic=True)
    vis.colour_height_quantity('stage', (0.0, 0.0, 1.0))
    vis.start()
#===============================================================================


#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------
i = 0
for t in domain.evolve(yieldstep = 0.01, finaltime = 0.25):
    #print domain.timestepping_statistics(track_speeds=True)
    print domain.timestepping_statistics()
    if vtk_visualiser: 
        vis.update()
        fileName = 'stage_%03d' % i + '.vtk'
        i = i+1
        vis.store_height_quantity('height',fileName)

#test against know data


    
if vtk_visualiser: vis.evolveFinished()