Changeset 5084 for anuga_work/debug/bad_time_step

Timestamp:

Feb 27, 2008, 10:02:55 AM (17 years ago)

Author:

ole

Message:

Simplified debug of transmissive_momentum_set_stage

Location:

anuga_work/debug/bad_time_step

Files:

• boundary.tms (moved) (moved from anuga_work/debug/bad_time_step/boundary.tsm)
• boundary_mellow.tms (moved) (moved from anuga_work/debug/bad_time_step/boundary_mellow.tsm)
• create_mesh.py (modified) (1 diff)
• project.py (modified) (1 diff)
• run_dam.py (modified) (3 diffs)

anuga_work/debug/bad_time_step/create_mesh.py

@@ -14 +14 @@
     
     """
+    
     #Basic geometry
-    
     xright  = 19.0
     ybottom = 0

anuga_work/debug/bad_time_step/project.py

@@ -12 +12 @@
 from anuga.utilities.system_tools import get_user_name
 
-boundary_file = 'boundary_mellow.tsm'
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
+     import File_boundary_time
+boundary_file = 'boundary.tms'
+#boundary_file = 'boundary_mellow.tms'
 
 class Project:    

anuga_work/debug/bad_time_step/run_dam.py

@@ -6 +6 @@
 
 Duncan Gray, GA - 2007 
-
-
-
 """
 
@@ -16 +13 @@
 #----------------------------------------------------------------------------
 
-# Standard modules
-import time
-from time import localtime, strftime
-import sys
-from shutil import copy
-from os import path, sep
-from os.path import dirname  #, basename
+# Related major packages
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Reflective_boundary
+from anuga.shallow_water import Dirichlet_boundary
+from anuga.shallow_water import Time_boundary
+from anuga.shallow_water import File_boundary
+from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
 
 
-# Related major packages
-from anuga.shallow_water import Domain, Reflective_boundary, \
-     Dirichlet_boundary,  Time_boundary, \
-     File_boundary, \
-     Transmissive_Momentum_Set_Stage_boundary
-from anuga.fit_interpolate.interpolate import interpolate_sww2csv
-from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, \
-     copy_code_files, file_function
-from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
-     import File_boundary_time
+from anuga.abstract_2d_finite_volumes.util import file_function
+#from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
+
+
 
 # Scenario specific imports
-import project                 # Definition of file names and polygons
 import create_mesh
+
 
 def elevation_function(x,y):
@@ -52 +44 @@
     return z 
 
-def main(friction=0.01, outputdir_name=None, is_trial_run=False):
-    basename = 'zz' + str(friction)
-    if is_trial_run is True:
-        outputdir_name += '_test'
-        yieldstep = 0.1
-        finaltime = 15.
-        maximum_triangle_area=0.01
-        #maximum_triangle_area=0.1 #OK
-    else:
-        yieldstep = 0.02
-        finaltime = 15.1
-        maximum_triangle_area=0.0001
-        
-        maximum_triangle_area=0.001
-        
-    pro_instance = project.Project([],
-                                   outputdir_name=outputdir_name,
-                                   home='.')
-    print 'The output dir is', pro_instance.outputdir
-    copy_code_files(pro_instance.outputdir,__file__,
-                    dirname(project.__file__) \
-                    + sep + project.__name__+'.py')
-    copy (pro_instance.codedir + 'run_dam.py',
-          pro_instance.outputdir + 'run_dam.py')
-    copy (pro_instance.codedir + 'create_mesh.py',
-          pro_instance.outputdir + 'create_mesh.py')
+
+#-------------------------------------------------------------------------
+# Model parameters
+#-------------------------------------------------------------------------
+friction=0.01
+maximum_triangle_area=0.01
+#maximum_triangle_area=0.1 #OK
+
+basename = 'debug_boundary'
+mesh_filename = basename + '.msh'
+
+boundary_file = 'boundary.tms'
+#boundary_file = 'boundary_mellow.tms'
+
+
+
+#-------------------------------------------------------------------------
+# Create the triangular mesh
+#-------------------------------------------------------------------------
+
+# Original unstructured mesh ()
+#create_mesh.generate(mesh_filename,
+#                     maximum_triangle_area=maximum_triangle_area)
+
+
+# Simple mesh with same dimensions
+
+xright  = 19.0
+ybottom = 0
+ytop    = 0.45
+xleft = 0.0
+#xslope = 4.0 
+
+#Outline 
+point_sw = [xleft, ybottom]
+point_se = [xright, ybottom]
+point_nw = [xleft, ytop]    
+point_ne = [xright, ytop]
+
+
+create_mesh_from_regions([point_sw,
+                          point_se,
+                          point_ne,
+                          point_nw],
+                         {'wall': [0,1,2], 'wave': [3]},
+                         maximum_triangle_area=maximum_triangle_area,
+                         filename=mesh_filename)
+
+
+#TODO
+#points, elements, boundary = rectangular_cross(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
+
+
+
+#-------------------------------------------------------------------------
+# Setup computational domain
+#-------------------------------------------------------------------------
+domain = Domain(mesh_filename, use_cache = False, verbose = True)
+print domain.statistics()
+
+domain.set_name(basename)
+domain.set_datadir('.')
+
+
+#-------------------------------------------------------------------------
+# Setup initial conditions
+#-------------------------------------------------------------------------
+
+domain.set_quantity('stage', 0.06)
+domain.set_quantity('friction', friction) 
+domain.set_quantity('elevation', elevation_function)
+
+
+print 'Available boundary tags', domain.get_boundary_tags()
+
+# Create boundary function from timeseries provided in file
+function = file_function(boundary_file, domain, verbose=True)
+
+
+Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
+Br = Reflective_boundary(domain)
+domain.set_boundary( {'wall': Br, 'wave': Bts} )
+
+
+#-------------------------------------------------------------------------
+# Evolve system through time
+#-------------------------------------------------------------------------
+
+for t in domain.evolve(0.1, 15):
+    print domain.timestepping_statistics(track_speeds=False)
+    print domain.boundary_statistics(tags=['wave'])
     
-    mesh_filename = pro_instance.meshdir + basename + '.msh'
+print 'finished'
 
-    #--------------------------------------------------------------------------
-    # Copy scripts to output directory and capture screen
-    # output to file
-    #--------------------------------------------------------------------------
-
-    # creates copy of code in output dir
-    if is_trial_run is False:
-        start_screen_catcher(pro_instance.outputdir, rank, pypar.size())
-
-    print 'USER:    ', pro_instance.user
-
-    
-    #-------------------------------------------------------------------------
-    # Create the triangular mesh
-    #-------------------------------------------------------------------------
-
-    # this creates the mesh
-    #gate_position = 12.0
-    create_mesh.generate(mesh_filename,
-                         maximum_triangle_area=maximum_triangle_area)
-
-    head,tail = path.split(mesh_filename)
-    copy (mesh_filename,
-          pro_instance.outputdir + tail )
-    #-------------------------------------------------------------------------
-    # Setup computational domain
-    #-------------------------------------------------------------------------
-    domain = Domain(mesh_filename, use_cache = False, verbose = True)
-   
-
-    print 'Number of triangles = ', len(domain)
-    print 'The extent is ', domain.get_extent()
-    print domain.statistics()
-
-    
-    domain.set_name(basename)
-    domain.set_datadir(pro_instance.outputdir)
-    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
-    domain.set_minimum_storable_height(0.001)
-    #domain.set_store_vertices_uniquely(True)  # for writting to sww
-
-    #-------------------------------------------------------------------------
-    # Setup initial conditions
-    #-------------------------------------------------------------------------
-
-    domain.set_quantity('stage', 0.06)
-    domain.set_quantity('friction', friction) 
-    domain.set_quantity('elevation', elevation_function)
-
-    
-    print 'Available boundary tags', domain.get_boundary_tags()
-
-    # Create boundary function from timeseries provided in file
-    function = file_function(project.boundary_file, domain, verbose=True)
-
-    #print dir(function)
-    #print function.time, function.quantities_range
-    #import sys; sys.exit()
-    
-    Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
-
-    Br = Reflective_boundary(domain)
-    #Bd = Dirichlet_boundary([0.3,0,0]) 
-    #Bts = Time_boundary(domain, function)
-    domain.set_boundary( {'wall': Br, 'wave': Bts} )
-
-    #-------------------------------------------------------------------------
-    # Evolve system through time
-    #-------------------------------------------------------------------------
-    t0 = time.time()
-
-    for t in domain.evolve(yieldstep, finaltime):
-        print domain.timestepping_statistics(track_speeds=False)
-        print domain.boundary_statistics(tags=['wave'])
-        
-    print 'That took %.2f seconds' %(time.time()-t0)
-    print 'finished'
-
-    points = [[2.8,0.225],  #-1.8m from SWL
-              [5.1,0.225],  #0.5m from SWL
-              [6.6,0.225],  #2m from SWL
-              [6.95,0.255], #2.35m from SWL
-              [7.6,0.255],  #3m from SWL
-              [8.2,0.255],  #3.5m from SWL
-              [9.2,0.255]  #4.5m from SWL
-              ]
-
-#     points = [[gate_position - 0.65,0.2],
-#               [gate_position - 0.55,0.2],
-#               [gate_position - 0.45,0.2],
-#               [gate_position - 0.35,0.2],
-#               [gate_position - 0.25,0.2]
-#               ]
-
-    #-------------------------------------------------------------------------
-    # Calculate gauge info
-    #-------------------------------------------------------------------------
-
-    if False:
-        interpolate_sww2csv(pro_instance.outputdir + basename +".sww",
-                            points,
-                            pro_instance.outputdir + "depth_manning_"+str(friction)+".csv",
-                            pro_instance.outputdir + "velocity_x.csv",
-                            pro_instance.outputdir + "velocity_y.csv")
-  
-    return pro_instance
-
-#-------------------------------------------------------------
-if __name__ == "__main__":
-    main( is_trial_run = True,
-         outputdir_name='Hinwood_draft')