Changeset 8158

Timestamp:

Mar 20, 2011, 11:32:14 AM (14 years ago)

Author:

steve

Message:

Using brent solver to solve for intermediate state in dam break calc

File:

• trunk/anuga_validation/validation_tests/dam_break_analytical.py (modified) (3 diffs)

trunk/anuga_validation/validation_tests/dam_break_analytical.py

@@ -21 +21 @@
 #       MA 02110-1301, USA.
 
-
     
-def wu(h2):
+def wu(h2, h0=1.0, h1=10.0):
     """ The intermediate state in dambreak problem
     must satisfy wu(h2) = 0.0
@@ -29 +28 @@
     Taken from formula (3.14) from Zoppou and Roberts SWW notes
     """
-    
+
     h21 = h2/h1
     h01 = h0/h1
@@ -36 +35 @@
 
 
-h0 = 1.0
-h1 = 2.0
+def h2(h0=1.0, h1=10.0):
 
-h21 = 1.0
-h01 = 0.5
+    assert h1 >  0.0
+    assert h0 <=  h1
+    assert h0 >= 0.0
 
-print h21**3 - 9.0*(h21**2)*h01 + 16.0*(h21**(1.5))*h01 - h21*h01*(h01+8.0) + h01**3
-
-h21 = 0.5
-h01 = 0.5
-
-print  - 8.0*h21**3 + 16.0*h21**(2.5) -8*h21**2 
+    from scipy.optimize import brentq
+    return brentq(wu, h0, h1, args=(h0,h1))
 
 
-n = 10
-import pylab as p
+h1 = 1.0
 
-import numpy
-
-x = h0 + numpy.arange(0,n+1)/float(n)*(h1-h0)
-
-new_wu = numpy.vectorize(wu)
-y = new_wu(x)
-
-print x
-print y
+for i in range(0,11):
+    h0 = float(i)/10.0*h1
+    print 'h2(%g,%g) = %g' % (h0, h1, h2(h0))
 
 
-print wu(h0), wu(h1)
 
-from scipy.optimize import brentq
 
-h2 = brentq(wu, h0,h1,disp=True)
 
-#print  h2, wu(h2) 
-