Changeset 4980

Timestamp:

Jan 30, 2008, 9:41:35 AM (17 years ago)

Author:

nick

Message:

used sqrt table to calculate elevation, to get a more accurate model.

File:

• anuga_validation/conical_island/run_circular.py (modified) (4 diffs)

anuga_validation/conical_island/run_circular.py

@@ -11 +11 @@
 from anuga.shallow_water import Dirichlet_boundary, Time_boundary
 from anuga.shallow_water.data_manager import get_maximum_inundation_data
-from anuga.abstract_2d_finite_volumes.util import file_function, sww2csv_gauges
+from anuga.abstract_2d_finite_volumes.util import file_function, sww2csv_gauges,square_root
 from anuga.pmesh.mesh_interface import create_mesh_from_regions
 from anuga.utilities.polygon import read_polygon, plot_polygons
 #from cmath import cos,sin,cosh,pi
 from math import cos,pi,sin,tan,sqrt
-from Numeric import array, zeros, Float, allclose
+from Numeric import array, zeros, Float, allclose,resize
+from anuga.shallow_water.data_manager import csv2dict
 
 #import project
@@ -98 +99 @@
 d=0.75
 res=.05
+res=1.
 
 #create polygons (circles) to have higher resolution
@@ -164 +166 @@
 domain.set_quantity('stage', water_height)
 
-def test():                
-    fileName = "test.csv"
-    file = open(fileName,"w")
-    file.write("easting,northing,elevation \n")
-    x = 12.96
-    y = 13.80
-    for i in range(0,30):
-        for j in range(0,25):
-            a = (x-i)**2
-            b = (y-j)**2
-            z = -sqrt((x-i)**2+(y-j)**2)
-            z = z*.25+(0.8975)
-            if z <0:
-                z=0
-            if z > .625:
-                z=0.625
-
-#            print 'x,y,f',i,j,a,b,z
-            file.write("%s, %s, %s\n" %(i, j, z))
-    file.close()
-
-test()
-
-domain.set_quantity('elevation',filename = "test.csv", alpha=0.1)
+#def square_root():
+#    fileName1 = "square_root_table.csv"
+#    file1 = open(fileName1,"w")
+##    file1.write("number,square_root \n")
+#    for i in range (0,37):
+#        for j in range (0,10):
+#            number = i+(j/10.)
+#            print "%s,%s, %s, %s" %(i,j, number, sqrt(number))
+#            
+#            file1.write("%s, %s\n" %(i, sqrt(number)))
+#    file1.close()
+#
+##    return sqrt(s)
+#
+#def test():                
+#    fileName = "test.csv"
+#    file = open(fileName,"w")
+#    file.write("easting,northing,elevation \n")
+#    
+#   # csv2dict()
+#    
+#    x = 12.96
+#    y = 13.80
+#    res=1
+#    for i in range(0,30*res):
+#        for j in range(0,25*res):
+##            a = (x-i)**2
+##            b = (y-j)**2
+#
+#          #  z = -sqrt((x-(i/res))**2+(y-(j/res))**2)
+#            du = (x-(i/res))**2+(y-(j/res))**2
+#            z = -sqrt(du)*.25+(0.8975)
+#            if z <0:
+#                z=0
+#            if z > .625:
+#                z=0.625
+#
+#            #print 'x,y,f',du,i,j,z
+#            file.write("%s, %s, %s\n" %(i/res, j/res, z))
+#    file.close()
+#
+##square_root()
+##test()
+#
+#square_root_data=array([],typecode=Float)
+#
+#square_root_data=resize(square_root_data,(400,2))
+#
+#
+#for i in range (0,37):
+#    for j in range (0,10):
+#        number = i+(j/10.)
+#        print "%s,%s, %s, %s" %(i,j, number, sqrt(number))
+#        
+#        square_root_data[i*10+j]= number, sqrt(number)
+#            
+##        file1.write("%s, %s\n" %(i, sqrt(number)))
+
+attribute_dic, title_index_dic = csv2dict('sqrt_table.csv')
+#print attribute_dic
+print attribute_dic['number'][10]
+
+#print title_index_dic
+
+
+def test(x,y):
+#    import anug square_root                
+    center_x = 12.96
+#    center_y = 13.80
+    center_y = 6.80
+
+    du = (center_x-x)**2+(center_y-y)**2
+#    n=int(round(du,0))
+    print 'x',min(x),max(x)
+    print 'y',min(y),max(y)
+    z1=[]
+    for d in du:
+        n1a = round(d,0)
+        n1=int(n1a)*10
+        n2=int(round(d - int(round(d,0)),1)*10)
+        nn=n1+n2
+        #print d,n1a,n1,n2,nn
+#    zz=square_root_data[nn]
+        zz=-float(attribute_dic['sqrt'][nn])
+        print nn,zz
+
+#    z = -square_root((center_x-x)**2+(center_y-y)**2)
+        z = zz*.25+0.8975
+
+        if z <0:
+            z=0
+        if z > .625:
+            z=0.625
+        print 'hello',z
+        z1.append(z)
+    return z1
+
+#domain.set_quantity('elevation',filename = "test.csv", alpha=0.1)
+#domain.set_quantity('elevation',filename = "test.csv", alpha=1.,verbose=True)
+domain.set_quantity('elevation',test, verbose=True)
 
 
@@ -246 +324 @@
 #    domain.write_time(track_speeds=False)
 
-#for t in domain.evolve(yieldstep = 0.1, finaltime = 36
-for t in domain.evolve(yieldstep = 0.5, finaltime = 36
+for t in domain.evolve(yieldstep = 0.1, finaltime = 36
+#for t in domain.evolve(yieldstep = 0.5, finaltime = 36
                        ,skip_initial_step = True): 
     domain.write_time()