Changeset 9158

Timestamp:

Jun 15, 2014, 5:22:06 PM (10 years ago)

Author:

steve

Message:

Combined plot_result scripts

Location:

trunk/anuga_core/validation_tests/analytical_exact/paraboloid_basin

Files:

• numerical_paraboloid_basin.py (modified) (4 diffs)
• plot_results.py (added)
• produce_results.py (modified) (1 diff)

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

@@ -19 +19 @@
 from time import localtime, strftime, gmtime
 from anuga.geometry.polygon import inside_polygon, is_inside_triangle
-
+from anuga import myid, finalize, distribute
 
 #-------------------------------------------------------------------------------
@@ -32 +32 @@
 
 
+args = anuga.get_args()
+alg = args.alg
+verbose = args.v
+
 #-------------------------------------------------------------------------------
 # Domain
@@ -41 +45 @@
 origin = (-4000.0, -4000.0)
 
-points, elements, boundary = anuga.rectangular_cross(m, n, lenx, leny, origin)
-domain = anuga.Domain(points, elements, boundary)
 
 
-#-------------------------------------------------------------------------------
-# Provide file name for storing output
-#-------------------------------------------------------------------------------
-domain.store = True #False
-domain.format = 'sww'
-domain.set_name(output_file)                
-domain.set_datadir(output_dir)
+#===============================================================================
+# Create Sequential Domain
+#===============================================================================
+if myid == 0:
+    points, elements, boundary = anuga.rectangular_cross(m, n, lenx, leny, origin)
+    domain = anuga.Domain(points, elements, boundary)
 
 
-#-------------------------------------------------------------------------------
-# Order of scheme
-# Good compromise between
-# limiting and CFL
-# CUT THIS FOR AUTOMATED TESTS
-#-------------------------------------------------------------------------------
-#domain.set_default_order(2)
-#domain.set_timestepping_method(2)
-#domain.set_beta(1.5)
-##domain.set_CFL(0.6)
-#domain.smooth = True
+    domain.set_name(output_file)                
+    domain.set_datadir(output_dir)
+    
+    domain.set_flow_algorithm(alg)
 
+    
+    #-------------------------------------------------------------------------------
+    # Initial conditions
+    #-------------------------------------------------------------------------------
+    t = 0.0
+    D0 = 1000.
+    L = 2500.
+    R0 = 2000.
+    
+    A = (L**4 - R0**4)/(L**4 + R0**4)
+    omega = 2./L*sqrt(2.*g*D0)
+    T = pi/omega
+    
+    #------------------
+    # Set bed elevation
+    def bed_elevation(x,y):
+        n = x.shape[0]
+        z = 0*x
+        for i in range(n):
+            r = sqrt(x[i]*x[i] + y[i]*y[i])
+            z[i] = -D0*(1.-r*r/L/L)
+        return z
+    domain.set_quantity('elevation', bed_elevation)
+    
+    #----------------------------
+    # Set the initial water level
+    def stage_init(x,y):
+        z = bed_elevation(x,y)
+        n = x.shape[0]
+        w = 0*x
+        for i in range(n):
+            r = sqrt(x[i]*x[i] + y[i]*y[i])
+            w[i] = D0*((sqrt(1-A*A))/(1.-A*cos(omega*t))
+                    -1.-r*r/L/L*((1.-A*A)/((1.-A*cos(omega*t))**2)-1.))
+        if w[i] < z[i]:
+            w[i] = z[i]
+        return w
+    domain.set_quantity('stage', stage_init)
 
-#-------------------------------------------------------------------------------
-# Decide which quantities are to be stored at each timestep
-#-------------------------------------------------------------------------------
-#domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
+else:
+    
+    domain = None
+ 
+#============================================================================
+# Create Parallel Domain
+#============================================================================
+domain = distribute(domain)
 
-
-#-------------------------------------------------------------------------------
-# Reduction operation for get_vertex_values
-from anuga.utilities.numerical_tools import mean
-domain.reduction = mean #domain.reduction = min  #Looks better near steep slopes
-
-
-#-------------------------------------------------------------------------------
-# Initial condition
-#-------------------------------------------------------------------------------
-t = 0.0
-D0 = 1000.
-L = 2500.
-R0 = 2000.
-
-A = (L**4 - R0**4)/(L**4 + R0**4)
-omega = 2./L*sqrt(2.*g*D0)
-T = pi/omega
-
-#------------------
-# Set bed elevation
-def bed_elevation(x,y):
-    n = x.shape[0]
-    z = 0*x
-    for i in range(n):
-        r = sqrt(x[i]*x[i] + y[i]*y[i])
-        z[i] = -D0*(1.-r*r/L/L)
-    return z
-domain.set_quantity('elevation', bed_elevation)
-
-#----------------------------
-# Set the initial water level
-def stage_init(x,y):
-    z = bed_elevation(x,y)
-    n = x.shape[0]
-    w = 0*x
-    for i in range(n):
-        r = sqrt(x[i]*x[i] + y[i]*y[i])
-        w[i] = D0*((sqrt(1-A*A))/(1.-A*cos(omega*t))
-                -1.-r*r/L/L*((1.-A*A)/((1.-A*cos(omega*t))**2)-1.))
-    if w[i] < z[i]:
-        w[i] = z[i]
-    return w
-domain.set_quantity('stage', stage_init)
 
 #---------
@@ -124 +116 @@
 domain.set_boundary({'left': D, 'right': D, 'top': D, 'bottom': D})
 
-#---------------------------------------------
-# 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
-    
-#t = domain.triangles[n_point]
-#print 't = ', t
-#tri = domain.get_vertex_coordinates(n)
-
-filename=domain.get_name()
-file = open(filename,'w')
-
-#----------
-# Evolution
-import time
-t0 = time.time()
-
-time_array = []
-stage_array = []
-Stage     = domain.quantities['stage']
-Xmomentum = domain.quantities['xmomentum']
-Ymomentum = domain.quantities['ymomentum']
-
-#interactive_visualisation = True
-#===============================================================================
-##if interactive_visualisation:
-##    from anuga.visualiser import RealtimeVisualiser
-##    vis = RealtimeVisualiser(domain)
-##    vis.render_quantity_height("stage", zScale =1000, dynamic=True)
-##    vis.colour_height_quantity('stage', (1.0, 0.5, 0.5))
-##    vis.start()
-#===============================================================================
-
 
 #------------------------------------------------------------------------------
 # 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()
 
+import time
+t0 = time.time()
 for t in domain.evolve(yieldstep = 1.0, finaltime = 200.0):
-    domain.write_time()
+    if myid == 0 and verbose : domain.write_time()
 
-    #tri_array = asarray(tri)
-    #t_array = asarray([[0,1,2]])
-    #interp = Interpolation(tri_array,t_array,[points])
+if myid==0 and verbose: print 'That took %s secs' % str(time.time()- t0)
 
-    stage     = Stage.get_values(location='centroids',indices=[n_point])
-    xmomentum = Xmomentum.get_values(location='centroids',indices=[n_point])
-    ymomentum = Ymomentum.get_values(location='centroids',indices=[n_point])
-    #print '%10.6f   %10.6f  %10.6f   %10.6f\n'%(t,stage,xmomentum,ymomentum)
-    file.write( '%10.6f   %10.6f  %10.6f   %10.6f\n'%(t,stage,xmomentum,ymomentum) )
+domain.sww_merge(delete_old=True)
 
-    #time_array.append(t)
-    #stage_array.append(stage)
-
-    #if interactive_visualisation:
-    #    vis.update()
-"""
-file.close()
-print 'That took %.2f seconds' %(time.time()-t0)
-from pylab import *
-#ion()
-hold(False)
-plot(time_array, stage_array, 'r.-')
-#title('Gauge %s' %name)
-xlabel('time (s)')
-ylabel('stage (m)')
-#legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
-#savefig(name, dpi = 300)
-
-#raw_input('Next')
-show()
-"""
-
-#if interactive_visualisation:
-#    vis.evolveFinished()
+finalize()

trunk/anuga_core/validation_tests/analytical_exact/paraboloid_basin/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()
+
+produce_report('numerical_paraboloid_basin.py', args=args)
 
 
-# Setup the python scripts which produce the output for this
-# validation test
-def build():
-    run_validation_script('numerical_paraboloid_basin.py')
-    run_validation_script('plot_results_cross_section.py')
-    run_validation_script('plot_results_origin_wrt_time.py')    
-    typeset_report()
-
-def clean():
-    autoclean()
-
-main()