Changeset 9174

Timestamp:

Jun 17, 2014, 9:46:06 PM (11 years ago)

Author:

steve

Message:

Commiting all the validation tests

Location:

trunk/anuga_core

Files:

• run_auto_validation_tests.py (modified) (1 diff)
• source/anuga/shallow_water/boundaries.py (modified) (1 diff)
• source/anuga/shallow_water/test_swb2_domain.py (modified) (2 diffs)
• source/anuga_parallel/distribute_mesh.py (modified) (6 diffs)
• source/anuga_parallel/parallel_api.py (modified) (7 diffs)
• source/anuga_parallel/sequential_distribute.py (modified) (7 diffs)
• source/anuga_parallel/test_distribute_mesh.py (modified) (1 diff)
• source/anuga_parallel/test_parallel_distribute_mesh.py (modified) (2 diffs)
• source/anuga_parallel/test_parallel_frac_op.py (modified) (7 diffs)
• source/anuga_parallel/test_parallel_inlet_operator.py (modified) (8 diffs)
• source/anuga_parallel/test_sequential_dist_sw_flow.py (modified) (9 diffs)
• validation_tests/analytical_exact/avalanche_dry/numerical_avalanche_dry.py (modified) (1 diff)
• validation_tests/analytical_exact/avalanche_dry/validate_avalanche_dry.py (modified) (1 diff)
• validation_tests/analytical_exact/carrier_greenspan_transient/numerical_cg_transient.py (modified) (1 diff)
• validation_tests/analytical_exact/dam_break_dry/numerical_dam_break_dry.py (modified) (1 diff)
• validation_tests/analytical_exact/dam_break_dry/validate_dam_break_dry.py (modified) (1 diff)
• validation_tests/analytical_exact/dam_break_wet/validate_dam_break_wet.py (modified) (1 diff)
• validation_tests/analytical_exact/lake_at_rest_immersed_bump/numerical_immersed_bump.py (modified) (1 diff)
• validation_tests/analytical_exact/lake_at_rest_steep_island/numerical_steep_island.py (modified) (1 diff)
• validation_tests/analytical_exact/landslide_tsunami/produce_results.py (modified) (1 diff)
• validation_tests/analytical_exact/landslide_tsunami/runup.py (modified) (8 diffs)
• validation_tests/analytical_exact/mac_donald_short_channel/numerical_MacDonald.py (modified) (7 diffs)
• validation_tests/analytical_exact/mac_donald_short_channel/produce_results.py (modified) (1 diff)
• validation_tests/analytical_exact/parabolic_basin/numerical_parabolic_basin.py (modified) (1 diff)
• validation_tests/analytical_exact/paraboloid_basin/numerical_paraboloid_basin.py (modified) (1 diff)
• validation_tests/analytical_exact/river_at_rest_varying_topo_width/numerical_varying_width.py (modified) (4 diffs)
• validation_tests/analytical_exact/river_at_rest_varying_topo_width/plot_results.py (modified) (2 diffs)
• validation_tests/analytical_exact/river_at_rest_varying_topo_width/produce_results.py (modified) (1 diff)
• validation_tests/analytical_exact/river_at_rest_varying_topo_width/results.tex (modified) (4 diffs)
• validation_tests/analytical_exact/rundown_mild_slope/numerical_rundown_channel.py (modified) (3 diffs)
• validation_tests/analytical_exact/rundown_mild_slope/plot_results.py (modified) (4 diffs)
• validation_tests/analytical_exact/rundown_mild_slope/produce_results.py (modified) (1 diff)
• validation_tests/analytical_exact/rundown_mild_slope/results.tex (modified) (3 diffs)
• validation_tests/analytical_exact/rundown_mild_slope_coarse/numerical_rundown_channel_coarse.py (modified) (3 diffs)
• validation_tests/analytical_exact/rundown_mild_slope_coarse/plot_results.py (modified) (3 diffs)
• validation_tests/analytical_exact/rundown_mild_slope_coarse/produce_results.py (modified) (1 diff)
• validation_tests/behaviour_only/lid_driven_cavity/numerical_lid_driven_cavity.py (modified) (1 diff)
• validation_tests/behaviour_only/weir_1/runup.py (modified) (1 diff)
• validation_tests/case_studies/okushiri/run_okushiri.py (modified) (1 diff)
• validation_tests/case_studies/patong/run_model.py (modified) (2 diffs)
• validation_tests/case_studies/towradgi/Towradgi_historic_flood.py (modified) (1 diff)
• validation_tests/experimental_data/dam_break_yeh_petroff/numerical_Yeh_Petroff.py (modified) (1 diff)
• validation_tests/other_references/radial_dam_break_dry/radial_dam_break.py (modified) (1 diff)
• validation_tests/other_references/radial_dam_break_wet/radial_dam_break.py (modified) (1 diff)
• validation_tests/reports/all_tests_report.tex (modified) (5 diffs)
• validation_tests/reports/local-defs.tex (modified) (1 diff)

trunk/anuga_core/run_auto_validation_tests.py

@@ -1 +1 @@
 import os
+
 
 

trunk/anuga_core/source/anuga/shallow_water/boundaries.py

@@ -346 +346 @@
         q[1] = q[2] = 0.0
         return q
+
+class Characteristic_stage_boundary(Boundary):
+    """Sets the stage via a function and the momentum is determined 
+    via assumption of simple incoming wave (uses Riemann invariant)
+
+
+    Example:
+
+    def waveform(t):
+        return sea_level + normalized_amplitude/cosh(t-25)**2
+
+    Bcs = Characteristic_stage_boundary(domain, waveform)
+
+    Underlying domain must be specified when boundary is instantiated
+    """
+
+    def __init__(self, domain=None, function=None, default_stage = 0.0):
+        """ Instantiate a
+            Characteristic_stage_boundary.
+            domain is the domain containing the boundary
+            function is the function to apply the wave
+            default_stage is the assumed stage pre the application of wave
+        """
+
+        Boundary.__init__(self)
+
+        if domain is None:
+            msg = 'Domain must be specified for this type boundary'
+            raise Exception, msg
+
+        if function is None:
+            msg = 'Function must be specified for this type boundary'
+            raise Exception, msg
+
+        self.domain = domain
+        self.function = function
+        self.default_stage = default_stage
+
+        self.Elev  = domain.quantitis['elevation']
+        self.Stage = domain.quantitis['stage']
+        self.Height = domain.quantitis['height']
+
+    def __repr__(self):
+        """ Return a representation of this instance. """
+        msg = 'Characteristic_stage_boundary '
+        msg += '(%s) ' % self.domain
+        msg += '(%s) ' % self.default_stage
+        return msg
+
+
+    def evaluate(self, vol_id, edge_id):
+        """Calculate reflections (reverse outward momentum).
+
+        vol_id   
+        edge_id  
+        """
+        
+        t = self.domain.get_time()
+
+
+        value = self.function(t)
+        try:
+            stage = float(value)
+        except:
+            stage = float(value[0])
+
+
+
+        q = self.conserved_quantities
+        #q[0] = self.stage[vol_id, edge_id]
+        q[0] = stage
+        q[1] = self.xmom[vol_id, edge_id]
+        q[2] = self.ymom[vol_id, edge_id]
+
+        normal = self.normals[vol_id, 2*edge_id:2*edge_id+2]
+
+        r = rotate(q, normal, direction = 1)
+        r[1] = -r[1]
+        q = rotate(r, normal, direction = -1)
+
+
+        return q
+
+
+
+
+
+
+    def evaluate_segment(self, domain, segment_edges):
+        """Apply reflective BC on the boundary edges defined by
+        segment_edges
+        """
+
+        if segment_edges is None:
+            return
+        if domain is None:
+            return
+
+        t = self.domain.get_time()
+
+
+        value = self.function(t)
+        try:
+            stage = float(value)
+        except:
+            stage = float(value[0])
+                       
+
+        ids = segment_edges
+        vol_ids  = domain.boundary_cells[ids]
+        edge_ids = domain.boundary_edges[ids]
+
+        Stage = domain.quantities['stage']
+        Elev  = domain.quantities['elevation']
+        Height= domain.quantities['height']
+        Xmom  = domain.quantities['xmomentum']
+        Ymom  = domain.quantities['ymomentum']
+        Xvel  = domain.quantities['xvelocity']
+        Yvel  = domain.quantities['yvelocity']
+
+        Normals = domain.normals
+
+        #print vol_ids
+        #print edge_ids
+        #Normals.reshape((4,3,2))
+        #print Normals.shape
+        #print Normals[vol_ids, 2*edge_ids]
+        #print Normals[vol_ids, 2*edge_ids+1]
+        
+        n1  = Normals[vol_ids,2*edge_ids]
+        n2  = Normals[vol_ids,2*edge_ids+1]
+
+        # Transfer these quantities to the boundary array
+        Stage.boundary_values[ids]  = Stage.edge_values[vol_ids,edge_ids]
+        Elev.boundary_values[ids]   = Elev.edge_values[vol_ids,edge_ids]
+        Height.boundary_values[ids] = Height.edge_values[vol_ids,edge_ids]
+
+        # Rotate and negate Momemtum
+        q1 = Xmom.edge_values[vol_ids,edge_ids]
+        q2 = Ymom.edge_values[vol_ids,edge_ids]
+
+        r1 = -q1*n1 - q2*n2
+        r2 = -q1*n2 + q2*n1
+
+        Xmom.boundary_values[ids] = n1*r1 - n2*r2
+        Ymom.boundary_values[ids] = n2*r1 + n1*r2
+
+        # Rotate and negate Velocity
+        q1 = Xvel.edge_values[vol_ids,edge_ids]
+        q2 = Yvel.edge_values[vol_ids,edge_ids]
+
+        r1 = q1*n1 + q2*n2
+        r2 = q1*n2 - q2*n1
+
+        Xvel.boundary_values[ids] = n1*r1 - n2*r2
+        Yvel.boundary_values[ids] = n2*r1 + n1*r2
+
+        
+
 
 class Dirichlet_discharge_boundary(Boundary):

trunk/anuga_core/source/anuga/shallow_water/test_swb2_domain.py

@@ -36 +36 @@
 
         domain=Domain(points,vertices,boundary)    # Create Domain
-        domain.set_flow_algorithm('tsunami')
-
+        domain.set_flow_algorithm('DE0')
+        
         domain.set_name('runup_sinusoid_v2')                         # Output to file runup.sww
         domain.set_datadir('.')                          # Use current folder
@@ -91 +91 @@
 
 
+        print vv.max()
 
         assert num.all(vv<1.0e-02)

trunk/anuga_core/source/anuga_parallel/distribute_mesh.py

@@ -645 +645 @@
 
 
-    #print ghosttri
-
-    # FIXME SR: For larger layers need to pass through the correct
-    # boundary tag!
-
-#    for t in ghosttri:
-#        ghost_list.append(t[0])
-#
-#
-#
-#    for t in ghosttri:
-#
-#        n = mesh.neighbours[t[0], 0]
-#        if not is_in_processor(ghost_list, tlower, tupper, n):
-#            if boundary.has_key( (t[0], 0) ):
-#                subboundary[t[0], 0] = boundary[t[0],0]
-#            else:
-#                subboundary[t[0], 0] = 'ghost'
-#
-#
-#        n = mesh.neighbours[t[0], 1]
-#        if not is_in_processor(ghost_list, tlower, tupper, n):
-#            if boundary.has_key( (t[0], 1) ):
-#                subboundary[t[0], 1] = boundary[t[0],1]
-#            else:
-#                subboundary[t[0], 1] = 'ghost'
-#
-#
-#        n = mesh.neighbours[t[0], 2]
-#        if not is_in_processor(ghost_list, tlower, tupper, n):
-#            if boundary.has_key( (t[0], 2) ):
-#                subboundary[t[0], 2] = boundary[t[0],2]
-#            else:
-#                subboundary[t[0], 2] = 'ghost'
-#
-
-
-
     new_ghost_list = ghosttri[:,0]
 
@@ -707 +669 @@
     values1 = ['ghost']*n1
 
-    # 1 edge boundary
+    # 2 edge boundary
     nghb2 = mesh.neighbours[new_ghost_list,2]
     gl2 = num.extract(num.logical_or(nghb2 < tlower, nghb2 >= tupper), new_ghost_list)
@@ -1278 +1240 @@
 
 
+
+    tri_l2g  = extract_l2g_map(tri_map)
+    node_l2g = extract_l2g_map(node_map)
+     
     return GAnodes, GAtriangles, GAboundary, quantities, ghost_rec, \
-           full_send, tri_map, node_map, ghost_layer_width
+           full_send, tri_map, node_map, tri_l2g, node_l2g, ghost_layer_width
 
 
@@ -1615 +1581 @@
 
     [GAnodes, GAtriangles, boundary, quantities, \
-     ghost_rec, full_send, tri_map, node_map, ghost_layer_width] = \
-              build_local_mesh(submesh_cell, lower_t, upper_t, \
-                               numproc)
+     ghost_rec, full_send, \
+     tri_map, node_map, tri_l2g, node_l2g, \
+     ghost_layer_width] = \
+     build_local_mesh(submesh_cell, lower_t, upper_t, numproc)
     
     return GAnodes, GAtriangles, boundary, quantities,\
            ghost_rec, full_send,\
            number_of_full_nodes, number_of_full_triangles, tri_map, node_map,\
-           ghost_layer_width
-
-
-
-
-
-
-
-
-
+           tri_l2g, node_l2g, ghost_layer_width
 
 
@@ -1651 +1609 @@
 #
 #########################################################
-def extract_submesh(submesh, triangles_per_proc, p=0):
+def extract_submesh(submesh, triangles_per_proc, p2s_map=None, p=0):
 
     
@@ -1679 +1637 @@
     numprocs = len(triangles_per_proc)
     points, vertices, boundary, quantities, ghost_recv_dict, \
-            full_send_dict, tri_map, node_map, ghost_layer_width = \
+            full_send_dict, tri_map, node_map, tri_l2g, node_l2g, \
+            ghost_layer_width = \
             build_local_mesh(submesh_cell, lower_t, upper_t, numprocs)
 
+    if p2s_map is None:
+        pass
+    else:
+        tri_l2g = p2s_map[tri_l2g]
 
     return  points, vertices, boundary, quantities, ghost_recv_dict, \
-           full_send_dict, tri_map, node_map, ghost_layer_width
+           full_send_dict, tri_map, node_map, tri_l2g, node_l2g, ghost_layer_width
            
 
 
-
-
+def extract_l2g_map(map):
+    # Extract l2g data  from corresponding map
+    # Maps
+
+    import numpy as num
+
+    b = num.arange(len(map))
+
+    l_ids = num.extract(map>-1,map)
+    g_ids = num.extract(map>-1,b)
+
+
+#    print len(g_ids)
+#    print len(l_ids)
+#    print l_ids
+#    print g_ids
+
+    l2g = num.zeros_like(g_ids)
+    l2g[l_ids] = g_ids
+
+    return l2g
+
+
+
+
+

trunk/anuga_core/source/anuga_parallel/parallel_api.py

@@ -37 +37 @@
 
 
-
 def distribute(domain, verbose=False, debug=False, parameters = None):
+    """ Distribute the domain to all processes
+
+    parameters allows user to change size of ghost layer
+    """
+
+    if not pypar_available or numprocs == 1 : return domain # Bypass
+
+        
+    if myid == 0:
+        from sequential_distribute import Sequential_distribute
+        partition = Sequential_distribute(domain, verbose, debug, parameters)
+
+        partition.distribute(numprocs)
+
+        kwargs, points, vertices, boundary, quantities, boundary_map, \
+                domain_name, domain_dir, domain_store, domain_store_centroids, \
+                domain_minimum_storable_height, domain_minimum_allowed_height, \
+                domain_flow_algorithm, georef = partition.extract_submesh(0)
+        
+        for p in range(1, numprocs):
+
+            tostore = partition.extract_submesh(p)
+
+            send(tostore,p)
+
+    else:
+
+        kwargs, points, vertices, boundary, quantities, boundary_map, \
+                       domain_name, domain_dir, domain_store, domain_store_centroids, \
+                       domain_minimum_storable_height, domain_minimum_allowed_height, \
+                       domain_flow_algorithm, georef = receive(0)
+
+    #---------------------------------------------------------------------------
+    # Now Create parallel domain
+    #---------------------------------------------------------------------------
+    parallel_domain = Parallel_domain(points, vertices, boundary, **kwargs)
+
+
+    #------------------------------------------------------------------------
+    # Copy in quantity data
+    #------------------------------------------------------------------------
+    for q in quantities:
+        parallel_domain.set_quantity(q, quantities[q])
+
+
+    #------------------------------------------------------------------------
+    # Transfer boundary conditions to each subdomain
+    #------------------------------------------------------------------------
+    boundary_map['ghost'] = None  # Add binding to ghost boundary
+    parallel_domain.set_boundary(boundary_map)
+
+
+    #------------------------------------------------------------------------
+    # Transfer other attributes to each subdomain
+    #------------------------------------------------------------------------
+    parallel_domain.set_name(domain_name)
+    parallel_domain.set_datadir(domain_dir)
+    parallel_domain.set_store(domain_store)
+    parallel_domain.set_store_centroids(domain_store_centroids)
+    parallel_domain.set_minimum_storable_height(domain_minimum_storable_height)
+    parallel_domain.set_minimum_allowed_height(domain_minimum_allowed_height)
+    parallel_domain.set_flow_algorithm(domain_flow_algorithm)
+    parallel_domain.geo_reference = georef
+
+
+    return parallel_domain
+
+
+        
+        
+
+
+
+def old_distribute(domain, verbose=False, debug=False, parameters = None):
     """ Distribute the domain to all processes
 
@@ -125 +198 @@
                 ghost_recv_dict, full_send_dict,\
                 number_of_full_nodes, number_of_full_triangles,\
-                s2p_map, p2s_map, tri_map, node_map, ghost_layer_width =\
+                s2p_map, p2s_map, tri_map, node_map, tri_l2g, node_l2g, \
+                ghost_layer_width =\
                 distribute_mesh(domain, verbose=verbose, debug=debug, parameters=parameters)
             
         # Extract l2g maps
-        tri_l2g  = extract_l2g_map(tri_map)
-        node_l2g = extract_l2g_map(node_map)
-
+        #tri_l2g  = extract_l2g_map(tri_map)
+        #node_l2g = extract_l2g_map(node_map)
+        #tri_l2g = p2s_map[tri_l2g]
+        
         if debug:
             print 'P%d' %myid
@@ -183 +258 @@
                 ghost_recv_dict, full_send_dict,\
                 number_of_full_nodes, number_of_full_triangles, \
-                tri_map, node_map, ghost_layer_width =\
+                tri_map, node_map, tri_l2g, node_l2g, ghost_layer_width =\
                 rec_submesh(0, verbose)
 
@@ -189 +264 @@
 
         # Extract l2g maps
-        tri_l2g  = extract_l2g_map(tri_map)
-        node_l2g = extract_l2g_map(node_map)
-        
-        # Recieve serial to parallel (s2p) and parallel to serial (p2s) triangle mapping
+        #tri_l2g  = extract_l2g_map(tri_map)
+        #node_l2g = extract_l2g_map(node_map)
+        #tri_l2g = p2s_map[tri_l2g]
+        
+        # Receive serial to parallel (s2p) and parallel to serial (p2s) triangle mapping
         if send_s2p :
             s2p_map_flat = receive(0)
@@ -264 +340 @@
 
 def distribute_mesh(domain, verbose=False, debug=False, parameters=None):
-
+    """ Distribute andsend the mesh info to all the processors.
+    Should only be run from processor 0 and will send info to the other
+    processors.
+    There should be a corresponding  rec_submesh called from all the other
+    processors
+    """
 
     if debug:
@@ -303 +384 @@
     if verbose: print 'Distribute submeshes'        
     for p in range(1, numprocs):
-        send_submesh(submesh, triangles_per_proc, p, verbose)
+        send_submesh(submesh, triangles_per_proc, p2s_map, p, verbose)
 
     # Build the local mesh for processor 0
     points, vertices, boundary, quantities, \
-            ghost_recv_dict, full_send_dict, tri_map, node_map, ghost_layer_width =\
-              extract_submesh(submesh, triangles_per_proc,0)
-
+            ghost_recv_dict, full_send_dict, \
+            tri_map, node_map, tri_l2g, node_l2g, ghost_layer_width =\
+              extract_submesh(submesh, triangles_per_proc, p2s_map, 0)
+
+
+              
     # Keep track of the number full nodes and triangles.
     # This is useful later if one needs access to a ghost-free domain
@@ -344 +428 @@
            ghost_recv_dict, full_send_dict,\
            number_of_full_nodes, number_of_full_triangles, \
-           s2p_map, p2s_map, tri_map, node_map, ghost_layer_width
-    
-
-
-def extract_l2g_map(map):
-    # Extract l2g data  from corresponding map
-    # Maps
-
-    import numpy as num
-    
-    # FIXME: this is where we loss the original order of 
-    # sequential domain
-    b = num.arange(len(map))
-
-    l_ids = num.extract(map>-1,map)
-    g_ids = num.extract(map>-1,b)
-
-#    print len(g_ids)
-#    print len(l_ids)
-#    print l_ids
-
-    l2g = num.zeros_like(g_ids)
-    l2g[l_ids] = g_ids
-
-    return l2g
-
+           s2p_map, p2s_map, tri_map, node_map, tri_l2g, node_l2g, \
+           ghost_layer_width
+    
+
+
+## def extract_l2g_map(map):
+##     # Extract l2g data  from corresponding map
+##     # Maps
+
+##     import numpy as num
+    
+##     # FIXME: this is where we loss the original order of 
+##     # sequential domain
+##     b = num.arange(len(map))
+
+##     l_ids = num.extract(map>-1,map)
+##     g_ids = num.extract(map>-1,b)
+
+## #    print len(g_ids)
+## #    print len(l_ids)
+## #    print l_ids
+
+##     l2g = num.zeros_like(g_ids)
+##     l2g[l_ids] = g_ids
+
+##     return l2g
+

trunk/anuga_core/source/anuga_parallel/sequential_distribute.py

@@ -18 +18 @@
 
 
-def sequential_distribute_dump(domain, numprocs=1, verbose=False, debug=False, parameters = None):
-    """ Distribute the domain, create parallel domain and pickle result
-    """
-
-
-    if debug:
-        verbose = True
-
-
-
-    # FIXME: Dummy assignment (until boundaries are refactored to
-    # be independent of domains until they are applied)
-    bdmap = {}
-    for tag in domain.get_boundary_tags():
-        bdmap[tag] = None
-
-    domain.set_boundary(bdmap)
-
-
-    domain_name = domain.get_name()
-    domain_dir = domain.get_datadir()
-    domain_store = domain.get_store()
-    domain_minimum_storable_height = domain.minimum_storable_height
-    domain_flow_algorithm = domain.get_flow_algorithm()
-    domain_minimum_allowed_height = domain.get_minimum_allowed_height()
-    georef = domain.geo_reference
-    number_of_global_triangles = domain.number_of_triangles
-    number_of_global_nodes = domain.number_of_nodes
-    boundary_map = domain.boundary_map
-
-
-    #sequential_distribute_mesh(domain, numprocs, verbose=verbose, debug=debug, parameters=parameters)
-
-
-    # Subdivide the mesh
-    if verbose: print 'sequential_distribute: Subdivide mesh'
-    new_nodes, new_triangles, new_boundary, triangles_per_proc, quantities, \
-           s2p_map, p2s_map = \
-           pmesh_divide_metis_with_map(domain, numprocs)
-
-    #PETE: s2p_map (maps serial domain triangles to parallel domain triangles)
-    #      sp2_map (maps parallel domain triangles to domain triangles)
-
-
-
-    # Build the mesh that should be assigned to each processor,
-    # this includes ghost nodes and the communication pattern
-    if verbose: print 'sequential_distribute: Build submeshes'
-    submesh = build_submesh(new_nodes, new_triangles, new_boundary, quantities, triangles_per_proc, parameters)
-
-    if debug:
-        for p in range(numprocs):
-            N = len(submesh['ghost_nodes'][p])
-            M = len(submesh['ghost_triangles'][p])
-            print 'There are %d ghost nodes and %d ghost triangles on proc %d'\
-                  %(N, M, p)
-
-    #if debug:
-    #    from pprint import pprint
-    #    pprint(submesh)
-
-
-    # extract data to create parallel domain
-    if verbose: print 'sequential_distribute: Distribute submeshes'
-    for p in range(0, numprocs):
-
-        # Build the local mesh for processor 0
+class Sequential_distribute(object):
+
+    def __init__(self, domain, verbose=False, debug=False, parameters=None):
+
+        if debug:
+            verbose = True
+            
+        self.domain = domain
+        self.verbose = verbose
+        self.debug = debug
+        self.parameters = parameters
+
+
+    def distribute(self, numprocs=1):
+
+        self.numprocs = numprocs
+        
+        domain = self.domain
+        verbose = self.verbose
+        debug = self.debug
+        parameters = self.parameters
+
+        # FIXME: Dummy assignment (until boundaries are refactored to
+        # be independent of domains until they are applied)
+        bdmap = {}
+        for tag in domain.get_boundary_tags():
+            bdmap[tag] = None
+
+        domain.set_boundary(bdmap)
+
+
+        self.domain_name = domain.get_name()
+        self.domain_dir = domain.get_datadir()
+        self.domain_store = domain.get_store()
+        self.domain_store_centroids = domain.get_store_centroids()
+        self.domain_minimum_storable_height = domain.minimum_storable_height
+        self.domain_flow_algorithm = domain.get_flow_algorithm()
+        self.domain_minimum_allowed_height = domain.get_minimum_allowed_height()
+        self.domain_georef = domain.geo_reference
+        self.number_of_global_triangles = domain.number_of_triangles
+        self.number_of_global_nodes = domain.number_of_nodes
+        self.boundary_map = domain.boundary_map
+
+
+        # Subdivide the mesh
+        if verbose: print 'sequential_distribute: Subdivide mesh'
+
+        new_nodes, new_triangles, new_boundary, triangles_per_proc, quantities, \
+               s2p_map, p2s_map = \
+               pmesh_divide_metis_with_map(domain, numprocs)
+
+
+        # Build the mesh that should be assigned to each processor,
+        # this includes ghost nodes and the communication pattern
+        if verbose: print 'sequential_distribute: Build submeshes'
+        if verbose: print parameters
+
+        submesh = build_submesh(new_nodes, new_triangles, new_boundary, \
+                                quantities, triangles_per_proc, parameters=parameters)
+
+        if verbose:
+            for p in range(numprocs):
+                N = len(submesh['ghost_nodes'][p])
+                M = len(submesh['ghost_triangles'][p])
+                print 'There are %d ghost nodes and %d ghost triangles on proc %d'\
+                      %(N, M, p)
+
+
+        self.submesh = submesh
+        self.triangles_per_proc = triangles_per_proc
+        self.p2s_map =  p2s_map
+
+
+    def extract_submesh(self, p=0):
+        """Build the local mesh for processor p
+        """
+
+        submesh = self.submesh
+        triangles_per_proc = self.triangles_per_proc 
+        p2s_map = self.p2s_map
+        verbose = self.verbose
+        debug = self.debug
+
+        assert p>=0
+        assert p<self.numprocs
+        
+        
         points, vertices, boundary, quantities, \
-            ghost_recv_dict, full_send_dict, tri_map, node_map, ghost_layer_width =\
-              extract_submesh(submesh, triangles_per_proc, p)
-
-
-#        from pprint import pprint
-#        print '='*80
-#        print p
-#        print '='*80
-#        pprint(tri_map)
-#        print len(tri_map)
-
-        # Keep track of the number full nodes and triangles.
-        # This is useful later if one needs access to a ghost-free domain
-        # Here, we do it for process 0. The others are done in rec_submesh.
+            ghost_recv_dict, full_send_dict, \
+            tri_map, node_map, tri_l2g, node_l2g, ghost_layer_width =\
+              extract_submesh(submesh, triangles_per_proc, p2s_map, p)
+              
+
         number_of_full_nodes = len(submesh['full_nodes'][p])
         number_of_full_triangles = len(submesh['full_triangles'][p])
 
-        # Extract l2g maps
-        tri_l2g  = extract_l2g_map(tri_map)
-        node_l2g = extract_l2g_map(node_map)
-
-
+
+        if debug:
+            import pprint
+            print  50*"="
+            print 'NODE_L2G'
+            pprint.pprint(node_l2g)
+        
+            pprint.pprint(node_l2g[vertices[:,0]])
+        
+            print 'VERTICES'
+            pprint.pprint(vertices[:,0])
+            pprint.pprint(new_triangles[tri_l2g,0])
+        
+            assert num.allclose(node_l2g[vertices[:,0]], new_triangles[tri_l2g,0])        
+            assert num.allclose(node_l2g[vertices[:,1]], new_triangles[tri_l2g,1]) 
+            assert num.allclose(node_l2g[vertices[:,2]], new_triangles[tri_l2g,2]) 
+        
+
+            print 'POINTS'
+            pprint.pprint(points)
+        
+            assert num.allclose(points[:,0], new_nodes[node_l2g,0])
+            assert num.allclose(points[:,1], new_nodes[node_l2g,1])
+
+
+            print 'TRI'
+            pprint.pprint(tri_l2g)
+            pprint.pprint(p2s_map[tri_l2g])
+        
+
+            assert num.allclose(original_triangles[tri_l2orig,0],node_l2g[vertices[:,0]])
+            assert num.allclose(original_triangles[tri_l2orig,1],node_l2g[vertices[:,1]])
+            assert num.allclose(original_triangles[tri_l2orig,2],node_l2g[vertices[:,2]])
+
+            print 'NODES'
+            pprint.pprint(node_map)
+            pprint.pprint(node_l2g)      
+        
+        #tri_l2orig = p2s_map[tri_l2g]        
+        
         s2p_map = None
         p2s_map = None
@@ -118 +164 @@
             print 'sequential_distribute: P%g, no_full_nodes = %g, no_full_triangles = %g' % (p, number_of_full_nodes, number_of_full_triangles)
 
-
-        #args = [points, vertices, boundary]
 
         kwargs = {'full_send_dict': full_send_dict,
@@ -125 +169 @@
                 'number_of_full_nodes': number_of_full_nodes,
                 'number_of_full_triangles': number_of_full_triangles,
-                'geo_reference': georef,
-                'number_of_global_triangles':  number_of_global_triangles,
-                'number_of_global_nodes':  number_of_global_nodes,
+                'geo_reference': self.domain_georef,
+                'number_of_global_triangles':  self.number_of_global_triangles,
+                'number_of_global_nodes':  self.number_of_global_nodes,
                 'processor':  p,
-                'numproc':  numprocs,
+                'numproc':  self.numprocs,
                 's2p_map':  s2p_map,
                 'p2s_map':  p2s_map, ## jj added this
@@ -136 +180 @@
                 'ghost_layer_width':  ghost_layer_width}
 
-        #-----------------------------------------------------------------------
-        # Now let's store the data for a  parallel_domain via cPickle
-        #-----------------------------------------------------------------------
-
-        #Looks like we reduce storage by a factor of 4 by just
-        # storing the data to create the parallel_domain instead of pickling
-        # a created domain
+
+        boundary_map = self.boundary_map
+        domain_name = self.domain_name
+        domain_dir = self.domain_dir
+        domain_store = self.domain_store
+        domain_store_centroids = self.domain_store_centroids
+        domain_minimum_storable_height = self.domain_minimum_storable_height
+        domain_minimum_allowed_height = self.domain_minimum_allowed_height
+        domain_flow_algorithm = self.domain_flow_algorithm
+        domain_georef = self.domain_georef
+            
+        tostore = (kwargs, points, vertices, boundary, quantities, \
+                   boundary_map, \
+                   domain_name, domain_dir, domain_store, domain_store_centroids, \
+                   domain_minimum_storable_height, \
+                   domain_minimum_allowed_height, domain_flow_algorithm, \
+                   domain_georef)
+
+
+        return tostore
+
+
+
+                       
+
+    
+def sequential_distribute_dump(domain, numprocs=1, verbose=False, debug=False, parameters = None):
+    """ Distribute the domain, create parallel domain and pickle result
+    """
+
+
+    partition = Sequential_distribute(domain, verbose, debug, parameters)
+
+    partition.distribute(numprocs)
+
+    
+    for p in range(0, numprocs):
+
+        tostore = partition.extract_submesh(p) 
+
         import cPickle
-        pickle_name = domain_name + '_P%g_%g.pickle'% (numprocs,p)
+        pickle_name = partition.domain_name + '_P%g_%g.pickle'% (numprocs,p)
         f = file(pickle_name, 'wb')
-        tostore = (kwargs, points, vertices, boundary, quantities, boundary_map, domain_name, domain_dir, domain_store, domain_minimum_storable_height, \
-                   domain_minimum_allowed_height, domain_flow_algorithm, georef)
         cPickle.dump( tostore, f, protocol=cPickle.HIGHEST_PROTOCOL)
 
@@ -165 +240 @@
     pickle_name = filename+'_P%g_%g.pickle'% (numprocs,myid)
     f = file(pickle_name, 'rb')
-    kwargs, points, vertices, boundary, quantities, boundary_map, domain_name, domain_dir, domain_store, domain_minimum_storable_height, \
-                   domain_minimum_allowed_height, domain_flow_algorithm, georef = cPickle.load(f)
+
+    kwargs, points, vertices, boundary, quantities, boundary_map, \
+                   domain_name, domain_dir, domain_store, domain_store_centroids, \
+                   domain_minimum_storable_height, domain_minimum_allowed_height, \
+                   domain_flow_algorithm, georef = cPickle.load(f)
     f.close()
 
@@ -194 +272 @@
     parallel_domain.set_name(domain_name)
     parallel_domain.set_datadir(domain_dir)
+    parallel_domain.set_flow_algorithm(domain_flow_algorithm)
     parallel_domain.set_store(domain_store)
+    parallel_domain.set_store_centroids(domain_store_centroids)
     parallel_domain.set_minimum_storable_height(domain_minimum_storable_height)
     parallel_domain.set_minimum_allowed_height(domain_minimum_allowed_height)
-    parallel_domain.set_flow_algorithm(domain_flow_algorithm)
     parallel_domain.geo_reference = georef
 
@@ -203 +282 @@
     return parallel_domain
 
-def extract_l2g_map(map):
-    # Extract l2g data  from corresponding map
-    # Maps
-
-    import numpy as num
-
-    b = num.arange(len(map))
-
-    l_ids = num.extract(map>-1,map)
-    g_ids = num.extract(map>-1,b)
-
-
-#    print len(g_ids)
-#    print len(l_ids)
-#    print l_ids
-#    print g_ids
-
-    l2g = num.zeros_like(g_ids)
-    l2g[l_ids] = g_ids
-
-    return l2g
-
-
-
-
-
+
+
+
+
+
+

trunk/anuga_core/source/anuga_parallel/test_distribute_mesh.py

@@ -1502 +1502 @@
         # Now test the extract_submesh for the 3 processors
 
-        submesh_cell_0 = extract_submesh(submesh,triangles_per_proc,0)
-        submesh_cell_1 = extract_submesh(submesh,triangles_per_proc,1)
-        submesh_cell_2 = extract_submesh(submesh,triangles_per_proc,2)
+        submesh_cell_0 = extract_submesh(submesh,triangles_per_proc,p=0)
+        submesh_cell_1 = extract_submesh(submesh,triangles_per_proc,p=1)
+        submesh_cell_2 = extract_submesh(submesh,triangles_per_proc,p=2)
 
 

trunk/anuga_core/source/anuga_parallel/test_parallel_distribute_mesh.py

@@ -196 +196 @@
         #----------------------------------------------------------------------------------
         points, vertices, boundary, quantities, \
-                ghost_recv_dict, full_send_dict, tri_map, node_map, ghost_layer_width =\
+                ghost_recv_dict, full_send_dict, tri_map, node_map, tri_l2g, node_l2g, \
+                ghost_layer_width =\
         extract_submesh(submesh, triangles_per_proc)
 
@@ -226 +227 @@
         points, vertices, boundary, quantities, \
                 ghost_recv_dict, full_send_dict, \
-                no_full_nodes, no_full_trigs, tri_map, node_map, ghost_layer_width  = \
+                no_full_nodes, no_full_trigs, tri_map, node_map, tri_l2g, node_l2g, \
+                ghost_layer_width  = \
                 rec_submesh(0, verbose=False)    
 

trunk/anuga_core/source/anuga_parallel/test_parallel_frac_op.py

@@ -38 +38 @@
 from anuga.geometry.polygon import inside_polygon, is_inside_polygon, line_intersect
 
-from parallel_operator_factory import Inlet_operator, Boyd_box_operator
+#from parallel_operator_factory import Inlet_operator, Boyd_box_operator
+from anuga import Inlet_operator, Boyd_box_operator
 
 import random
@@ -48 +49 @@
 This test exercises the parallel culvert and checks values
 """
+
 verbose = False
 nprocs = 3
@@ -72 +74 @@
     for i in range(N):
 
-       # Sloping Embankment Across Channel
+        # Sloping Embankment Across Channel
         if 5.0 < x[i] < 10.1:
             # Cut Out Segment for Culvert face                
             if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
-               z[i]=z[i]
+                z[i]=z[i]
             else:
-               z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
+                z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
         if 10.0 < x[i] < 12.1:
-           z[i] +=  2.5                    # Flat Crest of Embankment
+            z[i] +=  2.5                    # Flat Crest of Embankment
         if 12.0 < x[i] < 14.5:
             # Cut Out Segment for Culvert face                
             if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
-               z[i]=z[i]
+                z[i]=z[i]
             else:
-               z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
+                z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
                    
-        
     return z
 
@@ -191 +192 @@
     if boyd_box0 is not None and verbose: boyd_box0.print_statistics()
 
-#    if parallel:
-#        factory = Parallel_operator_factory(domain, verbose = True)
-#
-#        inlet0 = factory.inlet_operator_factory(line0, Q0)
-#        inlet1 = factory.inlet_operator_factory(line1, Q1)
-#        
-#        boyd_box0 = factory.boyd_box_operator_factory(end_points=[[9.0, 2.5],[19.0, 2.5]],
-#                                          losses=1.5,
-#                                          width=1.5,
-#                                          apron=5.0,
-#                                          use_momentum_jet=True,
-#                                          use_velocity_head=False,
-#                                          manning=0.013,
-#                                          verbose=False)
-#
-#    else:
-#        inlet0 = Inlet_operator(domain, line0, Q0)
-#        inlet1 = Inlet_operator(domain, line1, Q1)
-#
-#        # Enquiry point [ 19.    2.5] is contained in two domains in 4 proc case
-#        boyd_box0 = Boyd_box_operator(domain,
-#                          end_points=[[9.0, 2.5],[19.0, 2.5]],
-#                          losses=1.5,
-#                          width=1.5,
-#                          apron=5.0,
-#                          use_momentum_jet=True,
-#                          use_velocity_head=False,
-#                          manning=0.013,
-#                          verbose=False)
-    
-    #######################################################################
 
     ##-----------------------------------------------------------------------
@@ -276 +246 @@
                 if verbose: 
                     print 'P%d tri %d, control = %s, actual = %s, Success = %s' %(myid, i, control_data[i], stage.centroid_values[tri_ids[i]], local_success) 
+                assert local_success, 'Ouput P%d tri %d, control = %s, actual = %s, Success = %s' %(myid, i, control_data[i], stage.centroid_values[tri_ids[i]], local_success) 
+
+
                 
                 
@@ -299 +272 @@
 
 
-        assert(success)
+        #assert(success)
 
     return control_data
@@ -322 +295 @@
 
 if __name__=="__main__":
+
     if numprocs == 1:
         runner = unittest.TextTestRunner()

trunk/anuga_core/source/anuga_parallel/test_parallel_inlet_operator.py

@@ -72 +72 @@
     for i in range(N):
 
-       # Sloping Embankment Across Channel
+        # Sloping Embankment Across Channel
         if 5.0 < x[i] < 10.1:
             # Cut Out Segment for Culvert face                
             if  1.0+(x[i]-5.0)/5.0 <  y[i]  < 4.0 - (x[i]-5.0)/5.0: 
-               z[i]=z[i]
+                z[i]=z[i]
             else:
-               z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
+                z[i] +=  0.5*(x[i] -5.0)    # Sloping Segment  U/S Face
         if 10.0 < x[i] < 12.1:
-           z[i] +=  2.5                    # Flat Crest of Embankment
+            z[i] +=  2.5                    # Flat Crest of Embankment
         if 12.0 < x[i] < 14.5:
             # Cut Out Segment for Culvert face                
             if  2.0-(x[i]-12.0)/2.5 <  y[i]  < 3.0 + (x[i]-12.0)/2.5:
-               z[i]=z[i]
+                z[i]=z[i]
             else:
-               z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
+                z[i] +=  2.5-1.0*(x[i] -12.0) # Sloping D/S Face
                    
         
@@ -104 +104 @@
     success = True
 
-##-----------------------------------------------------------------------
-## Setup domain
-##-----------------------------------------------------------------------
+    ##-----------------------------------------------------------------------
+    ## Setup domain
+    ##-----------------------------------------------------------------------
 
     points, vertices, boundary = rectangular_cross(int(length/dx),
@@ -117 +117 @@
     domain.set_default_order(2)
 
-##-----------------------------------------------------------------------
-## Distribute domain
-##-----------------------------------------------------------------------
+    ##-----------------------------------------------------------------------
+    ## Distribute domain
+    ##-----------------------------------------------------------------------
 
     if parallel:
@@ -126 +126 @@
     
 
-##-----------------------------------------------------------------------
-## Setup boundary conditions
-##-----------------------------------------------------------------------
+    ##-----------------------------------------------------------------------
+    ## Setup boundary conditions
+    ##-----------------------------------------------------------------------
     
     domain.set_quantity('elevation', topography) 
@@ -141 +141 @@
 
 
-##-----------------------------------------------------------------------
-## Determine triangle index coinciding with test points
-##-----------------------------------------------------------------------
+    ##-----------------------------------------------------------------------
+    ## Determine triangle index coinciding with test points
+    ##-----------------------------------------------------------------------
 
     assert(test_points is not None)
@@ -173 +173 @@
     inlet1 = Inlet_operator(domain, line1, Q1, verbose = False)
     
-    # Enquiry point [ 19.    2.5] is contained in two domains in 4 proc case
-    
-    ## boyd_box0 = Boyd_box_operator(domain,
-    ##                               end_points=[[9.0, 2.5],[19.0, 2.5]],
-    ##                               losses=1.5,
-    ##                               width=5.0,
-    ##                               apron=5.0,
-    ##                               use_momentum_jet=True,
-    ##                               use_velocity_head=False,
-    ##                               manning=0.013,
-    ##                               verbose=False)
         
     if inlet0 is not None and verbose: inlet0.print_statistics()
@@ -189 +178 @@
     #if boyd_box0 is not None and verbose: boyd_box0.print_statistics()
 
-#    if parallel:
-#        factory = Parallel_operator_factory(domain, verbose = True)
-#
-#        inlet0 = factory.inlet_operator_factory(line0, Q0)
-#        inlet1 = factory.inlet_operator_factory(line1, Q1)
-#        
-#        boyd_box0 = factory.boyd_box_operator_factory(end_points=[[9.0, 2.5],[19.0, 2.5]],
-#                                          losses=1.5,
-#                                          width=1.5,
-#                                          apron=5.0,
-#                                          use_momentum_jet=True,
-#                                          use_velocity_head=False,
-#                                          manning=0.013,
-#                                          verbose=False)
-#
-#    else:
-#        inlet0 = Inlet_operator(domain, line0, Q0)
-#        inlet1 = Inlet_operator(domain, line1, Q1)
-#
-#        # Enquiry point [ 19.    2.5] is contained in two domains in 4 proc case
-#        boyd_box0 = Boyd_box_operator(domain,
-#                          end_points=[[9.0, 2.5],[19.0, 2.5]],
-#                          losses=1.5,
-#                          width=1.5,
-#                          apron=5.0,
-#                          use_momentum_jet=True,
-#                          use_velocity_head=False,
-#                          manning=0.013,
-#                          verbose=False)
-    
-    #######################################################################
 
     ##-----------------------------------------------------------------------
@@ -246 +204 @@
     success = True
 
-##-----------------------------------------------------------------------
-## Assign/Test Control data
-##-----------------------------------------------------------------------
+    ##-----------------------------------------------------------------------
+    ## Assign/Test Control data
+    ##-----------------------------------------------------------------------
 
     if not parallel:

trunk/anuga_core/source/anuga_parallel/test_sequential_dist_sw_flow.py

@@ -33 +33 @@
 from anuga_parallel.sequential_distribute import sequential_distribute_dump
 from anuga_parallel.sequential_distribute import sequential_distribute_load
+
+import anuga.utilities.plot_utils as util
 
 
@@ -45 +47 @@
 verbose = False
 
+
+new_parameters = {}
+new_parameters['ghost_layer_width'] = 2
+
 #---------------------------------
 # Setup Functions
@@ -61 +67 @@
     if myid == 0:
         domain = rectangular_cross_domain(M, N)
-        domain.set_name('sequential_dist_runup')                    # Set sww filename
+        domain.set_name('odomain')                    # Set sww filename
         domain.set_datadir('.')   
         domain.set_quantity('elevation', topography) # Use function for elevation
@@ -74 +80 @@
     if myid == 0:
         if verbose: print 'DUMPING PARTITION DATA'
-        sequential_distribute_dump(domain, numprocs, verbose=verbose)    
+        sequential_distribute_dump(domain, numprocs, verbose=verbose, parameters=new_parameters)    
 
     #--------------------------------------------------------------------------
@@ -82 +88 @@
         
         if myid == 0 and verbose : print 'DISTRIBUTING TO PARALLEL DOMAIN'
-        pdomain = distribute(domain, verbose=verbose)
+        pdomain = distribute(domain, verbose=verbose, parameters=new_parameters)
+        pdomain.set_name('pdomain')
         
         if myid == 0 and verbose : print 'LOADING IN PARALLEL DOMAIN'
-        sdomain = sequential_distribute_load(filename='sequential_dist_runup', verbose = verbose)
-        
+        sdomain = sequential_distribute_load(filename='odomain', verbose = verbose)
+        sdomain.set_name('sdomain')
         
     if myid == 0 and verbose: print 'EVOLVING pdomain'    
@@ -94 +101 @@
     setup_and_evolve(sdomain, verbose=verbose)
     
+    if myid == 0:
+        if verbose: print 'EVOLVING odomain'   
+        setup_and_evolve(domain, verbose=verbose)
+    
+
+    if myid == 0:
+        parameter_file=open('odomain.txt', 'w')
+        from pprint import pprint
+        pprint(domain.get_algorithm_parameters(),parameter_file,indent=4)
+        parameter_file.close()
+
+        parameter_file=open('sdomain.txt', 'w')
+        from pprint import pprint
+        pprint(sdomain.get_algorithm_parameters(),parameter_file,indent=4)
+        parameter_file.close()
+
+        parameter_file=open('pdomain.txt', 'w')
+        from pprint import pprint
+        pprint(pdomain.get_algorithm_parameters(),parameter_file,indent=4)
+        parameter_file.close()        
     
     assert num.allclose(pdomain.quantities['stage'].centroid_values, sdomain.quantities['stage'].centroid_values)
@@ -103 +130 @@
     
 
-
-
-
-
+    #---------------------------------
+    # Now compare the merged sww files
+    #---------------------------------
+    if myid == 0:
+        if verbose: print 'COMPARING SWW FILES'
+        
+        odomain_v = util.get_output('odomain.sww')
+        odomain_c = util.get_centroids(odomain_v)
+
+        pdomain_v = util.get_output('pdomain.sww')
+        pdomain_c = util.get_centroids(pdomain_v)
+        
+        sdomain_v = util.get_output('sdomain.sww')
+        sdomain_c = util.get_centroids(sdomain_v)
+
+        # Test some values against the original ordering
+        
+        if verbose:
+            
+            order = 0
+            print 'PDOMAIN CENTROID VALUES'
+            print num.linalg.norm(odomain_c.x-pdomain_c.x,ord=order)
+            print num.linalg.norm(odomain_c.y-pdomain_c.y,ord=order)
+            print num.linalg.norm(odomain_c.stage[-1]-pdomain_c.stage[-1],ord=order)
+            print num.linalg.norm(odomain_c.xmom[-1]-pdomain_c.xmom[-1],ord=order)
+            print num.linalg.norm(odomain_c.ymom[-1]-pdomain_c.ymom[-1],ord=order)
+            print num.linalg.norm(odomain_c.xvel[-1]-pdomain_c.xvel[-1],ord=order)
+            print num.linalg.norm(odomain_c.yvel[-1]-pdomain_c.yvel[-1],ord=order)        
+            
+             
+            print 'SDOMAIN CENTROID VALUES'        
+            print num.linalg.norm(odomain_c.x-sdomain_c.x,ord=order)
+            print num.linalg.norm(odomain_c.y-sdomain_c.y,ord=order)
+            print num.linalg.norm(odomain_c.stage[-1]-sdomain_c.stage[-1],ord=order)
+            print num.linalg.norm(odomain_c.xmom[-1]-sdomain_c.xmom[-1],ord=order)
+            print num.linalg.norm(odomain_c.ymom[-1]-sdomain_c.ymom[-1],ord=order)
+            print num.linalg.norm(odomain_c.xvel[-1]-sdomain_c.xvel[-1],ord=order)
+            print num.linalg.norm(odomain_c.yvel[-1]-sdomain_c.yvel[-1],ord=order)
+            
+            print 'PDOMAIN VERTEX VALUES'        
+            print num.linalg.norm(odomain_v.stage[-1]-pdomain_v.stage[-1],ord=order)
+            print num.linalg.norm(odomain_v.xmom[-1]-pdomain_v.xmom[-1],ord=order)
+            print num.linalg.norm(odomain_v.ymom[-1]-pdomain_v.ymom[-1],ord=order)
+            print num.linalg.norm(odomain_v.xvel[-1]-pdomain_v.xvel[-1],ord=order)
+            print num.linalg.norm(odomain_v.yvel[-1]-pdomain_v.yvel[-1],ord=order)
+            
+            print 'SDOMAIN VERTEX VALUES'     
+            print num.linalg.norm(odomain_v.stage[-1]-sdomain_v.stage[-1],ord=order)
+            print num.linalg.norm(odomain_v.xmom[-1]-sdomain_v.xmom[-1],ord=order)
+            print num.linalg.norm(odomain_v.ymom[-1]-sdomain_v.ymom[-1],ord=order)
+            print num.linalg.norm(odomain_v.xvel[-1]-sdomain_v.xvel[-1],ord=order)
+            print num.linalg.norm(odomain_v.yvel[-1]-sdomain_v.yvel[-1],ord=order)
+            
+            
+            
+
+        assert num.allclose(odomain_c.stage,pdomain_c.stage)
+        assert num.allclose(odomain_c.xmom,pdomain_c.xmom)
+        assert num.allclose(odomain_c.ymom,pdomain_c.ymom)
+        assert num.allclose(odomain_c.xvel,pdomain_c.xvel)
+        assert num.allclose(odomain_c.yvel,pdomain_c.yvel)
+        
+        assert num.allclose(odomain_v.x,pdomain_v.x)
+        assert num.allclose(odomain_v.y,pdomain_v.y)
+                
+        assert num.linalg.norm(odomain_v.x-pdomain_v.x,ord=0) == 0
+        assert num.linalg.norm(odomain_v.y-pdomain_v.y,ord=0) == 0
+        assert num.linalg.norm(odomain_v.stage[-1]-pdomain_v.stage[-1],ord=0) < 100
+        assert num.linalg.norm(odomain_v.xmom[-1]-pdomain_v.xmom[-1],ord=0) < 100 
+        assert num.linalg.norm(odomain_v.ymom[-1]-pdomain_v.ymom[-1],ord=0) < 100
+        assert num.linalg.norm(odomain_v.xvel[-1]-pdomain_v.xvel[-1],ord=0) < 100
+        assert num.linalg.norm(odomain_v.yvel[-1]-pdomain_v.yvel[-1],ord=0) < 100     
+        
+        assert num.allclose(odomain_c.x,sdomain_c.x)
+        assert num.allclose(odomain_c.y,sdomain_c.y)
+        assert num.allclose(odomain_c.stage,sdomain_c.stage)
+        assert num.allclose(odomain_c.xmom,sdomain_c.xmom)
+        assert num.allclose(odomain_c.ymom,sdomain_c.ymom)
+        assert num.allclose(odomain_c.xvel,sdomain_c.xvel)
+        assert num.allclose(odomain_c.yvel,sdomain_c.yvel)
+        
+        assert num.allclose(odomain_v.x,sdomain_v.x)
+        assert num.allclose(odomain_v.y,sdomain_v.y)
+        
+        order = 0
+        assert num.linalg.norm(odomain_v.x-sdomain_v.x,ord=order) == 0
+        assert num.linalg.norm(odomain_v.y-sdomain_v.y,ord=order) == 0
+        assert num.linalg.norm(odomain_v.stage[-1]-sdomain_v.stage[-1],ord=order) < 100
+        assert num.linalg.norm(odomain_v.xmom[-1]-sdomain_v.xmom[-1],ord=order) < 100
+        assert num.linalg.norm(odomain_v.ymom[-1]-sdomain_v.ymom[-1],ord=order) < 100
+        assert num.linalg.norm(odomain_v.xvel[-1]-sdomain_v.xvel[-1],ord=order) < 100
+        assert num.linalg.norm(odomain_v.yvel[-1]-sdomain_v.yvel[-1],ord=order) < 100        
+                
+        # COMPARE CENTROID PDOMAIN SDOMAIN  
+        assert num.allclose(pdomain_c.x,sdomain_c.x)
+        assert num.allclose(pdomain_c.y,sdomain_c.y)
+        assert num.allclose(pdomain_c.stage[-1],sdomain_c.stage[-1])
+        assert num.allclose(pdomain_c.xmom[-1],sdomain_c.xmom[-1])
+        assert num.allclose(pdomain_c.ymom[-1],sdomain_c.ymom[-1])
+        assert num.allclose(pdomain_c.xvel[-1],sdomain_c.xvel[-1])
+        assert num.allclose(pdomain_c.yvel[-1],sdomain_c.yvel[-1])
+            
+            
+        # COMPARE VERTEX PDOMAIN SDOMAIN
+        assert num.allclose(pdomain_v.x,sdomain_v.x)
+        assert num.allclose(pdomain_v.y,sdomain_v.y)
+        assert num.allclose(pdomain_v.stage[-1],sdomain_v.stage[-1])
+        assert num.allclose(pdomain_v.xmom[-1],sdomain_v.xmom[-1])
+        assert num.allclose(pdomain_v.ymom[-1],sdomain_v.ymom[-1])
+        assert num.allclose(pdomain_v.xvel[-1],sdomain_v.xvel[-1])
+        assert num.allclose(pdomain_v.yvel[-1],sdomain_v.yvel[-1])   
+            
+
+        
+        
 def setup_and_evolve(domain, verbose=False):
 
@@ -113 +251 @@
     #--------------------------------------------------------------------------
     domain.set_flow_algorithm('DE0')
-    domain.set_quantities_to_be_stored(None)
+    #domain.set_store_vertices_uniquely()
 
     #------------------------------------------------------------------------------
@@ -131 +269 @@
     for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
         if myid == 0 and verbose : domain.write_time()
-
-
+        #if myid == 0 and verbose : print domain.quantities['stage'].get_maximum_value()
+
+
+    domain.sww_merge()
+    
 
 # Test an nprocs-way run of the shallow water equations

trunk/anuga_core/validation_tests/analytical_exact/avalanche_dry/numerical_avalanche_dry.py

@@ -103 +103 @@
 alg = args.alg
 cfl = args.cfl
-verbose = args.v
+verbose = args.verbose
 
 if myid == 0:

trunk/anuga_core/validation_tests/analytical_exact/avalanche_dry/validate_avalanche_dry.py

@@ -11 +11 @@
 
 args = anuga.get_args()
-verbose = args.v
+verbose = args.verbose
 
 class Test_results(unittest.TestCase):

trunk/anuga_core/validation_tests/analytical_exact/carrier_greenspan_transient/numerical_cg_transient.py

@@ -29 +29 @@
 alg = args.alg
 cfl = args.cfl
-verbose = args.v
+verbose = args.verbose
 
 

trunk/anuga_core/validation_tests/analytical_exact/dam_break_dry/numerical_dam_break_dry.py

@@ -52 +52 @@
 alg = args.alg
 cfl = args.cfl
-verbose = args.v
+verbose = args.verbose
 
 print args

trunk/anuga_core/validation_tests/analytical_exact/dam_break_dry/validate_dam_break_dry.py

@@ -13 +13 @@
 indent = anuga.indent
 
-verbose = args.v
+verbose = args.verbose
 
 class Test_results(unittest.TestCase):

trunk/anuga_core/validation_tests/analytical_exact/dam_break_wet/validate_dam_break_wet.py

@@ -12 +12 @@
 indent = anuga.indent
 
-verbose = args.v
+verbose = args.verbose
 
 class Test_results(unittest.TestCase):

trunk/anuga_core/validation_tests/analytical_exact/lake_at_rest_immersed_bump/numerical_immersed_bump.py

@@ -39 +39 @@
 args = anuga.get_args()
 alg = args.alg
-verbose = args.v
+verbose = args.verbose
 
 if myid == 0:

trunk/anuga_core/validation_tests/analytical_exact/lake_at_rest_steep_island/numerical_steep_island.py

@@ -29 +29 @@
 args = anuga.get_args()
 alg = args.alg
-verbose = args.v
+verbose = args.verbose
 
 dx = 1.

trunk/anuga_core/validation_tests/analytical_exact/landslide_tsunami/produce_results.py

@@ -2 +2 @@
 # import modules
 #--------------------------------
-from anuga.validation_utilities.fabricate import *
-from anuga.validation_utilities import run_validation_script
-from anuga.validation_utilities import typeset_report
+import anuga
+from anuga.validation_utilities import produce_report
 
+args = anuga.get_args()
 
-# Setup the python scripts which produce the output for this
-# validation test
-def build():
-    run_validation_script('runup.py')
-    run_validation_script('plot_results.py')
-    typeset_report()
+produce_report('runup.py', args=args)
 
-def clean():
-    autoclean()
-
-main()

trunk/anuga_core/validation_tests/analytical_exact/landslide_tsunami/runup.py

@@ -5 +5 @@
 #--------
 import anuga
+from anuga import myid, finalize, distribute
 
 import numpy
@@ -11 +12 @@
 from math import sin, pi, exp
 
+
+args = anuga.get_args()
+alg = args.alg
+verbose = args.verbose
+
 #timer = strftime('%Y%m%d_%H%M%S',localtime())
 #---------
@@ -16 +22 @@
 # --- Very inefficient mesh!
 #---------
-print ' Building mesh (alternative non-uniform mesh could be much more efficient)'
+
 dx = dy = 25.
 l1 = 60.*1000.
@@ -22 +28 @@
 nx =int(l1/dx)
 ny =int(l2/dy)
-points, vertices, boundary = anuga.rectangular_cross(nx,ny, len1=l1,len2=l2, origin=(-200., 0.))
 
-print 'Creating Domain'
-domain=anuga.Domain(points,vertices,boundary)    # Create Domain
-domain.set_name('runup_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})
-
-#------------------------------------------------------------------------------
-# Setup Algorithm, either using command line arguments
-# or override manually yourself
-#------------------------------------------------------------------------------
-from anuga.utilities.argparsing import parse_standard_args
-alg, cfl = parse_standard_args()
-domain.set_flow_algorithm(alg)
-#domain.set_CFL(cfl)
-
-#------------------
-# Define topography
-#------------------
 
 # Beach slope of 1/10
@@ -47 +34 @@
         return -(x-200.)/10. 
 
-# Define initial condition using file
-initial_prof=scipy.genfromtxt('./DATA/initial_condition.txt', skip_header=13)
-initial_prof_fun=scipy.interpolate.interp1d(initial_prof[:,0], initial_prof[:,1], fill_value=0., bounds_error=False)
+#--------------------------------------------------------------------------------
+# Create Sequential Domain
+#--------------------------------------------------------------------------------
+if myid == 0:
+        print ' Building mesh (alternative non-uniform mesh could be much more efficient)'
+        points, vertices, boundary = anuga.rectangular_cross(nx,ny, len1=l1,len2=l2, origin=(-200., 0.))
+        
+        print 'Creating Domain'
+        domain=anuga.Domain(points,vertices,boundary)    # Create Domain
+        domain.set_name('runup_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})
+        
 
-def stagefun(x,y):
-    return initial_prof_fun(x-200.) 
+        domain.set_flow_algorithm(alg)
 
-domain.set_quantity('elevation',topography)     # Use function for elevation
+        #------------------
+        # Define Initial conditions
+        #------------------
+        
+        
+        # Define initial condition using file
+        initial_prof=scipy.genfromtxt('./DATA/initial_condition.txt', skip_header=13)
+        initial_prof_fun=scipy.interpolate.interp1d(initial_prof[:,0], initial_prof[:,1], fill_value=0., bounds_error=False)
+        
+        def stagefun(x,y):
+            return initial_prof_fun(x-200.) 
+        
+        domain.set_quantity('elevation',topography)     # Use function for elevation    
+        domain.set_quantity('friction',0.000)             # Constant friction
+        domain.set_quantity('stage', stagefun)              
 
-domain.set_quantity('friction',0.000)             # Constant friction
+else:
+        
+        domain = None
 
-domain.set_quantity('stage', stagefun)              
-
+#------------------------------------------------------------------------
+# Parallel Domain
+#------------------------------------------------------------------------       
+domain = distribute(domain)
+        
 #--------------------------
 # Setup boundary conditions
@@ -73 +88 @@
 # 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()
+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()
 
 #------------------------------
@@ -85 +101 @@
 
 for t in domain.evolve(yieldstep=5.0,finaltime=350.0):
-    print domain.timestepping_statistics()
+        if myid == 0: print domain.timestepping_statistics()
     #uh=domain.quantities['xmomentum'].centroid_values
     #vh=domain.quantities['ymomentum'].centroid_values
@@ -93 +109 @@
     #print 'peak speed is', vel.max()
 
+
+domain.sww_merge(delete_old=True)
+
+finalize()

trunk/anuga_core/validation_tests/analytical_exact/mac_donald_short_channel/numerical_MacDonald.py

@@ -10 +10 @@
 import anuga
 from anuga import Domain as Domain
+from anuga import myid, finalize, distribute
+from anuga import g
 from math import cos
 from numpy import zeros, ones, array
@@ -25 +27 @@
 output_file = 'MacDonald'
 
-#anuga.copy_code_files(output_dir,__file__)
-#start_screen_catcher(output_dir+'_')
+args = anuga.get_args()
+alg = args.alg
+verbose = args.verbose
 
 
-#------------------------------------------------------------------------------
-# Setup domain
-#------------------------------------------------------------------------------
 dx = 0.25
 dy = dx
@@ -37 +37 @@
 W = 3*dx
 
-# structured mesh
-points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, 0.0))
 
-#domain = anuga.Domain(points, vertices, boundary) 
-domain = Domain(points, vertices, boundary) 
-
-domain.set_name(output_file)                
-domain.set_datadir(output_dir) 
-
-#------------------------------------------------------------------------------
-# Setup Algorithm, either using command line arguments
-# or override manually yourself
-#------------------------------------------------------------------------------
-from anuga.utilities.argparsing import parse_standard_args
-from anuga import g
-alg, cfl = parse_standard_args()
-domain.set_flow_algorithm(alg)
-#domain.set_CFL(cfl)
-
-#------------------------------------------------------------------------------
-# Setup initial conditions
-#------------------------------------------------------------------------------
 L     = 100.  # Channel length
 q = q_s = 2.    # Steady discharge
@@ -96 +75 @@
     return elevation
 
-domain.set_quantity('stage', 2.87870797)
-domain.set_quantity('elevation', bed)
-domain.set_quantity('friction', n)
 
+#------------------------------------------------------------------------------
+# Setup sequential domain
+#------------------------------------------------------------------------------
+if myid == 0:
+    # structured mesh
+    points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, 0.0))
+    
+    #domain = anuga.Domain(points, vertices, boundary) 
+    domain = Domain(points, vertices, boundary) 
+    
+    domain.set_name(output_file)                
+    domain.set_datadir(output_dir) 
+    
+    domain.set_flow_algorithm(alg)
+
+    
+    #------------------------------------------------------------------------------
+    # Setup initial conditions
+    #------------------------------------------------------------------------------
+    
+    domain.set_quantity('stage', 2.87870797)
+    domain.set_quantity('elevation', bed)
+    domain.set_quantity('friction', n)
+else:
+    
+    domain = None
+    
+#--------------------------------------------------------------------
+# Parallel Domain
+#--------------------------------------------------------------------
+domain = distribute(domain)
 
 #-----------------------------------------------------------------------------
@@ -114 +121 @@
 
 
-#===============================================================================
-##from anuga.visualiser import RealtimeVisualiser
-##vis = RealtimeVisualiser(domain)
-##vis.render_quantity_height("stage", zScale =h0*500, dynamic=True)
-##vis.colour_height_quantity('stage', (0.0, 0.5, 1.0))
-##vis.start()
-#===============================================================================
 
 
@@ -126 +126 @@
 # 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()
+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()
 
 #------------------------------------------------------------------------------
@@ -138 +139 @@
 for t in domain.evolve(yieldstep = 1.0, finaltime = 100.):
     #print domain.timestepping_statistics(track_speeds=True)
-    print domain.timestepping_statistics()
+    if myid == 0: print domain.timestepping_statistics()
     #vis.update()
 
 
-#test against know data
-    
-#vis.evolveFinished()
+domain.sww_merge(delete_old=True)
 
+finalize()
+

trunk/anuga_core/validation_tests/analytical_exact/mac_donald_short_channel/produce_results.py

@@ -2 +2 @@
 # import modules
 #--------------------------------
-from anuga.validation_utilities.fabricate import *
-from anuga.validation_utilities import run_validation_script
-from anuga.validation_utilities import typeset_report
+import anuga
+from anuga.validation_utilities import produce_report
 
+args = anuga.get_args()
 
-# Setup the python scripts which produce the output for this
-# validation test
-def build():
-    run_validation_script('numerical_MacDonald.py')
-    run_validation_script('plot_results.py')
-    typeset_report()
+produce_report('numerical_MacDonald.py', args=args)
 
-def clean():
-    autoclean()
-
-main()
-

trunk/anuga_core/validation_tests/analytical_exact/parabolic_basin/numerical_parabolic_basin.py

@@ -25 +25 @@
 args = anuga.get_args()
 alg = args.alg
-verbose = args.v
+verbose = args.verbose
 
 m = 200

trunk/anuga_core/validation_tests/analytical_exact/paraboloid_basin/numerical_paraboloid_basin.py

@@ -34 +34 @@
 args = anuga.get_args()
 alg = args.alg
-verbose = args.v
+verbose = args.verbose
 
 #-------------------------------------------------------------------------------

trunk/anuga_core/validation_tests/analytical_exact/river_at_rest_varying_topo_width/numerical_varying_width.py

@@ -10 +10 @@
 import sys
 import anuga
+from anuga import myid, finalize, distribute
 from anuga import Domain as Domain
 from math import cos
@@ -33 +34 @@
 #start_screen_catcher(output_dir+'_')
 
+args = anuga.get_args()
+alg = args.alg
+verbose = args.verbose
 
 #------------------------------------------------------------------------------
@@ -42 +46 @@
 W = 60.
 
-# structured mesh
-points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (0.,-W/2.))
 
-#domain = anuga.Domain(points, vertices, boundary) 
-domain = Domain(points, vertices, boundary) 
+#===============================================================================
+# Create sequential domain
+#===============================================================================
+if myid == 0:
+    # structured mesh
+    points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (0.,-W/2.))
+    
+    #domain = anuga.Domain(points, vertices, boundary) 
+    domain = Domain(points, vertices, boundary) 
+    
+    domain.set_name(output_file)                
+    domain.set_datadir(output_dir) 
+    
+    #------------------------------------------------------------------------------
+    # Setup Algorithm, either using command line arguments
+    # or override manually yourself
+    #------------------------------------------------------------------------------
+    from anuga.utilities.argparsing import parse_standard_args
+    alg, cfl = parse_standard_args()
+    domain.set_flow_algorithm(alg)
+    #domain.set_CFL(cfl)
+    
+    #------------------------------------------------------------------------------
+    # Setup initial conditions
+    #------------------------------------------------------------------------------
+    domain.set_quantity('friction', 0.0)
+    domain.set_quantity('stage', 12.0)
+    XX = array([0.,50.,100.,150.,250.,300.,350.,400.,425.,435.,450.,470.,475.,500.,
+                505.,530.,550.,565.,575.,600.,650.,700.,750.,800.,820.,900.,950.,
+                1000.,1500.])
+    ZZ = array([0.,0.,2.5,5.,5.,3.,5.,5.,7.5,8.,9.,9.,9.,9.1,9.,9.,6.,5.5,5.5,5.,
+                4.,3.,3.,2.3,2.,1.2,0.4,0.,0.])
+    WW = array([40.,40.,30.,30.,30.,30.,25.,25.,30.,35.,35.,40.,40.,40.,45.,45.,50.,
+                45.,40.,40.,30.,40.,40.,5.,40.,35.,25.,40.,40.])/2.
+    
+    
+    depth = interp1d(XX, ZZ)
+    width = interp1d(XX, WW)
+    
+    def bed_elevation(x,y):
+        z = 25.0*ones_like(x)
+        wid = width(x)
+        dep = depth(x)
+        z = where( logical_and(y < wid, y>-wid), dep, z)
+        return z
+    
+    domain.set_quantity('elevation', bed_elevation)
+    
+else:
+    
+    domain = None
+    
 
-domain.set_name(output_file)                
-domain.set_datadir(output_dir) 
-
-#------------------------------------------------------------------------------
-# Setup Algorithm, either using command line arguments
-# or override manually yourself
-#------------------------------------------------------------------------------
-from anuga.utilities.argparsing import parse_standard_args
-alg, cfl = parse_standard_args()
-domain.set_flow_algorithm(alg)
-#domain.set_CFL(cfl)
-
-#------------------------------------------------------------------------------
-# Setup initial conditions
-#------------------------------------------------------------------------------
-domain.set_quantity('friction', 0.0)
-domain.set_quantity('stage', 12.0)
-XX = array([0.,50.,100.,150.,250.,300.,350.,400.,425.,435.,450.,470.,475.,500.,
-            505.,530.,550.,565.,575.,600.,650.,700.,750.,800.,820.,900.,950.,
-            1000.,1500.])
-ZZ = array([0.,0.,2.5,5.,5.,3.,5.,5.,7.5,8.,9.,9.,9.,9.1,9.,9.,6.,5.5,5.5,5.,
-            4.,3.,3.,2.3,2.,1.2,0.4,0.,0.])
-WW = array([40.,40.,30.,30.,30.,30.,25.,25.,30.,35.,35.,40.,40.,40.,45.,45.,50.,
-            45.,40.,40.,30.,40.,40.,5.,40.,35.,25.,40.,40.])/2.
-
-
-depth = interp1d(XX, ZZ)
-width = interp1d(XX, WW)
-
-def bed_elevation(x,y):
-    z = 25.0*ones_like(x)
-    wid = width(x)
-    dep = depth(x)
-    z = where( logical_and(y < wid, y>-wid), dep, z)
-    return z
-
-domain.set_quantity('elevation', bed_elevation)
-
+#===========================================================================
+# Create Parallel domain
+#===========================================================================
+domain = distribute(domain)
 
 #-----------------------------------------------------------------------------
@@ -99 +117 @@
 
 
-#===============================================================================
-##from anuga.visualiser import RealtimeVisualiser
-##vis = RealtimeVisualiser(domain)
-##vis.render_quantity_height("stage", zScale =h0*500, dynamic=True)
-##vis.colour_height_quantity('stage', (0.0, 0.5, 1.0))
-##vis.start()
-#===============================================================================
-
-
-#---------------------------------------------
-# Find triangle that contains the point Point
-# and print to file
-#---------------------------------------------
-##Point = (0.0, 0.0)
-##for n in range(len(domain.triangles)):
-##    tri = domain.get_vertex_coordinates(n)
-##    if is_inside_triangle(Point,tri):
-##        #print 'Point is within triangle with vertices '+'%s'%tri
-##        n_point = n
-##        break
-##print 'n_point = ',n_point
-##bla
-
-
-
-
 #------------------------------------------------------------------------------
 # 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()
+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
 #------------------------------------------------------------------------------
+import time
+t0 = time.time()
 for t in domain.evolve(yieldstep = 0.1, finaltime = 5.0):
     #print domain.timestepping_statistics(track_speeds=True)
-    print domain.timestepping_statistics()
+    if myid == 0 and verbose: print domain.timestepping_statistics()
     #vis.update()
 
+if myid == 0 and verbose : print 'That took %s sec' % str(time.time()-t0)
+domain.sww_merge(delete_old=True)
 
-#test against know data
-    
-#vis.evolveFinished()
+finalize()
 
+

trunk/anuga_core/validation_tests/analytical_exact/river_at_rest_varying_topo_width/plot_results.py

@@ -27 +27 @@
 pyplot.ylabel('Stage')
 pyplot.xlim([0,1500])
+pyplot.ylim([0,25])
 pyplot.savefig('stage_plot.png')
 
@@ -52 +53 @@
 pyplot.xlim([0,1500])
 pyplot.savefig('xvel_plot.png')
+
+
+
+
+# Sweep in y direction
+
+v = p2_st.x[116]
+v2=(p2_st.x==v)
+
+
+#p_dev = util.get_output('dam_break.sww', 0.001)
+#p2_dev=util.get_centroids(p_dev, velocity_extrapolation=True)
+
+#Plot stages
+pyplot.clf()
+pyplot.plot(p2_st.y[v2], p2_st.stage[-1,v2],'b.', label='numerical stage')
+pyplot.plot(p2_st.y[v2], 12.*ones_like(p2_st.x[v2]),'r--', label='analytical stage')
+pyplot.plot(p2_st.y[v2], p2_st.elev[v2],'k-', label='discretised bed')
+pyplot.title('Stage at an instant in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Yposition')
+pyplot.ylabel('Stage')
+pyplot.xlim([-30,30])
+pyplot.ylim([-1,26])
+pyplot.savefig('stage_plot_y.png')
+
+
+#Plot xmomentum
+pyplot.clf()
+pyplot.plot(p2_st.y[v2], p2_st.xmom[-1,v2], 'b.', label='numerical')
+pyplot.plot(p2_st.y[v2], zeros(len(p2_st.x[v2])),'r-', label='analytical')
+pyplot.title('Xmomentum at an instant in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Yposition')
+pyplot.ylabel('Xmomentum')
+pyplot.xlim([-30,30])
+pyplot.savefig('xmom_plot_y.png')
+
+
+#Plot velocities
+pyplot.clf()
+pyplot.plot(p2_st.y[v2], p2_st.xvel[-1,v2], 'b.', label='numerical')
+pyplot.plot(p2_st.y[v2], zeros(len(p2_st.x[v2])),'r-', label='analytical')
+pyplot.title('Xvelocity at an instant in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Yposition')
+pyplot.ylabel('Xvelocity')
+pyplot.xlim([-30,30])
+pyplot.savefig('xvel_plot_y.png')

trunk/anuga_core/validation_tests/analytical_exact/river_at_rest_varying_topo_width/produce_results.py

@@ -2 +2 @@
 # import modules
 #--------------------------------
-from anuga.validation_utilities.fabricate import *
-from anuga.validation_utilities import run_validation_script
-from anuga.validation_utilities import typeset_report
+import anuga
+from anuga.validation_utilities import produce_report
 
+args = anuga.get_args()
 
-# Setup the python scripts which produce the output for this
-# validation test
-def build():
-    run_validation_script('numerical_varying_width.py')
-    run_validation_script('plot_results.py')
-    typeset_report()
+produce_report('numerical_varying_width.py', args=args)
 
-def clean():
-    autoclean()
-
-main()
-

trunk/anuga_core/validation_tests/analytical_exact/river_at_rest_varying_topo_width/results.tex

@@ -50 +50 @@
 
 
-Current version of \anuga{} might not treat the wet/dry interface appropriately. The following three figures show the stage, $x$-momentum, and $x$-velocity respectively, after we run the simulation for some time. We should see excellent agreement between the analytical and numerical solutions if the method is well-balanced and if the wet/dry interface has been correctly treated.
+Current version of \anuga{} might not treat the wet/dry interface appropriately. 
+The following figures show the stage, $x$-momentum, and $x$-velocity respectively, 
+after we run the simulation for some time. The first three show a slice in the x direction (down the river) and 
+the last three figures show a cross section across the river. We should see excellent agreement between 
+the analytical and numerical solutions if the method is well-balanced and if the wet/dry 
+interface has been correctly treated.
 
 \begin{figure}
@@ -56 +61 @@
 \includegraphics[width=0.9\textwidth]{stage_plot.png}
 \end{center}
-\caption{Stage results}
+\caption{Stage results down the river}
 \end{figure}
 
@@ -64 +69 @@
 \includegraphics[width=0.9\textwidth]{xmom_plot.png}
 \end{center}
-\caption{Xmomentum results}
+\caption{Xmomentum results down the river}
 \end{figure}
 
@@ -72 +77 @@
 \includegraphics[width=0.9\textwidth]{xvel_plot.png}
 \end{center}
-\caption{Xvelocity results}
+\caption{Xvelocity results down the river}
+\end{figure}
+
+\begin{figure}
+\begin{center}
+\includegraphics[width=0.9\textwidth]{stage_plot_y.png}
+\end{center}
+\caption{Stage results across the river}
 \end{figure}
 
 
+\begin{figure}
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xmom_plot_y.png}
+\end{center}
+\caption{Xmomentum results across the river}
+\end{figure}
+
+
+\begin{figure}
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xvel_plot_y.png}
+\end{center}
+\caption{Xvelocity results accros the river}
+\end{figure}
+
 \endinput

trunk/anuga_core/validation_tests/analytical_exact/rundown_mild_slope/numerical_rundown_channel.py

@@ -12 +12 @@
 from anuga.structures.inlet_operator import Inlet_operator
 from anuga import Domain
+from anuga import myid, finalize, distribute
 
-#------------------------------------------------------------------------------
-# Setup computational domain
-#------------------------------------------------------------------------------
-points, vertices, boundary = rectangular_cross(50, 50, len1=100.0, len2=100.0)
-domain = Domain(points, vertices, boundary) # Create domain
-domain.set_name('channel') # Output name
-
-#------------------------------------------------------------------------------
-# Setup Algorithm, either using command line arguments
-# or override manually yourself
-#------------------------------------------------------------------------------
-from anuga.utilities.argparsing import parse_standard_args
-alg, cfl = parse_standard_args()
-domain.set_flow_algorithm(alg)
-#domain.set_CFL(cfl)
-
-#------------------------------------------------------------------------------
-# Setup initial conditions
-#------------------------------------------------------------------------------
+#===============================================================================
+# Setup flow conditions
+#===============================================================================
 Qin=20.
 fluxin=Qin/100. #The momentum flux at the upstream boundary ( = discharge / width)
 uana= 2.15843634571 # analytical Xvelocity
 dana= 0.0926596702273 # analytical water depth
+slope = -0.1
+mannings = 0.03
 
-def topography(x, y):
-        return -x/10. # linear bed slope
+args = anuga.get_args()
+alg = args.alg
+verbose = args.verbose
 
-def init_stage(x,y):
-    stg= -x/10.+0.01 # Constant depth: 1 cm.
-    return stg
-#line0=[ [0.,0.], [0., 100.] ]
-#Uin=[uana, 0.0]
-#Inlet_operator(domain, line0, Q=Qin, velocity=Uin)
+#------------------------------------------------------------------------------
+# Setup sequential computational domain
+#------------------------------------------------------------------------------
+if myid == 0:
+        points, vertices, boundary = rectangular_cross(50, 50, len1=100.0, len2=100.0)
+        domain = Domain(points, vertices, boundary) # Create domain
+        domain.set_name('channel') # Output name
 
-domain.set_quantity('elevation', topography) # Use function for elevation
-domain.set_quantity('friction', 0.03) # Constant friction
-domain.set_quantity('stage', init_stage)
+        domain.set_flow_algorithm(alg)
+        domain.set_store_centroids(True)
+        
+        #------------------------------------------------------------------------------
+        # Setup initial conditions
+        #------------------------------------------------------------------------------
+        
+        def topography(x, y):
+                return x*slope # linear bed slope
+        
+        def init_stage(x,y):
+                stg= x*slope+0.01 # Constant depth: 1 cm.
+                return stg
+        
+        domain.set_quantity('elevation', topography) # Use function for elevation
+        domain.set_quantity('friction', mannings) # Constant friction
+        domain.set_quantity('stage', init_stage)
 
+else:
+        
+        domain = None
+
+#===============================================================================
+# Parallel Domain
+#===============================================================================
+domain = distribute(domain)
+
+#domain.set_store_centroids(True)
 #------------------------------------------------------------------------------
 # Setup boundary conditions
@@ -64 +77 @@
 # 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()
+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()
 
 #------------------------------------------------------------------------------
@@ -75 +89 @@
 #------------------------------------------------------------------------------
 for t in domain.evolve(yieldstep=2.0, finaltime=200.0):
-        print domain.timestepping_statistics()
+        if myid ==0 and verbose: print domain.timestepping_statistics()
+
+
+domain.sww_merge(delete_old=True)
+
+finalize()

trunk/anuga_core/validation_tests/analytical_exact/rundown_mild_slope/plot_results.py

@@ -24 +24 @@
 #------------------
 #v=(y==2.5)
-v=(p2.y==p2.y[3])
+#v=(p2.y==p2.y[3])
 
 #-------------------------------
@@ -45 +45 @@
 # Make plot animation
 #--------------------
+# Find an y value close to y==50
+tmp=(abs(p2.y-50.)).argmin()
+#vx=(p2.y==p2.y[tmp])
+vx=(abs(p2.y - p2.y[tmp])<1.5)
+
 pyplot.clf()
+pyplot.plot(p2.x[vx],p2.stage[-1,vx]-p2.elev[vx], 'o', label='numerical')
+pyplot.plot((0,100),(dana,dana),label='analytical')
+pyplot.ylim([0.05,0.1])
+pyplot.xlabel('Xposition m')
+pyplot.ylabel('Depth m')
+pyplot.title('Depth down the slope (along y=50.)')
+pyplot.legend(loc='best')
+pyplot.savefig('depth_x.png')
 
-line, = pyplot.plot( (p2.x[v].min(),p2.x[v].max()) ,( (p2.stage[:,v]-p2.elev[v]).max(),(p2.stage[:,v]-p2.elev[v]).min() ) )
-line.set_label('numerical')
-pyplot.plot( (0,100),(dana,dana), 'r',label='analytical' )
-pyplot.legend()
-#pyplot.plot(x[v],elev[v],'green')
-for i in range(p2.xmom.shape[0]):
-    line.set_xdata(p2.x[v])
-    line.set_ydata(p2.stage[i,v]-p2.elev[v])
-    pyplot.draw()
-    pyplot.title('Flow depth: should be converging to steady uniform flow ' )
-    pyplot.xlabel('Xposition')
-    pyplot.ylabel('Depth')
-pyplot.ylim([0.085,0.095])
-pyplot.legend(loc='best')
-pyplot.title('Flow depth at 400s')# -- there should be a central region with steady uniform flow ' )
-pyplot.savefig('final_depth.png')
+
+
+
+
 
 #--------------------------------------------
@@ -69 +70 @@
 # Find an x value close to x==50
 tmp=(abs(p2.x-50.)).argmin()
-v=(p2.x==p2.x[tmp])
+v=(abs(p2.x - p2.x[tmp])<1.5)
 
 pyplot.clf()
-pyplot.plot(p2.y[v],p2.xvel[100,v], 'o', label='numerical')
+pyplot.plot(p2.y[v],p2.stage[-1,v]-p2.elev[v], 'o', label='numerical')
+pyplot.plot((0,100),(dana,dana),label='analytical')
+pyplot.ylim([0.05,0.1])
+pyplot.xlabel('Yposition m')
+pyplot.ylabel('Depth m')
+pyplot.title('Depth across the slope (x=50.)')
+pyplot.legend(loc='best')
+pyplot.savefig('depth_y.png')
+
+
+pyplot.clf()
+pyplot.plot(p2.y[v],p2.xvel[-1,v], 'o', label='numerical')
 pyplot.plot((0,100),(uana,uana),label='analytical')
 pyplot.ylim([1.0,3.0])
@@ -82 +94 @@
 
 pyplot.clf()
-pyplot.plot(p2.y[v],p2.yvel[100,v],'o', label='numerical')
+pyplot.plot(p2.y[v],p2.yvel[-1,v],'o', label='numerical')
 pyplot.plot((0,100),(0.0, 0.0),label='analytical')
 pyplot.xlabel('Yposition along the line x=50')

trunk/anuga_core/validation_tests/analytical_exact/rundown_mild_slope/produce_results.py

@@ -3 +3 @@
 It was called "steep_slope" in an old validation test.
 """
-#--------------------------------
-# import modules
-#--------------------------------
 
-from anuga.validation_utilities.fabricate import *
-from anuga.validation_utilities import run_validation_script
-from anuga.validation_utilities import typeset_report
+import anuga
+from anuga.validation_utilities import produce_report
 
-# Setup the python scripts which produce the output for this
-# validation test
-def build():
-    run_validation_script('numerical_rundown_channel.py')
-    run_validation_script('plot_results.py')
-    typeset_report()
+args = anuga.get_args()
 
-def clean():
-    autoclean()
+produce_report('numerical_rundown_channel.py', args=args)
 
-main()
+
+

trunk/anuga_core/validation_tests/analytical_exact/rundown_mild_slope/results.tex

@@ -9 +9 @@
 Suppose that we are given a one dimensional domain. The steady state conditions with a contant water depth everywhere make the shallow water equations to the single identity
 \begin{equation}
-z_x = - S_f.
+z_x = - S_f, 
 \end{equation}
-Here $q=uh$ is the momentum or water discharge and $S_f$ is the symbol for the force of bottom friction involving Manning's coefficient $n$. We take 
+where $z_x$ is the bed slope, and $S_f$ is the symbol for the force of bottom friction.  We take Manning's friction
 \begin{equation}
-S_f = n^2 \frac{q|q|}{h^{10/3}}.
+S_f = n^2 \frac{q|q|}{h^{10/3}}
 \end{equation}
-If $q$, $n$, and $z_x$ are given, then the analytical solution is
+where $n$ is the Manning's coefficient and $q$ is the discharge $uh$. 
+If $q$, $n$, and $z_x$ are given, then the analytical solution for $u$ and $h$ is
 \begin{equation}
 u(x)= \left[- n^{-2} q^{4/3} z_x\right]^{3/10},
@@ -35 +36 @@
 
 
-Some simualtion results are as follows.
-Figures~\ref{fig:depthdownchan} shows the steady state depth in the downstream direction. There should be a good agreement with the analytical solution, at least away from the boundaries.  
+Some simulation results are as follows.
+Figures~\ref{fig:depthdownchan} shows the steady state depth in the downstream direction. 
+There should be a good agreement with the analytical solution, at least away from the boundaries. 
+
+Figures~\ref{fig:depthacrosschan} shows the steady state depth across the slope around the line x = 50m. 
+There should be a good agreement with the analytical solution, at least away from the boundaries. 
+
 Figures~\ref{fig:xvelscrosschan} and~\ref{fig:yvelscroschan} show the steady state $x$- and $y$-velocities, along a slice in the cross slope direction (near $x=50$). The $x$-velocities should agree well with the analytical solution, and the $y$-velocities should be zero. 
 
 \begin{figure}
 \begin{center}
-\includegraphics[width=0.8\textwidth]{final_depth.png}
+\includegraphics[width=0.8\textwidth]{depth_x.png}
 \caption{Depth in the downstream direction}
 \label{fig:depthdownchan}
@@ -47 +53 @@
 \end{figure}
  
+\begin{figure}
+\begin{center}
+\includegraphics[width=0.8\textwidth]{depth_y.png}
+\caption{Depth across the slope around $x$ = 50m}
+\label{fig:depthacrosschan}
+\end{center}
+\end{figure}
 
 \begin{figure}

trunk/anuga_core/validation_tests/analytical_exact/rundown_mild_slope_coarse/numerical_rundown_channel_coarse.py

@@ -9 +9 @@
 #------------------------------------------------------------------------------
 import anuga
-from anuga import rectangular_cross as rectangular_cross
-from anuga.structures.inlet_operator import Inlet_operator
+from anuga import rectangular_cross
+from anuga import Inlet_operator
 from anuga import Domain
+from anuga import myid, finalize, distribute
+
+Qin = 0.1
+fluxin=Qin/100. #The momentum flux at the upstream boundary ( = discharge / width)
+mann=0.03 # Manning's coef
+bedslope=-0.1
+uana= ( mann**(-2.)*abs(bedslope)*fluxin**(4./3.) )**(3./10.) # Velocity
+dana= fluxin/uana # Depth
+
+
+
+args = anuga.get_args()
+alg = args.alg
+verbose = args.verbose
 
 #------------------------------------------------------------------------------
-# Setup computational domain
+# Setup sequential computational domain
 #------------------------------------------------------------------------------
-points, vertices, boundary = rectangular_cross(14, 10, len1=140.0, len2=100.0)
-domain = Domain(points, vertices, boundary) # Create domain
-domain.set_name('channel') # Output name
+if myid == 0:
+        points, vertices, boundary = rectangular_cross(40, 10, len1=400.0, len2=100.0)
+        domain = Domain(points, vertices, boundary) # Create domain
+        domain.set_name('channel') # Output name
+        
+        domain.set_flow_algorithm(alg)
 
-#------------------------------------------------------------------------------
-# Setup Algorithm, either using command line arguments
-# or override manually yourself
-#------------------------------------------------------------------------------
-from anuga.utilities.argparsing import parse_standard_args
-alg, cfl = parse_standard_args()
-domain.set_flow_algorithm(alg)
-#domain.set_CFL(cfl)
+        
+        #------------------------------------------------------------------------------
+        # Setup initial conditions
+        #------------------------------------------------------------------------------
 
-#------------------------------------------------------------------------------
-# Setup initial conditions
-#------------------------------------------------------------------------------
-Qin=0.1
-#fluxin=Qin/100. #The momentum flux at the upstream boundary ( = discharge / width)
-#uana= 2.15843634571 # analytical Xvelocity
-#dana= 0.0926596702273 # analytical water depth
+        
+        def topography(x, y):
+                return -x/10. # linear bed slope
+        
+        def init_stage(x,y):
+                stg= -x/10.+0.004 # Constant depth: 10 cm.
+                return stg
+        #line0=[ [0.,0.], [0., 100.] ]
+        #Uin=[uana, 0.0]
+        #Inlet_operator(domain, line0, Q=Qin, velocity=Uin)
+        
 
-def topography(x, y):
-        return -x/10. # linear bed slope
-
-def init_stage(x,y):
-    stg= -x/10.+0.10 # Constant depth: 10 cm.
-    return stg
-#line0=[ [0.,0.], [0., 100.] ]
-#Uin=[uana, 0.0]
-#Inlet_operator(domain, line0, Q=Qin, velocity=Uin)
-
-line1=[ [0.,0.], [0., 100.] ]
-#Qin=0.1
-Inlet_operator(domain, line1,Qin)
-
-
-domain.set_quantity('elevation', topography) # Use function for elevation
-domain.set_quantity('friction', 0.03) # Constant friction
-domain.set_quantity('stage', init_stage)
+        
+        domain.set_quantity('elevation', topography) # Use function for elevation
+        domain.set_quantity('friction', mann) # Constant friction
+        domain.set_quantity('stage', init_stage)
+        domain.set_quantity('xmomentum', dana*uana)
+else:
+        
+        domain = None
+        
+#===========================================================================
+# Create Parallel Domain
+#===========================================================================
+domain = distribute(domain)
 
 #------------------------------------------------------------------------------
@@ -60 +72 @@
 #------------------------------------------------------------------------------
 Bt = anuga.Transmissive_boundary(domain)
+#Bts = anuga.Transmissive_momentum_set_stage_boundary(domain, dana-160.0)
+Bts = anuga.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain, lambda t: dana-160.0)
 #BdIN = anuga.Dirichlet_boundary([dana, fluxin, 0.0])
-#BdOUT = anuga.Dirichlet_boundary([dana-10., fluxin, 0.0])
+BdOUT = anuga.Dirichlet_boundary([dana-40., dana*uana, 0.0])
+
+print dana-40.
+
 Br = anuga.Reflective_boundary(domain) # Solid reflective wall
-domain.set_boundary({'left': Br, 'right': Bt, 'top': Br, 'bottom': Br})
+domain.set_boundary({'left': Br, 'right': BdOUT, 'top': Br, 'bottom': Br})
+
+
+line1=[ [0.0, 0.], [0.0, 100.] ]
+Qin=0.1
+inlet = Inlet_operator(domain, line1, Q = Qin)
+
+
+#if inlet: print inlet.statistics()
+
+
+stage = domain.quantities['stage']
+elev  = domain.quantities['elevation']
+
+print (stage-elev).get_integral()
+
+
+
+        
 
 
@@ -69 +104 @@
 # 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()
+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.0, finaltime=2000.0):
-    print domain.timestepping_statistics()
-    print (domain.areas*(domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values)).sum()
-    s3 = domain.get_flow_through_cross_section([[30., 0.0], [30., 100.]])
-    s4 = domain.get_flow_through_cross_section([[32., 0.0], [32., 100.]])
-    s5 = domain.get_flow_through_cross_section([[34., 0.0], [34., 100.]])
-    s2 = domain.get_flow_through_cross_section([[45., 0.0], [45., 100.]])
-    s1 = domain.get_flow_through_cross_section([[53., 0.0], [53., 100.]])
-    s0 = domain.get_flow_through_cross_section([[60., 0.0], [60., 100.]])
-    print 'Xsectional flow:', s0, s1, s2, s3, s4, s5
+for t in domain.evolve(yieldstep=10.0, finaltime=3000.0):
+        if myid == 0 and verbose: print domain.timestepping_statistics()
+        #print (stage-elev).get_integral()
+    #print (domain.areas*(domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values)).sum()
+    #s3 = domain.get_flow_through_cross_section([[30., 0.0], [30., 100.]])
+    #s4 = domain.get_flow_through_cross_section([[32., 0.0], [32., 100.]])
+    #s5 = domain.get_flow_through_cross_section([[34., 0.0], [34., 100.]])
+    #s2 = domain.get_flow_through_cross_section([[45., 0.0], [45., 100.]])
+    #s1 = domain.get_flow_through_cross_section([[53., 0.0], [53., 100.]])
+    #s0 = domain.get_flow_through_cross_section([[60., 0.0], [60., 100.]])
+    #print 'Xsectional flow:', s0, s1, s2, s3, s4, s5
+
+
+domain.sww_merge(delete_old=True)
+
+finalize()

trunk/anuga_core/validation_tests/analytical_exact/rundown_mild_slope_coarse/plot_results.py

@@ -26 +26 @@
 #v=(p2.y==p2.y[3])
 
-tmp = abs(p2.y-50.)
-n = tmp.argmin()
-v=(p2.y==p2.y[n])
-
-ll=p2.vel.shape[0]-1
-
 #-------------------------------
 # Define variables of case study
@@ -51 +45 @@
 # Make plot animation
 #--------------------
+# Find an y value close to y==50
+#tmp=(abs(p2.y-50.)).argmin()
+#vx=(p2.y==p2.y[tmp])
+vx=(abs(p2.y - 50.0)<10.0)
+
 pyplot.clf()
+pyplot.plot(p2.x[vx],p2.stage[-1,vx]-p2.elev[vx], 'o', label='numerical')
+pyplot.plot((0,400),(dana,dana),label='analytical')
+pyplot.ylim([0.0,0.01])
+pyplot.xlabel('Xposition m')
+pyplot.ylabel('Depth m')
+pyplot.title('Depth down the slope (along y=50.)')
+pyplot.legend(loc='best')
+pyplot.savefig('depth_x.png')
 
-line, = pyplot.plot( (p2.x[v].min(),p2.x[v].max()) ,( (p2.stage[:,v]-p2.elev[v]).max(),(p2.stage[:,v]-p2.elev[v]).min() ) )
-line.set_label('numerical')
-pyplot.plot( (0,140),(dana,dana), 'r',label='analytical' )
-pyplot.yscale('log')
-pyplot.legend()
-#pyplot.plot(x[v],elev[v],'green')
-for i in range(p2.xmom.shape[0]):
-    line.set_xdata(p2.x[v])
-    line.set_ydata(p2.stage[i,v]-p2.elev[v])
-    pyplot.draw()
-    pyplot.title('Flow depth: should be converging to steady uniform flow ' )
-    pyplot.xlabel('Xposition')
-    pyplot.ylabel('Depth')
-#pyplot.ylim([0.085,0.095])
+
+
+pyplot.clf()
+pyplot.plot(p2.x[vx],p2.xvel[-1,vx], 'o', label='numerical')
+pyplot.plot((0,400),(uana,uana),label='analytical')
+#pyplot.ylim([0.0,0.05])
+pyplot.xlabel('Xposition m')
+pyplot.ylabel('Depth m')
+pyplot.title('Velocity down the slope (along y=50.)')
 pyplot.legend(loc='best')
-pyplot.title('Flow depth at 800s')# -- there should be a central region with steady uniform flow ' )
-pyplot.savefig('final_depth.png')
+pyplot.savefig('xvelocity_x.png')
+
+
+
 
 #--------------------------------------------
@@ -75 +79 @@
 #--------------------------------------------
 # Find an x value close to x==50
-tmp=(abs(p2.x-50.)).argmin()
-v=(p2.x==p2.x[tmp])
+#tmp=(abs(p2.x-50.)).argmin()
+v=(abs(p2.x - 50.0)<10.0)
 
 pyplot.clf()
-pyplot.plot(p2.y[v],p2.xvel[ll,v], 'o', label='numerical')
+pyplot.plot(p2.y[v],p2.stage[-1,v]-p2.elev[v], 'o', label='numerical')
+pyplot.plot((0,100),(dana,dana),label='analytical')
+pyplot.ylim([0.0,0.005])
+pyplot.xlabel('Yposition m')
+pyplot.ylabel('Depth m')
+pyplot.title('Depth across the slope (x=50.)')
+pyplot.legend(loc='best')
+pyplot.savefig('depth_y.png')
+
+
+pyplot.clf()
+pyplot.plot(p2.y[v],p2.xvel[-1,v], 'o', label='numerical')
 pyplot.plot((0,100),(uana,uana),label='analytical')
-#pyplot.ylim([1.0,3.0])
+pyplot.ylim([0.0,0.35])
 pyplot.xlabel('Yposition along the line x=50')
 pyplot.ylabel('Xvelocity m/s')
 pyplot.title('Final Xvelocity around the line x=50.')
 pyplot.legend(loc='best')
-pyplot.ylim((0,1.5*uana))
 pyplot.savefig('x_velocity.png')
 
 pyplot.clf()
-pyplot.plot(p2.y[v],p2.yvel[ll,v],'o', label='numerical')
+pyplot.plot(p2.y[v],p2.yvel[-1,v],'o', label='numerical')
 pyplot.plot((0,100),(0.0, 0.0),label='analytical')
+#pyplot.ylim([0.0,0.3])
 pyplot.xlabel('Yposition along the line x=50')
 pyplot.ylabel('Yvelocity')

trunk/anuga_core/validation_tests/analytical_exact/rundown_mild_slope_coarse/produce_results.py

@@ -3 +3 @@
 It was called "steep_slope" in an old validation test.
 """
-#--------------------------------
-# import modules
-#--------------------------------
 
-from anuga.validation_utilities.fabricate import *
-from anuga.validation_utilities import run_validation_script
-from anuga.validation_utilities import typeset_report
+import anuga
+from anuga.validation_utilities import produce_report
 
-# Setup the python scripts which produce the output for this
-# validation test
-def build():
-    run_validation_script('numerical_rundown_channel_coarse.py')
-    run_validation_script('plot_results.py')
-    typeset_report()
+args = anuga.get_args()
 
-def clean():
-    autoclean()
+produce_report('numerical_rundown_channel_coarse.py', args=args)
 
-main()
+
+

trunk/anuga_core/validation_tests/behaviour_only/lid_driven_cavity/numerical_lid_driven_cavity.py

@@ -23 +23 @@
 alg, cfl = parse_standard_args()
 domain.set_flow_algorithm(alg)
-domain.set_CFL(cfl)
+#domain.set_CFL(cfl)
 
 domain.set_name('dimensional_lid_driven')   # Output to file runup.sww

trunk/anuga_core/validation_tests/behaviour_only/weir_1/runup.py

@@ -60 +60 @@
 alg, cfl = parse_standard_args()
 domain.set_flow_algorithm(alg)
-domain.set_CFL(cfl)
+#domain.set_CFL(cfl)
 
 domain.set_name('runup_riverwall')                         

trunk/anuga_core/validation_tests/case_studies/okushiri/run_okushiri.py

@@ -84 +84 @@
         alg, cfl = parse_standard_args()
         domain.set_flow_algorithm(alg)
-        domain.set_CFL(cfl)
+        #domain.set_CFL(cfl)
     
     

trunk/anuga_core/validation_tests/case_studies/patong/run_model.py

@@ -38 +38 @@
 import shutil
 
-
-verbose = True
+#--------------------------------------------------
+# Pick up useful command line arguments (which over rides 
+# values set before
+#--------------------------------------------------
+cfl = anuga.args.cfl
+alg = anuga.args.alg
+verbose = anuga.args.verbose
+np = anuga.args.np
+
+if myid == 0 and verbose:
+    print 80*'#'
+    print '#'
+    print '# Long Validation Test, takes 10 minutes on my desktop'
+    print '#'
+    print '# Consider running in parallel'
+    print '#'
+    print 80*'#'
+
 #-------------------------------------------------------------------------------
 # Copy scripts to time stamped output directory and capture screen
@@ -133 +149 @@
     domain.set_name(project.scenario_name)
     domain.set_datadir(project.output_run) 
-    print 'WARNING: FORCING FLOW ALGORITHM TO DE0 -- NEEDS TO BE CHANGED FOR INTEGRATION INTO VALIDATION TESTS'
-    #domain.set_flow_algorithm('tsunami')
-    domain.set_flow_algorithm('DE0')
+
+    domain.set_flow_algorithm(alg)
+    #domain.set_CFL(cfl)
 
     #-------------------------------------------------------------------------------

trunk/anuga_core/validation_tests/case_studies/towradgi/Towradgi_historic_flood.py

@@ -225 +225 @@
     alg, cfl = parse_standard_args()
     domain.set_flow_algorithm(alg)
-    domain.set_CFL(cfl)
+    #domain.set_CFL(cfl)
     #domain.set_flow_algorithm('DE1') 
     #domain.set_store_vertices_smoothly()

trunk/anuga_core/validation_tests/experimental_data/dam_break_yeh_petroff/numerical_Yeh_Petroff.py

@@ -58 +58 @@
 
 domain.set_flow_algorithm(alg)
-domain.set_CFL(cfl)
+#domain.set_CFL(cfl)
 #domain.set_minimum_allowed_height(0.01) # Avoid such statements in the validation tests
 

trunk/anuga_core/validation_tests/other_references/radial_dam_break_dry/radial_dam_break.py

@@ -50 +50 @@
 alg, cfl = parse_standard_args()
 domain.set_flow_algorithm(alg)
-domain.set_CFL(cfl)
+#domain.set_CFL(cfl)
 #domain.set_minimum_allowed_height(0.002)
 

trunk/anuga_core/validation_tests/other_references/radial_dam_break_wet/radial_dam_break.py

@@ -50 +50 @@
 alg, cfl = parse_standard_args()
 domain.set_flow_algorithm(alg)
-domain.set_CFL(cfl)
+#domain.set_CFL(cfl)
 
 #------------------------------------------------------------------------------

trunk/anuga_core/validation_tests/reports/all_tests_report.tex

@@ -91 +91 @@
 %======================
 
-\inputresults{analytical_exact/dam_break_dry}
-\inputresults{analytical_exact/dam_break_wet}
-\inputresults{analytical_exact/avalanche_dry}
-\inputresults{analytical_exact/avalanche_wet}
-\inputresults{analytical_exact/carrier_greenspan_periodic}
-\inputresults{analytical_exact/carrier_greenspan_transient}
-\inputresults{analytical_exact/deep_wave}
-\inputresults{analytical_exact/mac_donald_short_channel}
-\inputresults{analytical_exact/parabolic_basin}
-\inputresults{analytical_exact/paraboloid_basin}
-
-\inputresults{analytical_exact/runup_on_beach}
-\inputresults{analytical_exact/runup_on_sinusoid_beach}
-\inputresults{analytical_exact/landslide_tsunami}
-
-\inputresults{analytical_exact/lake_at_rest_immersed_bump}
-\inputresults{analytical_exact/lake_at_rest_steep_island}
-\inputresults{analytical_exact/river_at_rest_varying_topo_width}
-\inputresults{analytical_exact/rundown_mild_slope}
-\inputresults{analytical_exact/rundown_mild_slope_coarse}
-\inputresults{analytical_exact/subcritical_over_bump}
-\inputresults{analytical_exact/transcritical_with_shock}
-\inputresults{analytical_exact/transcritical_without_shock}
-\inputresults{analytical_exact/trapezoidal_channel}
+\inputresults{../analytical_exact/dam_break_dry}
+\inputresults{../analytical_exact/dam_break_wet}
+\inputresults{../analytical_exact/avalanche_dry}
+\inputresults{../analytical_exact/avalanche_wet}
+\inputresults{../analytical_exact/carrier_greenspan_periodic}
+\inputresults{../analytical_exact/carrier_greenspan_transient}
+\inputresults{../analytical_exact/deep_wave}
+\inputresults{../analytical_exact/mac_donald_short_channel}
+\inputresults{../analytical_exact/parabolic_basin}
+\inputresults{../analytical_exact/paraboloid_basin}
+
+\inputresults{../analytical_exact/runup_on_beach}
+\inputresults{../analytical_exact/runup_on_sinusoid_beach}
+\inputresults{../analytical_exact/landslide_tsunami}
+
+\inputresults{../analytical_exact/lake_at_rest_immersed_bump}
+\inputresults{../analytical_exact/lake_at_rest_steep_island}
+\inputresults{../analytical_exact/river_at_rest_varying_topo_width}
+\inputresults{../analytical_exact/rundown_mild_slope}
+\inputresults{../analytical_exact/rundown_mild_slope_coarse}
+\inputresults{../analytical_exact/subcritical_over_bump}
+\inputresults{../analytical_exact/transcritical_with_shock}
+\inputresults{../analytical_exact/transcritical_without_shock}
+\inputresults{../analytical_exact/trapezoidal_channel}
 
 %%======================
 \chapter{Tests against reference data or solutions} \label{ch:ref}
 %%======================
-\inputresults{experimental_data/dam_break_yeh_petroff}
-\inputresults{behaviour_only/lid_driven_cavity}
-
-\inputresults{other_references/radial_dam_break_dry}
-\inputresults{other_references/radial_dam_break_wet}
-
-\inputresults{case_studies/okushiri}
+\inputresults{../experimental_data/dam_break_yeh_petroff}
+\inputresults{../behaviour_only/lid_driven_cavity}
+
+\inputresults{../other_references/radial_dam_break_dry}
+\inputresults{../other_references/radial_dam_break_wet}
+
+\inputresults{../case_studies/okushiri}
 
 
@@ -151 +151 @@
 \section{Algorithm Parameters}
 Note that parameters can be communicated from the \verb|local_parameters.py|
-file in the \verb|anuga_validation_tests| directory. If there is no file
+file in the \verb|validation_tests/reports| directory. If there is no file
 \verb|local_parameters.py| then the parameters are taken from the
-\verb|parameters.py| file. You should not change this file.
+\verb|anuga.validation_utilities.parameters|. 
 
 In particular the
@@ -160 +160 @@
 test directories.
 
-You can pass though the parameters straight from the \verb|parameters.py| file as follows
+You can pass though the standard parameters as follows
 \begin{verbatim}
-from anuga_validation_tests.parameters import alg
-from anuga_validation_tests.parameters import cfl
+from anuga.validation_utilities.parameters import alg
+from anuga.validation_utilities.parameters import cfl
 \end{verbatim}
 
@@ -172 +172 @@
 
 \begin{verbatim}
-from anuga_validation_test_utilities.fabricate import *
-from anuga_validation_test_utilities import run_validation_script
+from anuga.validation_utilities.fabricate import *
+from anuga.validation_utilities import run_validation_script
 
 # Setup the python scripts which produce the output for this
@@ -180 +180 @@
     run_validation_script('run_problem.py')
     run_validation_script('plot_problem.py')
-    run('python','latex_report.py')
+    run('python','typeset_report.py')
     pass
 

trunk/anuga_core/validation_tests/reports/local-defs.tex

@@ -3 +3 @@
 
 
-\newcommand{\inputresults}[1]{\graphicspath{{#1/}}\input{#1/results}}
+\newcommand{\inputresults}[1]{\graphicspath{{../}{#1/}}\input{#1/results}}
+%\newcommand{\inputresults}[1]{\graphicspath{{../}}\input{#1/results}}