Changeset 4457


Ignore:
Timestamp:
May 16, 2007, 4:36:20 PM (17 years ago)
Author:
nick
Message:

working version of alpha validation

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/alpha_validation/find_alpha.py

    r4427 r4457  
    66from anuga.pmesh.mesh_interface import create_mesh_from_regions
    77from anuga.shallow_water import Reflective_boundary
    8 
     8from anuga.utilities.numerical_tools import cov
     9from Numeric import array, resize,shape,Float,zeros
    910
    1011
     
    1819                         maximum_triangle_area=100000,
    1920                         filename=mesh_dir_name,
    20                          use_cache=False,
     21                         use_cache=True,
    2122                         verbose=True)
    2223
     
    3031G_small, G_other = G.split(0.1,True)
    3132
    32 #print 'start export'
    33 #G_small.export_points_file(sample)
    34 #G_other.export_points_file(remainder)
    3533
    3634
    37 #alphas = [0.001, 0.01, 0.1, 1.0, 10.0, 100.0]
    38 alphas = [1]
     35#alphas = [0.0001, 0.001, 0.01, 0.1, 1.0, 10.0, 100.0,1000.0,10000.0]
     36alphas = [0.001,1,100]
    3937domains = {}
     38#domains = []
    4039
    4140print 'Setup computational domain'
     
    4342   
    4443    #think about making domain absolute when it is created
    45     geo_ref=Geo_reference(zone = 0, xllcorner = 0.0, yllcorner = 0.0)
    46    
    47 #    domain = Domain(mesh_dir_name,geo_reference=geo_ref, use_cache=False, verbose=True)
    4844    domain = Domain(mesh_dir_name, use_cache=False, verbose=True)
    4945    print domain.statistics()
     
    5450                    alpha = alpha)
    5551    domains[alpha]=domain
    56     domain.set_name('temp'+str(alpha))
    57     domain.set_datadir('.')
    58     print'set_boundary'
    59 #    points=[[665000,7753000],[665100,7753000]]
    60     points=G_small.get_data_points()
    61     print'geo',domain.geo_reference,points[0:10]
    62     #gets relative point from sample
    63     points=domain.geo_reference.change_points_geo_ref(points)
     52
     53points=G_small.get_data_points()
     54data=array([],typecode=Float)
     55data=resize(data,(len(points),3+len(alphas)))
     56#gets relative point from sample
     57
     58data[:,0]=points[:,0]
     59data[:,1]=points[:,1]
     60print'geo',domain.geo_reference,points[0:10],points[0:10,0]
     61elevation_sample=G_small.get_attributes(attribute_name='elevation')
     62
     63data[:,2]=elevation_sample
     64
     65print'data',data[0:10],len(alphas)
     66normal_cov=array(zeros([len(alphas),2]),typecode=Float)
     67i=0
     68#print 'domains',domains.keys()
     69#keys_list=domains.keys()
     70#print 'domains',keys_list.sort()
     71for domain in domains:
     72#    print'domain',domain
     73    points_geo=domains[domain].geo_reference.change_points_geo_ref(points)
    6474   
    65     point_elevation=G_small.get_attributes(attribute_name='elevation')
    6675
    67     print'poin',points[0:10],point_elevation[0:10]
    68     print'quantities',domain.get_quantity_names()
    69     hmm=domain.quantities['elevation'].get_values(interpolation_points=points)
    70     print'hmm',hmm[0:10],len(points),len(hmm)
     76    print'poin',points_geo[0:10],elevation_sample[0:10]
     77    print'quantities',domains[domain].get_quantity_names()
     78    elevation_predicted=domains[domain].quantities['elevation'].get_values(interpolation_points=points_geo)
     79    print'elevation_predicted',elevation_predicted[0:10],len(points),len(elevation_predicted)
    7180   
    72 #    for i in si
     81    data[:,i+3]=elevation_predicted
     82    sample_cov= cov(elevation_sample)
     83
     84    ele_cov= cov(elevation_sample,elevation_predicted)
     85    normal_cov[i,:]= [domain,ele_cov/sample_cov]
     86
     87    print'cov',alpha,'= ',normal_cov, ele_cov, sample_cov
     88    i+=1
    7389   
    74 #    Br = Reflective_boundar[y(domain)
    7590   
    76 #    domain.set_boundary({'back': Br,
    77 #                     'side': Br,
    78 #                     'ocean': Br})
    79 
    80 #    for t in domain.evolve(yieldstep = 1, finaltime = 1):
    81 #        domain.write_time()
    82        
    83 
    84 
     91print'data',data[100:200]
     92print'normal cov',normal_cov
     93   
    8594
    8695# from file load x,y,z data
Note: See TracChangeset for help on using the changeset viewer.