Changeset 2450

Timestamp:

Feb 27, 2006, 12:23:18 PM (19 years ago)

Author:

duncan

Message:

Further changes to second example

File:

• documentation/user_manual/anuga_user_manual.tex (modified) (19 diffs)

documentation/user_manual/anuga_user_manual.tex

@@ -423 +423 @@
     momentum quantities assumed to be the second and third conserved
     quantities.
-    
+
     \item[Transmissive boundary]Returns same conserved quantities as
     those present in its neighbour volume.
-    
+
     \item[Dirichlet boundary]Specifies a fixed value at the
 boundary.
@@ -458 +458 @@
 series of steps indicated by the values of \code{yieldstep} and
 \code{finaltime}, which can be altered as required.
-The value of \code{yieldstep} controls the time interval between successive model outputs. 
+The value of \code{yieldstep} controls the time interval between successive model outputs.
 Behind the scenes more time steps are generally taken.
 
@@ -489 +489 @@
 \section{An Example with Real Data}
 
-The following discussion builds on the concepts introduced through \code{bedslope.py} 
-example and introduces a second example, \code{run_sydney_smf.py}, that follows
-the same basic outline, but incorporates many more
-complex features.
-
-The chief difference is that data is taken from a file and represents an actual physical 
-scenario rather than an illustrative example. However a different 
-method is used to create the mesh. Instead of imposing a mesh
-structure on a rectangular grid, the technique used for this example involves building
-mesh structures inside polygons.
-
-In its simplest form, the mesh is created within a single polygon
-whose vertices are at geographical locations specified by the user.  A
-triangular mesh is created using points inside the polygon selected
-through a random process, the user specifying the
-\emph{resolution}---that is, the maximal area of a triangle used for
-triangulation.
+The following discussion builds on the concepts introduced through
+the \code{bedslope.py} example and introduces a second example,
+\code{run_sydney_smf.py}, that follows the same basic outline, but
+incorporates many more complex features.
+
+Here is the code for \code{run_sydney_smf.py}:
+
+{\scriptsize \begin{verbatim}
+"""Script for running a tsunami
+inundation scenario for Sydney, NSW, Australia.
+
+Source data such as elevation and boundary data is assumed to be
+available in directories specified by project.py The output sww file
+is stored in project.outputdir
+
+The scenario is defined by a triangular mesh created from
+project.polygon, the elevation data and a simulated submarine
+landslide.
+
+Ole Nielsen and Duncan Gray, GA - 2005 and Adrian Hitchman and Jane
+Sexton, GA - 2006 """
+
+
+#-------------------------------------------------------------------------------#
+Import necessary modules
+#-------------------------------------------------------------------------------
+
+# Standard modules
+import os
+import time
+
+# Related major packages
+from pyvolution.shallow_water import
+Domain, Reflective_boundary
+from pyvolution.data_manager import
+convert_dem_from_ascii2netcdf, dem2pts
+from pyvolution.combine_pts
+import combine_rectangular_points_files
+from pyvolution.pmesh2domain
+import pmesh_to_domain_instance
+
+# Application specific imports
+import project                 #
+Definition of file names and polygons
+from smf import slump_tsunami
+# Function for submarine mudslide
+
+
+#-------------------------------------------------------------------------------
+# Preparation of topographic data #
+# Convert ASC 2 DEM 2 PTS using
+source data and store result in source data
+# Do for coarse and fine
+data
+
+# Fine pts file to be clipped to area of interest
+#-------------------------------------------------------------------------------
+
+# filenames
+coarsedemname = project.coarsedemname
+finedemname =
+project.finedemname meshname = project.meshname+'.msh'
+
+# coarse data convert_dem_
+from_ascii2netcdf(coarsedemname,
+use_cache=True, verbose=True) dem2pts(coarsedemname, use_cache=True,
+verbose=True)
+
+# fine data (clipping the points file to smaller area)
+convert_dem_from_ascii2netcdf(finedemname, use_cache=True,
+verbose=True) dem2pts(finedemname,
+        easting_min=project.eastingmin,
+        easting_max=project.eastingmax,
+        northing_min=project.northingmin,
+        northing_max= project.northingmax,
+        use_cache=True,
+        verbose=True)
+
+
+# combining the coarse and fine data
+combine_rectangular_points_files(project.finedemname + '.pts',
+                                 project.coarsedemname + '.pts',
+                                 project.combineddemname + '.pts')
+
+
+#-------------------------------------------------------------------------------
+# Create the triangular mesh based on overall clipping polygon with
+a tagged # boundary and interior regions defined in project.py along
+with # resolutions (maximal area of per triangle) for each polygon
+#-------------------------------------------------------------------------------
+
+#def get_polygon_from_file(filename): #    fid = open(filename) #
+lines = fid.readlines() #    fid.close()
+
+#    polygon = [] #    for line in lines[1:]: #       fields =
+line.split(',') #        x = float(fields[3]) #        y =
+float(fields[4]) #        polygon.append([x, y])
+
+#    return interior_polygon
+
+from pmesh.create_mesh import create_mesh_from_regions
+
+num_polygons = 1 filelist = [] fileext = '.csv' filename =
+project.polygonptsfile interior_res = 50000
+
+for p in range(1, num_polygons+1):
+    thefilename = filename + str(p) + fileext
+    interior_polygon = get_polygon_from_file(thefilename)
+    interior_regions.append([interior_polygon, interior_res])
+
+# test #interior_regions = [[project.poly1, interior_res], #
+[project.poly2, interior_res]] # original #interior_regions =
+[[project.harbour_polygon_2, interior_res], #
+[project.botanybay_polygon_2, interior_res]]
+
+#FIXME: Fix caching of this one once the mesh_interface is ready
+from caching import cache _ = cache(create_mesh_from_regions,
+          project.diffpolygonall,
+          {'boundary_tags': {'bottom': [0],
+                             'right1': [1], 'right0': [2],
+                             'right2': [3], 'top': [4], 'left1': [5],
+                             'left2': [6], 'left3': [7]},
+           'resolution': 100000,
+           'filename': meshname,
+           'interior_regions': interior_regions,
+           'UTM': True,
+           'refzone': project.refzone},
+          verbose = True)
+
+
+#-------------------------------------------------------------------------------
+# Setup computational domain
+#-------------------------------------------------------------------------------
+
+domain = pmesh_to_domain_instance(meshname, Domain,
+                                  use_cache = True,
+                                  verbose = True)
+
+print 'Number of triangles = ', len(domain) print 'The extent is ',
+domain.get_extent()
+
+domain.set_name(project.basename)
+domain.set_datadir(project.outputdir)
+domain.set_quantities_to_be_stored(['stage', 'xmomentum',
+'ymomentum'])
+
+
+#-------------------------------------------------------------------------------
+# Set up scenario (tsunami_source is a callable object used with
+set_quantity)
+#-------------------------------------------------------------------------------
+
+tsunami_source = slump_tsunami(length=30000.0,
+                               depth=400.0,
+                               slope=6.0,
+                               thickness=176.0,
+                               radius=3330,
+                               dphi=0.23,
+                               x0=project.slump_origin[0],
+                               y0=project.slump_origin[1],
+                               alpha=0.0,
+                               domain=domain)
+
+
+#-------------------------------------------------------------------------------
+# Setup initial conditions
+#-------------------------------------------------------------------------------
+
+domain.set_quantity('stage', tsunami_source)
+domain.set_quantity('friction', 0.03) # supplied by Benfield
+domain.set_quantity('elevation',
+                    filename = project.combineddemname + '.pts',
+                    use_cache = True,
+                    verbose = True)
+
+
+#-------------------------------------------------------------------------------
+# Setup boundary conditions (all reflective)
+#-------------------------------------------------------------------------------
+
+print 'Available boundary tags', domain.get_boundary_tags()
+
+Br = Reflective_boundary(domain) domain.set_boundary( {'bottom': Br,
+'right1': Br, 'right0': Br,
+                      'right2': Br, 'top': Br, 'left1': Br,
+                      'left2': Br, 'left3': Br} )
+
+
+#-------------------------------------------------------------------------------
+# Evolve system through time
+#-------------------------------------------------------------------------------
+
+import time t0 = time.time()
+
+for t in domain.evolve(yieldstep = 120, finaltime = 18000):
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'bottom')
+
+print 'That took %.2f seconds' %(time.time()-t0)
+
+\end{verbatim}}
+
+Before discussing the details of this example, let us take a more
+high-level perspective of the various tasks it undertakes.
+
+If we compare this example with \code{bedslope.py}, the most
+noticeable difference is that data is taken from a file. This is to
+be expected, since the example now represents an actual physical
+scenario rather than an illustrative example.
+
+However, another important difference is that this example uses a
+different method to create the mesh---one more suited to use for a
+physical scenario. Instead of imposing a mesh structure on a
+rectangular grid, the technique used for this example involves
+building mesh structures inside polygons specified by the user.
+
+In its simplest form, this technique creates the mesh within a
+single polygon whose vertices are at geographical locations
+specified by the user. A triangular mesh is created using points
+inside the polygon selected through a random process, the user
+specifying the \emph{resolution}---that is, the maximal area of a
+triangle used for triangulation.
 
 Figure XXX shows a simple example, in which the triangulation is
@@ -523 +728 @@
 
 
+
+
+
+
+
 \chapter{\anuga Public Interface}
 
-This chapter lists the functions and classes available at the public interface. 
+This chapter lists the functions and classes available at the public interface.
 
 \section{Functions and Classes}
@@ -531 +741 @@
 \indexedcodeheader{create_mesh_from_region}
 
-  Creates a triangular mesh based on a bounding polygon and 
-  a number of internal polygons. For each polygon the user specifies a resolution---that is, the maximal area 
-  of triangles in the mesh. The bounding polygon also has symbolic \code{tags} associated with it.   
-  
-  \textbf{Arguments:} 
-  
+  Creates a triangular mesh based on a bounding polygon and
+  a number of internal polygons. For each polygon the user specifies a resolution---that is, the maximal area
+  of triangles in the mesh. The bounding polygon also has symbolic \code{tags} associated with it.
+
+  \textbf{Arguments:}
+
   \begin{itemize}
     \item the bounding polygon, %do we need to spell out how a polygon is specified?
-    \item a dictionary of boundary tags, for all segments of the bounding polygon, 
+    \item a dictionary of boundary tags, for all segments of the bounding polygon,
     \emph{[not clear what the keys and values of this dictionary are]}
     \item the resolution for the bounding polygon,
     \item (optional) a filename,  \emph{[what is the file for?]}
-    \item a list of 2-tuples \code{(polygon, resolution)}, specifying the interior polygons and 
+    \item a list of 2-tuples \code{(polygon, resolution)}, specifying the interior polygons and
     their associated resolutions.
   \end{itemize}
-  
-    
+
+
 \indexedcodeheader{pmesh_to_domain_instance}
 
- Converts a generated mesh file to a domain object. 
-  
-  \textbf{Arguments:} 
-  
+ Converts a generated mesh file to a domain object.
+
+  \textbf{Arguments:}
+
   \begin{itemize}
     \item \code{file_name} is the name of the mesh file to convert, including the extension
@@ -560 +770 @@
     \item \code{use_cache}: \code{True} means that caching is attempted for the computed domain.
   \end{itemize}
-  
-  
+
+
   \begin{itemize}
     \item Mesh file name
-    \item Class name, specifying the domain class to be instantiated. 
+    \item Class name, specifying the domain class to be instantiated.
   \end{itemize}
-  
-\indexedcodeheader{file_function} %in util.py "High priority" 
+
+\indexedcodeheader{file_function} %in util.py "High priority"
 
   Reads the time history of spatial data from NetCDF file and returns a callable object.
 
   \textbf{Input variables:}
-    
+
     \code{filename} - Name of \code{sww} or \code{tms} file (see
     Section \ref{sec:file formats} on page \pageref{sec:file formats}
     for details about file formats).
-       
+
        \begin{quote}
        If the file has extension \code{sww} then it is assumed to be spatio-temporal
@@ -589 +799 @@
        \end{quote}
 
-    \code{domain} - Associated domain object   
+    \code{domain} - Associated domain object
        If domain is specified, model time (\code{domain.starttime})
        will be checked and possibly modified.
-    
+
        \begin{quote}
        All times are assumed to be in UTC.
-       
+
        All spatial information is assumed to be in absolute UTM coordinates.
        \end{quote}
@@ -601 +811 @@
     \code{quantities} - the name of the quantity to be interpolated or a
                  list of quantity names. The resulting function will return
-                 a tuple of values -- one for each quantity.  
+                 a tuple of values -- one for each quantity.
 
     \code{interpolation_points} - list of absolute UTM coordinates for points at
     which values are sought
-    
+
     \code{use_cache}: \code{True} means that caching of intermediate result of
                \code{Interpolation_function} is attempted
 
-    
-  %  See Interpolation function for further documentation 
+
+  %  See Interpolation function for further documentation
 \indexedcodeheader{Interpolation_function}
  Creates a callable object \code{f(t, id)} or \code{f(t,x,y)}
@@ -617 +827 @@
 
     Let $m$ be the number of vertices, $n$ the number of triangles
-    and $p$ the number of time steps. 
+    and $p$ the number of time steps.
 
     \textbf{Mandatory input:}
-    
+
         \begin{tabular}{ll}
         \code{time}: & $p \times 1$ array of monotonously increasing times (Float)\\
-        
+
         \code{quantities}: & Dictionary of arrays or one array (Float). The arrays must  \\
         & have dimensions either $p \times m$ or $m \times 1$. The resulting function \\
@@ -629 +839 @@
         & in the latter case.\\
         \end{tabular}
-        
-        
+
+
     \textbf{Optional input:}
-    
+
         \begin{tabular}{ll}
         \code{quantity_names}: & List of keys into the quantities dictionary\\
-        
+
         \code{vertex_coordinates}: & $m \times 2$ array of coordinates (Float)\\
-        
+
         \code{triangles}: & $n \times 3$ array of indices into \code{vertex_coordinates} (Int)\\
-        
+
         \code{interpolation_points}: & $N \times 2$ array of coordinates to be interpolated to \\
-        
+
         \code{verbose}: & Level of reporting\\
         \end{tabular}
-    
+
     The quantities returned by the callable object are specified by
     the list quantities which must contain the names of the
@@ -654 +864 @@
     quantities are to be computed whenever the object is called.
     If \code{None}, returns average value.
-   
-
-\indexedcodeheader{set_region} ``Low priority. Will be merged into set\_quantity'' 
-   
+
+
+\indexedcodeheader{set_region} ``Low priority. Will be merged into set\_quantity''
+
 \indexedcodeheader{set_quantity} ``Pretty mature''
-  
+
 \indexedcodeheader{set_boundary} ``Pretty mature''
-  
+
 
 
@@ -755 +965 @@
   \begin{array}{ccr}
   2 & 4 & 4\\
-  1 & 1 & 1  
-  \end{array}  
+  1 & 1 & 1
+  \end{array}
   \right]
 \]
@@ -825 +1035 @@
 \section{caching}
 
-  The \code{cache} function is used to provide supervised caching of function results. A Python 
+  The \code{cache} function is used to provide supervised caching of function results. A Python
   function call of the form
 
@@ -838 +1048 @@
       result = cache(func,(arg1,...,argn))
       \end{verbatim}}
-  
+
   which returns the same output but reuses cached
   results if the function has been computed previously in the same context.
   \code{result} and the arguments can be simple types, tuples, list, dictionaries or
-  objects, but not unhashable types such as functions or open file objects. 
+  objects, but not unhashable types such as functions or open file objects.
   The function \code{func} may be a member function of an object or a module.
 
@@ -852 +1062 @@
   If the function definition changes after a result has been cached it will be
   detected by examining the functions \code{bytecode (co_code, co_consts,
-  func_defualts, co_argcount)} and it will be recomputed. 
-  
-  Options are set 
-  by means of the function \code{set_option(key, value)}, where \code{key} is a key associated with a 
-  Python dictionary \code{options} that stores settings such as the name of the directory used, the maximum 
+  func_defualts, co_argcount)} and it will be recomputed.
+
+  Options are set
+  by means of the function \code{set_option(key, value)}, where \code{key} is a key associated with a
+  Python dictionary \code{options} that stores settings such as the name of the directory used, the maximum
   number of cached files allowed, and so on.
-  
+
   The \code{cache} function allows the user also to specify a list of dependent files. If any of these
   have been changed, the function is recomputed and the results stored again.
-  
+
   Other features include support for compression and a capability to \ldots
 
-  
+
    \textbf{USAGE:}
-  
+
     {\scriptsize \begin{verbatim}
     result = cache(func, args, kwargs, dependencies, cachedir, verbose,
@@ -873 +1083 @@
 
   \textbf{ARGUMENTS:}
-  
+
   \begin{tabular}{ll}
     \code{func} & Function object (Required)\\
     \code{args} & Arguments to func (Default: ())\\
-    \code{kwargs} & Keyword arguments to func (Default: {})  \\  
+    \code{kwargs} & Keyword arguments to func (Default: {})  \\
     \code{dependencies} & Filenames that func depends on (Default: \code{None})\\
     \code{cachedir} & Directory for cache files (Default: \code{options['cachedir']})\\
@@ -885 +1095 @@
     \code{evaluate} & Flag forced evaluation of func (Default: 0)\\
     \code{test} & Flag test for cached results (Default: 0)\\
-    \code{clear} & Flag delete cached results (Default: 0)\\    
-    \code{return_filename} & Flag return of cache filename (Default: 0)\\    
+    \code{clear} & Flag delete cached results (Default: 0)\\
+    \code{return_filename} & Flag return of cache filename (Default: 0)\\
   \end{tabular}
- 
+
 
   \textbf{LIMITATIONS:}
-  
+
   \begin{itemize}
    \item Caching uses the apply function and will work with anything that can be
-      pickled, so any limitation in apply or pickle extends to caching. 
-      
+      pickled, so any limitation in apply or pickle extends to caching.
+
    \item A function to be cached should not depend on global variables
       as wrong results may occur if globals are changed after a result has
@@ -901 +1111 @@
    \end{itemize}
 
- 
+
 
 
@@ -1055 +1265 @@
     \item \indexedbold{resolution}   refers to the maximal area of each triangular cell in the mesh
 
-    \item \indexedbold{polygon} A sequence of points in the plane. (Arbitrary polygons can be created 
+    \item \indexedbold{polygon} A sequence of points in the plane. (Arbitrary polygons can be created
     in this way )
-    ANUGA represents polygons as either a list of 2-tuples, where the latter are either Python tuples 
-    or Python lists of length 2. The unit square, for example, would be represented by the polygon 
-    [ [0,0], [1,0], [1,1], [0,1] ]. Alternatively, polygons can be represented as $N \times 2$ Numeric 
+    ANUGA represents polygons as either a list of 2-tuples, where the latter are either Python tuples
+    or Python lists of length 2. The unit square, for example, would be represented by the polygon
+    [ [0,0], [1,0], [1,1], [0,1] ]. Alternatively, polygons can be represented as $N \times 2$ Numeric
     arrays, where $N$ is the number of points.