Changeset 8189

Timestamp:

Jul 1, 2011, 2:52:36 PM (14 years ago)

Author:

steve

Message:

Added in Sudi's most recent 1d codes

Location:

trunk/anuga_work/development/2010-projects/anuga_1d

Files:

• avalanche-sudi (added)
• avalanche-sudi/comp_flux_ext_include_entropy.c (added)
• avalanche-sudi/config.py (added)
• avalanche-sudi/debris_avalanche_RUN_to_the_left.py (added)
• avalanche-sudi/domain.py (added)
• avalanche-sudi/generic_boundary_conditions.py (added)
• avalanche-sudi/parameters.py (added)
• avalanche-sudi/quantity.py (added)
• avalanche-sudi/sww_domain_shm.py (added)
• avalanche-sudi/util.py (added)
• avalanche-sudi/util_ext.h (added)
• base/limiters_python.py (modified) (2 diffs)
• base/quantity_ext.c (modified) (2 diffs)
• sww-sudi (added)
• sww-sudi/analytical_prescription.py (added)
• sww-sudi/bisect_function.py (added)
• sww-sudi/comp_flux_ext_wellbalanced.c (added)
• sww-sudi/config.py (added)
• sww-sudi/domain.py (added)
• sww-sudi/domain_sudi.py (added)
• sww-sudi/error.py (added)
• sww-sudi/gaussPivot.py (added)
• sww-sudi/generic_boundary_conditions.py (added)
• sww-sudi/parameter.py (added)
• sww-sudi/quantity.py (added)
• sww-sudi/rootsearch.py (added)
• sww-sudi/run-discrepancy.py (added)
• sww-sudi/run-program-numeric-sudi.py (added)
• sww-sudi/swap.py (added)
• sww-sudi/util.py (added)
• sww/run_parabolic_canal.py (modified) (1 diff)
• sww/sww_comp_flux_ext.c (modified) (1 diff)
• sww/sww_domain.py (modified) (2 diffs)
• test_all.py (modified) (1 diff)

trunk/anuga_work/development/2010-projects/anuga_1d/base/limiters_python.py

@@ -12 +12 @@
     from numpy import abs, where
 
-    phi = where((abs(a) < abs(b)) & (a*b > 0.0), a, 0.0)
-    phi = where((abs(b) < abs(a)) & (a*b > 0.0), b, phi)
+    phi = where((abs(a) < abs(b)) & (a*b >= 0.0), a, 0.0)
+    phi = where((abs(b) < abs(a)) & (a*b >= 0.0), b, phi)
 
     return phi
@@ -20 +20 @@
     from numpy import sign, abs, minimum, where
 
-    return where( (sign(a)*sign(b) > 0.0) & (sign(a)*sign(c)>0.0),
+    return where( (sign(a)*sign(b) >= 0.0) & (sign(a)*sign(c)>0.0),
         sign(a)*minimum(minimum(abs(a),abs(b)),abs(c)), 0.0 )
 

trunk/anuga_work/development/2010-projects/anuga_1d/base/quantity_ext.c

@@ -88 +88 @@
 
           phi = 0.0;
-          if ((fabs(a) < fabs(b)) & (a*b > 0.0 )) {
+          if ((fabs(a) < fabs(b)) & (a*b >= 0.0 )) {
               phi = a;
           }
-          if ((fabs(b) < fabs(a)) & (a*b > 0.0 )) {
+          if ((fabs(b) < fabs(a)) & (a*b >= 0.0 )) {
               phi = b;
           }
@@ -191 +191 @@
 
           phi = 0.0;
-          if ((sign(a)*sign(b) > 0.0) & (sign(a)*sign(c) > 0.0 )) {
+          if ((sign(a)*sign(b) > 0.0) & (sign(a)*sign(c) >= 0.0 )) {
               phi = sign(a)*min(theta*min(fabs(a),fabs(b)),fabs(c));
           }

trunk/anuga_work/development/2010-projects/anuga_1d/sww/run_parabolic_canal.py

@@ -1 +1 @@
 import os
 from math import sqrt, pi, sin, cos
-from shallow_water_vel_domain import *
-from Numeric import allclose, array, zeros, ones, Float, take, sqrt
-from config import g, epsilon
-
-def analytic_cannal(C,t):
-    N = len(C)
-    u = zeros(N,Float)    ## water velocity
-    h = zeros(N,Float)    ## water depth
-    x = C
-    g = 9.81
-    ## Define Basin Bathymetry
-    z_b = zeros(N,Float) ## elevation of basin
-    z = zeros(N,Float)   ## elevation of water surface
-    z_infty = 10.0       ## max equilibrium water depth at lowest point.
-    L_x = 2500.0         ## width of channel
-    A0 = 0.5*L_x                  ## determines amplitudes of oscillations
-    omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
-    for i in range(N):
-        z_b[i] = z_infty*(x[i]**2/L_x**2)
-        u[i] = -A0*omega*sin(omega*t)
-        z[i] = z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x[i]/L_x-0.5*A0/(L_x)*cos(omega*t))
-    h = z-z_b
-    T = 2.0*pi/omega
-    return u,h,z,z_b, T
+import numpy
+import time
 
 
-L_x = 2500.0     # Length of channel (m)
-N = 1000         # Number of compuational cells
-cell_len = 4*L_x/N   # Origin = 0.0
+from anuga_1d.sww.sww_domain import *
+from anuga_1d.config import g, epsilon
+from anuga_1d.base.generic_mesh import uniform_mesh
 
-points = zeros(N+1,Float)
-for i in range(N+1):
-        points[i] = -2*L_x +i*cell_len
+#===============================================================================
+# setup problem
+#===============================================================================
+
+
+z_infty = 5.0       ## max equilibrium water depth at lowest point.
+L_x = 2500.0         ## width of channel
+A0 = 0.5*L_x                  ## determines amplitudes of oscillations
+omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
+
+def analytic_canal(x,t):
+    u = numpy.zeros_like(x)    ## water velocity
+    h = numpy.zeros_like(x)    ## water depth
+
+    ## Define Basin Bathymetry
+    z = numpy.zeros_like(x) ## elevation of basin
+    w = numpy.zeros_like(x)   ## elevation of water surface
+
+    z[:] = z_infty*(x**2/L_x**2)
+    u[:] = -A0*omega*sin(omega*t)
+    w[:] = numpy.maximum(z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x/L_x-0.5*A0/(L_x)*cos(omega*t)),z)
+    h[:] = numpy.maximum(w-z, 0.0)
+
+    T = 2.0*pi/omega
+
+    return u,h,w,z, T
+
 
 def stage(x):
-    z_infty = 10.0       ## max equilibrium water depth at lowest point.
-    L_x = 2500.0         ## width of channel
-    A0 = 0.5*L_x                  ## determines amplitudes of oscillations
-    omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
     t=0.0
-    y = zeros(len(x),Float)
-    for i in range(len(x)):
-        y[i] = z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x[i]/L_x-0.5*A0/(L_x)*cos(omega*t))
-        #y[i] = 12.0
-    return y
+    u,h,w,z,T = analytic_canal(x,t)
+
+    return w
 
 def elevation(x):
-    N = len(x)
-    z_infty = 10.0
-    z = zeros(N,Float)
-    L_x = 2500.0
-    A0 = 0.5*L_x
-    omega = sqrt(2*g*z_infty)/L_x
-    for i in range(N):
-        z[i] = z_infty*(x[i]**2/L_x**2)
+    t=0.0
+    u,h,w,z,T = analytic_canal(x,t)
+
     return z
 
+
 def height(x):
-    z_infty = 10.0       ## max equilibrium water depth at lowest point.
-    L_x = 2500.0         ## width of channel
-    A0 = 0.5*L_x                  ## determines amplitudes of oscillations
-    omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
     t=0.0
-    y = zeros(len(x),Float)
-    for i in range(len(x)):
-        y[i] = max(z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x[i]/L_x-0.5*A0/(L_x)*cos(omega*t))-z_infty*(x[i]**2/L_x**2),0.0)
-    return y
+    u,h,w,z,T = analytic_canal(x,t)
 
-domain = Domain(points)
-domain.order = 2
-domain.set_timestepping_method('euler')
+    return h
+
+
+#===============================================================================
+finaltime = 1000.0
+yieldstep = 10.0
+
+
+
+N = 100
+print "Evaluating domain with %d cells" %N
+domain = Domain(*uniform_mesh(N, x_0 = -2.0*L_x, x_1 = 2.0*L_x))
+
+domain.set_spatial_order(2)
+domain.set_timestepping_method('rk2')
 domain.set_CFL(1.0)
-domain.beta = 1.0
-domain.set_limiter("vanleer")
+domain.set_limiter("minmod_kurganov")
 
-    
+domain.set_beta(1.0)
+
 domain.set_quantity('stage', stage)
 domain.set_quantity('elevation',elevation)
-domain.set_boundary({'exterior': Reflective_boundary(domain)})
- 
-X = domain.vertices
-u,h,z,z_b,T = analytic_cannal(X.flat,domain.time)
-print 'T = ',T
 
-yieldstep = finaltime = 10.0 #T/2.0
-StageQ = domain.quantities['stage'].vertex_values
-XmomQ = domain.quantities['xmomentum'].vertex_values
+domain.quantities['elevation'].extrapolate_second_order()
+domain.quantities['stage'].extrapolate_second_order()
 
-import time
+Br = Reflective_boundary(domain)
+
+domain.set_boundary({'left': Br, 'right' : Br})
+
+#domain.h0=0.0001
+
 t0 = time.time()
+
 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
     domain.write_time()
-    print "integral", domain.quantities['stage'].get_integral()
-    if t>0.0:
-        
-        x = X.flat
-        print 'domain.time=',domain.time
-        
-        w_V = domain.quantities['stage'].vertex_values.flat
-        uh_V = domain.quantities['xmomentum'].vertex_values.flat
-        u_V = domain.quantities['velocity'].vertex_values.flat
-        h_V = domain.quantities['height'].vertex_values.flat
-        
-        u,h,z,z_b,T = analytic_cannal(x,domain.time)
-        w = z
-        for k in range(len(h)):
-            if h[k] < 0.0:
-                h[k] = 0.0
-                w[k] = z_b[k]
-        
-        from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
-        hold(False)
 
-        #print 'size X',X.shape
-        #print 'size w',w.shape
-        signh = (h_V>0)
+print 'That took %.2f seconds' %(time.time()-t0)
 
-        plot1 = subplot(311)
-        #plot(x,w,  x,w_V,  x,z_b)
-        plot(x,signh)
-        plot1.set_xlim([-6000,6000])
-        xlabel('Position')
-        ylabel('Stage')
-        legend(('Analytic Solution', 'Numerical Solution', 'Bed'),
-               'upper center', shadow=True)
 
-        
-        plot2 = subplot(312)
-        plot(x,u*h,x,uh_V)
-        #plot2.set_ylim([-1,25])
-        xlabel('Position')
-        ylabel('Xmomentum')
+N = float(N)
+HeightC    = domain.quantities['height'].centroid_values
+StageC     = domain.quantities['stage'].centroid_values
+BedC       = domain.quantities['elevation'].centroid_values
+C          = domain.centroids
 
-        print u_V
-        
-        plot3 = subplot(313)
-        plot(x,u, x,u_V)
-        #plot2.set_ylim([-1,25])
-        xlabel('Position')
-        ylabel('Velocity')
-        show()
-raw_input("Press the return key!")
+HeightV    = domain.quantities['height'].vertex_values
+StageV     = domain.quantities['stage'].vertex_values
+BedV       = domain.quantities['elevation'].vertex_values
+VelocityV  = domain.quantities['velocity'].vertex_values
+X          = domain.vertices
+
+
+import pylab
+
+pylab.hold(False)
+
+plot1 = pylab.subplot(311)
+
+pylab.plot(X.flat,BedV.flat,X.flat,StageV.flat)
+
+plot1.set_ylim([-1,40])
+
+pylab.xlabel('Position')
+pylab.ylabel('Stage')
+pylab.legend(('Analytical Solution', 'Numerical Solution'),
+           'upper right', shadow=True)
+
+
+plot2 = pylab.subplot(312)
+
+pylab.plot(X.flat,HeightV.flat)
+
+plot2.set_ylim([-1,10])
+
+pylab.xlabel('Position')
+pylab.ylabel('Height')
+
+plot3 = pylab.subplot(313)
+
+pylab.plot(X.flat,VelocityV.flat)
+plot3.set_ylim([-15,15])
+
+pylab.xlabel('Position')
+pylab.ylabel('Velocity')
+
+pylab.show()
+
+
+
+
+
+

trunk/anuga_work/development/2010-projects/anuga_1d/sww/sww_comp_flux_ext.c

@@ -3 +3 @@
 #include "math.h"
 #include <stdio.h>
+
 const double pi = 3.14159265358979;
 

trunk/anuga_work/development/2010-projects/anuga_1d/sww/sww_domain.py

@@ -237 +237 @@
     h_C[:] = w_C - z_C
     u_C[:]  = uh_C/(h_C + h0/h_C)
-        
+
+    #print 'domain.order', domain.order
+    
     for name in [ 'velocity', 'stage' ]:
         Q = domain.quantities[name]
@@ -258 +260 @@
 
     h_V[:,:]    = stage_V - bed_V
+
+
+    #print 'any' ,numpy.any( h_V[:,0] < 0.0)
+    #print 'any' ,numpy.any( h_V[:,1] < 0.0)
+
+
+    h_0 = numpy.where(h_V[:,0] < 0.0, 0.0, h_V[:,0])
+    h_1 = numpy.where(h_V[:,0] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,1])
+
+    h_V[:,0] = h_0
+    h_V[:,1] = h_1
+
+
+    h_0 = numpy.where(h_V[:,1] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,0])
+    h_1 = numpy.where(h_V[:,1] < 0.0, 0.0, h_V[:,1])
+
+    h_V[:,0] = h_0
+    h_V[:,1] = h_1
+    
+    #print 'any' ,numpy.any( h_V[:,0] < 0.0)
+    #print 'any' ,numpy.any( h_V[:,1] < 0.0)
+
+    #h00 = 1e-12
+    #print h00
+
+    h_V[:,:] = numpy.where (h_V <= h0, 0.0, h_V)
+    u_V[:,:] = numpy.where (h_V <= h0, 0.0, u_V)
 
 #    # protect from edge values going negative

trunk/anuga_work/development/2010-projects/anuga_1d/test_all.py

@@ -23 +23 @@
 
 # Directories that should not be searched for test files.
-exclude_dirs = ['channel',  '.svn'] 
+exclude_dirs = [ '.svn']