Changeset 8616

Timestamp:

Nov 13, 2012, 9:11:35 PM (11 years ago)

Author:

steve

Message:

Small changes to inlet operators

Location:

trunk/anuga_core

Files:

• source/anuga/abstract_2d_finite_volumes/mesh_factory_ext.c (modified) (1 diff)
• source/anuga/structures/inlet.py (modified) (1 diff)
• source/anuga/structures/inlet_operator.py (modified) (4 diffs)
• source/anuga_parallel/run_sequential_dist_distribute_merimbula.py (modified) (2 diffs)
• source/anuga_parallel/run_sequential_dist_evolve_merimbula.py (modified) (3 diffs)
• validation_tests/Tests/Benchmarks/Okushiri/create_okushiri.py (modified) (2 diffs)
• validation_tests/Tests/Benchmarks/Okushiri/run_okushiri.py (modified) (1 diff)

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/mesh_factory_ext.c

@@ -118 +118 @@
                 }
 
+        free(vertices);
+
         return Py_BuildValue("O", pyobj_boundary);
 }

trunk/anuga_core/source/anuga/structures/inlet.py

@@ -233 +233 @@
         summed_volume = num.zeros_like(areas)       
         summed_volume[1:] = num.cumsum(summed_areas[:-1]*num.diff(stages[stages_order]))
-        
-        # find the number of cells which will be filled
-        index = num.nonzero(summed_volume<=volume)[0][-1]
-
-        # calculate stage needed to fill chosen cells with given volume of water
-        depth = (volume - summed_volume[index])/summed_areas[index]
-        stages[stages_order[0:index+1]] = stages[stages_order[index]]+depth
-
-        self.set_stages(stages) 
+
+
+        if volume>= 0.0 :
+            # find the number of cells which will be filled
+            #print summed_volume
+            #print summed_volume<=volume
+
+            #Test = summed_volume<=volume
+
+            #print num.any(Test)
+            #print num.nonzero(summed_volume<=volume)
+            #print num.nonzero(summed_volume<=volume)[0]
+
+            index = num.nonzero(summed_volume<=volume)[0][-1]
+
+            # calculate stage needed to fill chosen cells with given volume of water
+            depth = (volume - summed_volume[index])/summed_areas[index]
+            stages[stages_order[0:index+1]] = stages[stages_order[index]]+depth
+
+        else:
+            if summed_volume[-1]>= -volume :
+                depth = (summed_volume[-1] + volume)/summed_areas[-1]
+                stages[:] = depth
+            else:
+                #assert summed_volume[-1] >= -volume
+                import warnings
+                warnings.warn('summed_volume < volume to be extracted')
+                stages[:] = 0.0
+
+        self.set_stages(stages)
+
+
+
+
+    def set_stages_evenly_new(self,volume):
+        """ Distribute volume of water over
+        inlet exchange region so that stage is level
+        """
+        # WARNING: requires synchronization, must be called by all procs associated
+        # with this inlet
+
+        centroid_coordinates = self.domain.get_full_centroid_coordinates(absolute=True)
+        areas = self.get_areas()
+        stages = self.get_stages()
+
+        stages_order = stages.argsort()
+
+        # PETE: send stages and areas, apply merging procedure
+
+        total_stages = len(stages)
+
+
+        s_areas = areas
+        s_stages = stages
+        s_stages_order = stages_order
+
+
+        pos = 0
+        summed_volume = 0.
+        summed_areas = 0.
+        prev_stage = 0.
+        num_stages = 0.
+        first = True
+
+        for i in self.procs:
+            pos[i] = 0
+
+        while num_stages < total_stages:
+            # Determine current minimum stage of all the processors in s_stages
+            num_stages = num_stages + 1
+            current_stage = num.finfo(num.float32).max
+            index = -1
+
+
+            if pos >= len(s_stages):
+                continue
+
+            if s_stages[s_stages_order[pos]] < current_stage:
+                current_stage = s_stages[i][s_stages_order[i][pos[i]]]
+                index = i
+
+            # If first iteration, then only update summed_areas, position, and prev|current stage
+
+            if first:
+                first = False
+                summed_areas = s_areas[index][s_stages_order[index][pos[index]]]
+                pos[index] = pos[index] + 1
+                prev_stage = current_stage
+                continue
+
+            assert index >= 0, "Index out of bounds"
+
+            # Update summed volume and summed areas
+            tmp_volume = summed_volume + (summed_areas * (current_stage - prev_stage))
+
+            # Terminate if volume exceeded
+            if tmp_volume >= volume:
+                break
+
+            summed_areas = summed_areas + s_areas[index][s_stages_order[index][pos[index]]]
+            pos[index] = pos[index] + 1
+            summed_volume = tmp_volume
+
+        # Update position of index processor and current stage
+        prev_stage = current_stage
+
+        # Calculate new stage
+        new_stage = prev_stage + (volume - summed_volume) / summed_areas
+
+
+        # Update own depth
+        stages[stages_order[0:pos]] = new_stage
+
+
+
+        self.set_stages(stages)
+
+        stages = self.get_stages()
+        stages_order = stages.argsort()
+
 
     def set_depths_evenly(self,volume):

trunk/anuga_core/source/anuga/structures/inlet_operator.py

@@ -36 +36 @@
 
 
+        #print self.Q
+
         self.enquiry_point = 0.5*(self.line[0] + self.line[1])
         self.outward_vector = self.line
@@ -59 +61 @@
         Q2 = self.update_Q(t + timestep)
 
+
+        #print Q1,Q2
         Q = 0.5*(Q1+Q2)
         volume = Q*timestep
+
+        #print volume
         
         #print Q, volume
+
+        # store last discharge
+        self.applied_Q = Q
 
         msg =  'Requesting too much water to be removed from an inlet! \n'
         msg += 'current_water_volume = %5.2e Increment volume = %5.2e' % (current_volume, volume)
-        assert current_volume + volume >= 0.0, msg
-
-        # store last discharge
-        self.applied_Q = Q
-
+        import warnings
+        if current_volume + volume < 0.0:
+            #warnings.warn(msg)
+            volume = -current_volume
+            self.applied_Q = volume/timestep
+
+
+        #print 'applied_Q', self.applied_Q
+        
         # Distribute volume so as to obtain flat surface
         self.inlet.set_stages_evenly(volume)
@@ -76 +89 @@
         # Distribute volume evenly over all cells
         #self.inlet.set_depths_evenly(volume)
-        
+
+
+
     def update_Q(self, t):
         """Allowing local modifications of Q
@@ -84 +99 @@
         if callable(self.Q):
             try:
-                Q = self.Q(t)[0]
+                Q = self.Q(t)
             except Modeltime_too_early, e:
                 raise Modeltime_too_early(e)

trunk/anuga_core/source/anuga_parallel/run_sequential_dist_distribute_merimbula.py

@@ -35 +35 @@
 
 
-from anuga_parallel.sequential_distribute import sequential_distribute
+from anuga_parallel.sequential_distribute import sequential_distribute_dump
 
 
@@ -95 +95 @@
 
 if verbose: print 'DISTRIBUTING DOMAIN'
-sequential_distribute(domain, 4, verbose=True)
+sequential_distribute_dump(domain, 4, verbose=True)
 
 

trunk/anuga_core/source/anuga_parallel/run_sequential_dist_evolve_merimbula.py

@@ -37 +37 @@
 from anuga_parallel import distribute, myid, numprocs, finalize, barrier
 
+from anuga_parallel.sequential_distribute import sequential_distribute_load
 
 #--------------------------------------------------------------------------
 # Setup parameters
 #--------------------------------------------------------------------------
-
-#mesh_filename = "merimbula_10785_1.tsh" ; x0 = 756000.0 ; x1 = 756500.0
-mesh_filename = "merimbula_17156.tsh"   ; x0 = 756000.0 ; x1 = 756500.0
-#mesh_filename = "merimbula_43200_1.tsh"   ; x0 = 756000.0 ; x1 = 756500.0
-#mesh_filename = "test-100.tsh" ; x0 = 0.25 ; x1 = 0.5
-#mesh_filename = "test-20.tsh" ; x0 = 250.0 ; x1 = 350.0
-yieldstep = 50
+yieldstep = 10
 finaltime = 1500
 verbose = True
 
-#--------------------------------------------------------------------------
-# Setup procedures
-#--------------------------------------------------------------------------
-class Set_Stage:
-    """Set an initial condition with constant water height, for x0<x<x1
-    """
-
-    def __init__(self, x0=0.25, x1=0.5, h=1.0):
-        self.x0 = x0
-        self.x1 = x1
-        self.h  = h
-
-    def __call__(self, x, y):
-        return self.h*((x>self.x0)&(x<self.x1))+1.0
-
-
-class Set_Elevation:
-    """Set an elevation
-    """
-
-    def __init__(self, h=1.0):
-        self.x0 = x0
-        self.x1 = x1
-        self.h  = h
-
-    def __call__(self, x, y):
-        return x/self.h
-    
-
-##--------------------------------------------------------------------------
-## Setup Domain only on processor 0
-##--------------------------------------------------------------------------
-#if myid == 0:
-#    domain = create_domain_from_file(mesh_filename)
-#    domain.set_quantity('stage', Set_Stage(x0, x1, 2.0))
-#    #domain.set_datadir('.')
-#    domain.set_name('merimbula_new')
-#    domain.set_store(True)
-#    #domain.set_quantity('elevation', Set_Elevation(500.0))
-#
-#    #print domain.statistics()
-#    #print domain.get_extent()
-#    #print domain.get_extent(absolute=True)
-#    #print domain.geo_reference
-#else:
-#    domain = None
-#
-##--------------------------------------------------------------------------
-## Distribute sequential domain on processor 0 to other processors
-##--------------------------------------------------------------------------
-#
-#if myid == 0 and verbose: print 'DISTRIBUTING DOMAIN'
-#domain = distribute(domain)
 
 
 #-------------------------------------------------------------------------------
-# Read in cPickled parallel domains
+# Read in cPickled distribute mesh data
 #-------------------------------------------------------------------------------
-
-
-import cPickle
-pickle_name = 'merimbula_new_P%g_%g.pickle'% (numprocs,myid)
-f = file(pickle_name, 'rb')
-domain = cPickle.load(f)
-f.close()
+domain = sequential_distribute_load(filename='merimbula_new', verbose = True)
 
 #--------------------------------------------------------------------------
@@ -122 +58 @@
 #--------------------------------------------------------------------------
 
-domain.set_flow_algorithm('2_0')
+domain.set_flow_algorithm('tsunami')
 
-#domain.smooth = False
-#domain.set_default_order(2)
-#domain.set_timestepping_method('rk2')
-#domain.set_CFL(0.7)
-#domain.set_beta(1.5)
-
-
-#for p in range(numprocs):
-#    if myid == p:
-#        print 'P%d'%p
-#        print domain.get_extent()
-#        print domain.get_extent(absolute=True)
-#        print domain.geo_reference
-#        print domain.s2p_map
-#        print domain.p2s_map
-#        print domain.tri_l2g
-#        print domain.node_l2g
-#    else:
-#        pass
-#
-#    barrier()
 
 
@@ -185 +100 @@
 # Merge the individual sww files into one file
 #--------------------------------------------------
-domain.sww_merge(delete_old=False)
+
+
+
 
 finalize()

trunk/anuga_core/validation_tests/Tests/Benchmarks/Okushiri/create_okushiri.py

@@ -28 +28 @@
 
 
-base_resolution = 1 # Use this to coarsen or refine entire mesh. 
+base_resolution = 0.25 # Use this to coarsen or refine entire mesh.
 
 # Basic geometry and bounding polygon
@@ -165 +165 @@
     
     if elevation_in_mesh is True:
+        from anuga.fit_interpolate.fit import fit_to_mesh_file
         fit_to_mesh_file(project.mesh_filename, project.bathymetry_filename,
-                         project.mesh_filename)
+                         project.mesh_filename, verbose = True)
 
 
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    create_mesh()
+    create_mesh(elevation_in_mesh=True)
 
 

trunk/anuga_core/validation_tests/Tests/Benchmarks/Okushiri/run_okushiri.py

@@ -101 +101 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    main()
+    main(elevation_in_mesh=True)