Changeset 8582 for trunk/anuga_core

Timestamp:

Sep 18, 2012, 3:42:23 PM (13 years ago)

Author:

steve

Message:

Moved Gareth's tsunami algorithm over to main trunk

Location:

trunk/anuga_core

Files:

• source/anuga/abstract_2d_finite_volumes/generic_domain.py (modified) (1 diff)
• source/anuga/shallow_water/shallow_water_domain.py (modified) (3 diffs)
• source/anuga/shallow_water/swb2_domain_ext.c (modified) (6 diffs)
• source/anuga/shallow_water/test_swb2_domain.py (added)
• source/anuga_parallel/parallel_api.py (modified) (4 diffs)
• source/anuga_parallel/run_parallel_sw_rectangular_cross.py (modified) (4 diffs)
• validation_tests/Tests/Simple/runup_sinusoid/runup_sinusoid.py (modified) (3 diffs)
• validation_tests/Tests/Simple/runup_sinusoid/runup_sinusoidplot.py (modified) (1 diff)
• validation_tests/parameters.py (modified) (1 diff)

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_domain.py

@@ -141 +141 @@
         self.geo_reference = self.mesh.geo_reference
 
+        self.verbose = verbose
 
         if verbose: log.critical('Domain: Expose quantity names and types')

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

@@ -344 +344 @@
 
         # We demand that vertex values are stored uniquely
-        self.set_store_vertices_smoothly(False)
+        #self.set_store_vertices_smoothly(False)
 
         self.maximum_allowed_speed=0.0
+        #self.minimum_allowed_height=0.01
+        
         #self.forcing_terms.append(manning_friction_explicit)
         #self.forcing_terms.remove(manning_friction_implicit)
 
-        print '##########################################################################'
-        print '#'
-        print '# Using anuga_tsunami solver in anuga_work/development/gareth/experimental/anuga_tsunami/'
-        print "#"
-        print "# This solver is experimental. Here are some tips on using it"
-        print "#"
-        print "# 1) When plotting outputs, I strongly suggest you examine centroid values, not vertex values"
-        print "# , as the latter can be completely misleading near strong gradients in the flow. "
-        print "# There is a plot_util.py script in anuga_core/utilities/ which might help you extract"
-        print "# quantities at centroid values from sww files."
-        print "# Note that to accuractely compute centroid values from sww files, the files need to store "
-        print "# vertices uniquely. This makes for large sww files (3x), but is the price to pay for the right answer"
-        print "# (unless we alter IO to allow centroids to be written to sww files, which would then affect"
-        print "# ANUGA viewer as well -- I expect this would be lots of work)"
-        print "#"
-        print "# 2) In field scale applications (where the depth is typically > 1m), I suggest you set"
-        print "# domain.minimum_allowed_height=0.01 (the default is 1.0e-3). "
-        print "#"
-        print "# 3) This solver is not expected to perform well in problems with very"
-        print "# shallow water flowing down steep slopes (such that the stage_centroid_value "
-        print "# is less than the maximum bed_edge_value on a given triangle). However, analytical tests"
-        print "# suggest it can do typical wetting/drying situations very well (parabolic oscillations test case) "
-        print "#"
-        print "# 4) This solver allows the stage_centroid_value to drop to slightly below the minimum bed_vertex_value"
-        print "# on it's triangle. In other ANUGA versions (e.g. 1.2.1), the limit would be the"
-        print "# bed_centroid_value. This means that triangles store slightly more water than they are"
-        print "# typically interpreted to store, which might have significance in some applications."
-        print "#"
-        print "# You will probably be able to tell this is causing you problems by convergence testing"
-        print "#"
-        print '# 5) Note that many options in config.py have been overridden by the solver -- I have '
-        print '# deliberately attempted to get the solver to perform well with consistent values of '
-        print '# these parameters -- so I would advise against changing them unless you at least check that '
-        print '# it really does improve things'
-        print '#'
-        print '##########################################################################'
+#        print '##########################################################################'
+#        print '#'
+#        print '# Using anuga_tsunami solver in anuga_work/development/gareth/experimental/anuga_tsunami/'
+#        print "#"
+#        print "# This solver is experimental. Here are some tips on using it"
+#        print "#"
+#        print "# 1) When plotting outputs, I strongly suggest you examine centroid values, not vertex values"
+#        print "# , as the latter can be completely misleading near strong gradients in the flow. "
+#        print "# There is a plot_util.py script in anuga_core/utilities/ which might help you extract"
+#        print "# quantities at centroid values from sww files."
+#        print "# Note that to accuractely compute centroid values from sww files, the files need to store "
+#        print "# vertices uniquely. This makes for large sww files (3x), but is the price to pay for the right answer"
+#        print "# (unless we alter IO to allow centroids to be written to sww files, which would then affect"
+#        print "# ANUGA viewer as well -- I expect this would be lots of work)"
+#        print "#"
+#        print "# 2) In field scale applications (where the depth is typically > 1m), I suggest you set"
+#        print "# domain.minimum_allowed_height=0.01 (the default is 1.0e-3). "
+#        print "#"
+#        print "# 3) This solver is not expected to perform well in problems with very"
+#        print "# shallow water flowing down steep slopes (such that the stage_centroid_value "
+#        print "# is less than the maximum bed_edge_value on a given triangle). However, analytical tests"
+#        print "# suggest it can do typical wetting/drying situations very well (parabolic oscillations test case) "
+#        print "#"
+#        print "# 4) This solver allows the stage_centroid_value to drop to slightly below the minimum bed_vertex_value"
+#        print "# on it's triangle. In other ANUGA versions (e.g. 1.2.1), the limit would be the"
+#        print "# bed_centroid_value. This means that triangles store slightly more water than they are"
+#        print "# typically interpreted to store, which might have significance in some applications."
+#        print "#"
+#        print "# You will probably be able to tell this is causing you problems by convergence testing"
+#        print "#"
+#        print '# 5) Note that many options in config.py have been overridden by the solver -- I have '
+#        print '# deliberately attempted to get the solver to perform well with consistent values of '
+#        print '# these parameters -- so I would advise against changing them unless you at least check that '
+#        print '# it really does improve things'
+#        print '#'
+#        print '##########################################################################'
 
 
@@ -1123 +1125 @@
             areas = self.areas
 
-            protect(self.minimum_allowed_height, self.maximum_allowed_speed,
+            mass_error = protect(self.minimum_allowed_height, self.maximum_allowed_speed,
                 self.epsilon, wc, wv, zc,zv, xmomc, ymomc, areas)
+
+            if mass_error > 0.0 and self.verbose :
+                print 'Cumulative mass protection: %g m^3 '% mass_error
 
 
@@ -1278 +1283 @@
             areas = self.areas
 
-            protect(self.minimum_allowed_height, self.maximum_allowed_speed,
+            mass_error = protect(self.minimum_allowed_height, self.maximum_allowed_speed,
                 self.epsilon, wc, wv, zc,zv, xmomc, ymomc, areas)
+
+            print 'Cumulative mass protection: python'+str(mass_error)+'m^3 '
+
 
         else:

trunk/anuga_core/source/anuga/shallow_water/swb2_domain_ext.c

@@ -558 +558 @@
 
 // Protect against the water elevation falling below the triangle bed
-int _protect(int N,
+double _protect(int N,
          double minimum_allowed_height,
          double maximum_allowed_speed,
@@ -573 +573 @@
   double hc, bmin, bmax;
   double u, v, reduced_speed;
-  static double mass_error=0.; 
+  double mass_error=0.; 
   // This acts like minimum_allowed height, but scales with the vertical
   // distance between the bed_centroid_value and the max bed_edge_value of
@@ -626 +626 @@
     }
 
+/*
   if(mass_added==1){
      printf("Cumulative mass protection: %f m^3 \n", mass_error);
   }
-  return 0;
+*/
+  return mass_error;
 }
 
@@ -1734 +1736 @@
   int N;
   double minimum_allowed_height, maximum_allowed_speed, epsilon;
+  double mass_error;
 
   // Convert Python arguments to C
@@ -1747 +1750 @@
   N = wc -> dimensions[0];
 
-  _protect(N,
+  mass_error = _protect(N,
        minimum_allowed_height,
        maximum_allowed_speed,
@@ -1759 +1762 @@
        (double*) areas -> data);
 
-  return Py_BuildValue("");
+  return Py_BuildValue("d", mass_error);
 }
 

trunk/anuga_core/source/anuga_parallel/parallel_api.py

@@ -68 +68 @@
         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
@@ -77 +79 @@
             # FIXME SR: Creates cPickle dump
             send((domain_name, domain_dir, domain_store, \
-                  domain_minimum_storable_height, georef, \
+                  domain_minimum_storable_height, domain_flow_algorithm, \
+                  domain_minimum_allowed_height, georef, \
                   number_of_global_triangles, number_of_global_nodes), p)
     else:
@@ -83 +86 @@
 
         domain_name, domain_dir, domain_store, \
-                  domain_minimum_storable_height, \
-                  georef, number_of_global_triangles, \
+                  domain_minimum_storable_height, domain_flow_algorithm, \
+                  domain_minimum_allowed_height, georef, number_of_global_triangles, \
                   number_of_global_nodes = receive(0)
 
@@ -236 +239 @@
     domain.set_store(domain_store)
     domain.set_minimum_storable_height(domain_minimum_storable_height)
+    domain.set_minimum_allowed_height(domain_minimum_allowed_height)
+    domain.set_flow_algorithm(domain_flow_algorithm)
     domain.geo_reference = georef   
 

trunk/anuga_core/source/anuga_parallel/run_parallel_sw_rectangular_cross.py

@@ -21 +21 @@
 # Sequential interface
 #---------------------------
-from anuga import Domain
 from anuga import Transmissive_boundary, Reflective_boundary
 from anuga import rectangular_cross_domain
@@ -44 +43 @@
     #dx = dy = 0.00125
     #dx = dy  = 0.5
-    seq_domain = rectangular_cross_domain(int(length/dx), int(width/dy),
+    domain = rectangular_cross_domain(int(length/dx), int(width/dy),
                                               len1=length, len2=width, verbose=verbose)
 
-    seq_domain.set_store(False)
-    seq_domain.set_quantity('elevation', lambda x,y : -1.0-x )
-    seq_domain.set_quantity('stage', 1.0)
-    seq_domain.print_statistics()
+
+    domain.set_store(True)
+    domain.set_quantity('elevation', lambda x,y : -1.0-x )
+    domain.set_quantity('stage', 1.0)
+    domain.set_flow_algorithm('tsunami')
+    domain.print_statistics()
 else:
-    seq_domain = None
+    domain = None
 
 t1 = time.time()
@@ -65 +66 @@
 barrier()
 
-# setup parameters to test using different ghost_layer_widths
-parameters = dict(ghost_layer_width = 2)
-domain = distribute(seq_domain,verbose=verbose, parameters=parameters)
-
-del seq_domain
+#-------------------------------------------------------------------------
+# Distribute domain
+#-------------------------------------------------------------------------
+domain = distribute(domain,verbose=verbose)
 
 
@@ -121 +121 @@
 if myid == 0 : print 'after set quantity'
 
-# Set Evolve parameters
-domain.set_flow_algorithm('2_0')
-
-
-
-
-#domain.update_ghosts()
-
-
-#print 'after evolve parameters'
-
-
-#import pdb; pdb.set_trace()
-
 
 

trunk/anuga_core/validation_tests/Tests/Simple/runup_sinusoid/runup_sinusoid.py

@@ -9 +9 @@
 
 from math import sin, pi, exp
-from anuga.shallow_water.shallow_water_domain import Domain as Domain
-#from anuga.shallow_water_balanced2.swb2_domain import Domain as Domain
-#path.append('/home/gareth/storage/anuga_clean/anuga_jan12/trunk/anuga_work/development/gareth/balanced_basic')
-#from swb2_domain import *
-#from balanced_basic import *
-#from balanced_dev import *
+
+from anuga import Domain
+
 #---------
 #Setup computational domain
 #---------
-points, vertices, boundary = anuga.rectangular_cross(40,40, len1=1., len2=1.)
+domain = anuga.rectangular_cross_domain(40,40, len1=1., len2=1.)
 
-domain=Domain(points,vertices,boundary)    # Create Domain
-domain.set_name('runup_sinusoid_v2')                         # Output to file runup.sww
-domain.set_datadir('.')                          # Use current folder
-domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
-#domain.set_store_vertices_uniquely(True)
-#domain.minimum_allowed_height=0.001
 
 #------------------------------------------------------------------------------
@@ -35 +26 @@
 domain.set_flow_algorithm(alg)
 domain.set_CFL(cfl)
+
+
+domain.set_name('runup_sinusoid_v2')                         # Output to file runup.sww
+domain.set_datadir('.')                          # Use current folder
+domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
+#domain.set_store_vertices_uniquely(True)
+domain.set_minimum_allowed_height(0.01)
+
 
 
@@ -87 +86 @@
 #------------------------------
 
-for t in domain.evolve(yieldstep=0.1,finaltime=20.0):
+for t in domain.evolve(yieldstep=0.5,finaltime=20.0):
     print domain.timestepping_statistics()
     xx = domain.quantities['xmomentum'].centroid_values

trunk/anuga_core/validation_tests/Tests/Simple/runup_sinusoid/runup_sinusoidplot.py

@@ -30 +30 @@
         pyplot.draw()
         pyplot.plot( (0,1),(0,0), 'r' )
-        pyplot.title(str(i)+'/200') # : velocity does not converge to zero' )
+        pyplot.title(str(i)+'/40') # : velocity does not converge to zero' )
         pyplot.xlabel('x')
         pyplot.ylabel('Velocity (m/s)')

trunk/anuga_core/validation_tests/parameters.py

@@ -8 +8 @@
 
 
-alg = '2_0'
+alg = 'tsunami'
 cfl = 1.0