Changeset 4980


Ignore:
Timestamp:
Jan 30, 2008, 9:41:35 AM (16 years ago)
Author:
nick
Message:

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_validation/conical_island/run_circular.py

    r4979 r4980  
    1111from anuga.shallow_water import Dirichlet_boundary, Time_boundary
    1212from anuga.shallow_water.data_manager import get_maximum_inundation_data
    13 from anuga.abstract_2d_finite_volumes.util import file_function, sww2csv_gauges
     13from anuga.abstract_2d_finite_volumes.util import file_function, sww2csv_gauges,square_root
    1414from anuga.pmesh.mesh_interface import create_mesh_from_regions
    1515from anuga.utilities.polygon import read_polygon, plot_polygons
    1616#from cmath import cos,sin,cosh,pi
    1717from math import cos,pi,sin,tan,sqrt
    18 from Numeric import array, zeros, Float, allclose
     18from Numeric import array, zeros, Float, allclose,resize
     19from anuga.shallow_water.data_manager import csv2dict
    1920
    2021#import project
     
    9899d=0.75
    99100res=.05
     101res=1.
    100102
    101103#create polygons (circles) to have higher resolution
     
    164166domain.set_quantity('stage', water_height)
    165167
    166 def test():               
    167     fileName = "test.csv"
    168     file = open(fileName,"w")
    169     file.write("easting,northing,elevation \n")
    170     x = 12.96
    171     y = 13.80
    172     for i in range(0,30):
    173         for j in range(0,25):
    174             a = (x-i)**2
    175             b = (y-j)**2
    176             z = -sqrt((x-i)**2+(y-j)**2)
    177             z = z*.25+(0.8975)
    178             if z <0:
    179                 z=0
    180             if z > .625:
    181                 z=0.625
    182 
    183 #            print 'x,y,f',i,j,a,b,z
    184             file.write("%s, %s, %s\n" %(i, j, z))
    185     file.close()
    186 
    187 test()
    188 
    189 domain.set_quantity('elevation',filename = "test.csv", alpha=0.1)
     168#def square_root():
     169#    fileName1 = "square_root_table.csv"
     170#    file1 = open(fileName1,"w")
     171##    file1.write("number,square_root \n")
     172#    for i in range (0,37):
     173#        for j in range (0,10):
     174#            number = i+(j/10.)
     175#            print "%s,%s, %s, %s" %(i,j, number, sqrt(number))
     176#           
     177#            file1.write("%s, %s\n" %(i, sqrt(number)))
     178#    file1.close()
     179#
     180##    return sqrt(s)
     181#
     182#def test():               
     183#    fileName = "test.csv"
     184#    file = open(fileName,"w")
     185#    file.write("easting,northing,elevation \n")
     186#   
     187#   # csv2dict()
     188#   
     189#    x = 12.96
     190#    y = 13.80
     191#    res=1
     192#    for i in range(0,30*res):
     193#        for j in range(0,25*res):
     194##            a = (x-i)**2
     195##            b = (y-j)**2
     196#
     197#          #  z = -sqrt((x-(i/res))**2+(y-(j/res))**2)
     198#            du = (x-(i/res))**2+(y-(j/res))**2
     199#            z = -sqrt(du)*.25+(0.8975)
     200#            if z <0:
     201#                z=0
     202#            if z > .625:
     203#                z=0.625
     204#
     205#            #print 'x,y,f',du,i,j,z
     206#            file.write("%s, %s, %s\n" %(i/res, j/res, z))
     207#    file.close()
     208#
     209##square_root()
     210##test()
     211#
     212#square_root_data=array([],typecode=Float)
     213#
     214#square_root_data=resize(square_root_data,(400,2))
     215#
     216#
     217#for i in range (0,37):
     218#    for j in range (0,10):
     219#        number = i+(j/10.)
     220#        print "%s,%s, %s, %s" %(i,j, number, sqrt(number))
     221#       
     222#        square_root_data[i*10+j]= number, sqrt(number)
     223#           
     224##        file1.write("%s, %s\n" %(i, sqrt(number)))
     225
     226attribute_dic, title_index_dic = csv2dict('sqrt_table.csv')
     227#print attribute_dic
     228print attribute_dic['number'][10]
     229
     230#print title_index_dic
     231
     232
     233def test(x,y):
     234#    import anug square_root               
     235    center_x = 12.96
     236#    center_y = 13.80
     237    center_y = 6.80
     238
     239    du = (center_x-x)**2+(center_y-y)**2
     240#    n=int(round(du,0))
     241    print 'x',min(x),max(x)
     242    print 'y',min(y),max(y)
     243    z1=[]
     244    for d in du:
     245        n1a = round(d,0)
     246        n1=int(n1a)*10
     247        n2=int(round(d - int(round(d,0)),1)*10)
     248        nn=n1+n2
     249        #print d,n1a,n1,n2,nn
     250#    zz=square_root_data[nn]
     251        zz=-float(attribute_dic['sqrt'][nn])
     252        print nn,zz
     253
     254#    z = -square_root((center_x-x)**2+(center_y-y)**2)
     255        z = zz*.25+0.8975
     256
     257        if z <0:
     258            z=0
     259        if z > .625:
     260            z=0.625
     261        print 'hello',z
     262        z1.append(z)
     263    return z1
     264
     265#domain.set_quantity('elevation',filename = "test.csv", alpha=0.1)
     266#domain.set_quantity('elevation',filename = "test.csv", alpha=1.,verbose=True)
     267domain.set_quantity('elevation',test, verbose=True)
    190268
    191269
     
    246324#    domain.write_time(track_speeds=False)
    247325
    248 #for t in domain.evolve(yieldstep = 0.1, finaltime = 36
    249 for t in domain.evolve(yieldstep = 0.5, finaltime = 36
     326for t in domain.evolve(yieldstep = 0.1, finaltime = 36
     327#for t in domain.evolve(yieldstep = 0.5, finaltime = 36
    250328                       ,skip_initial_step = True):
    251329    domain.write_time()
Note: See TracChangeset for help on using the changeset viewer.