Changeset 4457
- Timestamp:
- May 16, 2007, 4:36:20 PM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/development/alpha_validation/find_alpha.py
r4427 r4457 6 6 from anuga.pmesh.mesh_interface import create_mesh_from_regions 7 7 from anuga.shallow_water import Reflective_boundary 8 8 from anuga.utilities.numerical_tools import cov 9 from Numeric import array, resize,shape,Float,zeros 9 10 10 11 … … 18 19 maximum_triangle_area=100000, 19 20 filename=mesh_dir_name, 20 use_cache= False,21 use_cache=True, 21 22 verbose=True) 22 23 … … 30 31 G_small, G_other = G.split(0.1,True) 31 32 32 #print 'start export'33 #G_small.export_points_file(sample)34 #G_other.export_points_file(remainder)35 33 36 34 37 #alphas = [0.00 1, 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] 36 alphas = [0.001,1,100] 39 37 domains = {} 38 #domains = [] 40 39 41 40 print 'Setup computational domain' … … 43 42 44 43 #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)48 44 domain = Domain(mesh_dir_name, use_cache=False, verbose=True) 49 45 print domain.statistics() … … 54 50 alpha = alpha) 55 51 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 53 points=G_small.get_data_points() 54 data=array([],typecode=Float) 55 data=resize(data,(len(points),3+len(alphas))) 56 #gets relative point from sample 57 58 data[:,0]=points[:,0] 59 data[:,1]=points[:,1] 60 print'geo',domain.geo_reference,points[0:10],points[0:10,0] 61 elevation_sample=G_small.get_attributes(attribute_name='elevation') 62 63 data[:,2]=elevation_sample 64 65 print'data',data[0:10],len(alphas) 66 normal_cov=array(zeros([len(alphas),2]),typecode=Float) 67 i=0 68 #print 'domains',domains.keys() 69 #keys_list=domains.keys() 70 #print 'domains',keys_list.sort() 71 for domain in domains: 72 # print'domain',domain 73 points_geo=domains[domain].geo_reference.change_points_geo_ref(points) 64 74 65 point_elevation=G_small.get_attributes(attribute_name='elevation')66 75 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) 71 80 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 73 89 74 # Br = Reflective_boundar[y(domain)75 90 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 91 print'data',data[100:200] 92 print'normal cov',normal_cov 93 85 94 86 95 # from file load x,y,z data
Note: See TracChangeset
for help on using the changeset viewer.