source: trunk/anuga_work/development/stochastic_study/create_realisations.py @ 7884

Last change on this file since 7884 was 3514, checked in by duncan, 18 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

File size: 2.9 KB
Line 
1"""Create several realisations of bathymetries for use with run_model.
2
3This script will create a number of random samples of the input bathymetry and
4then fit all of them to arrays that can be used in set_quantity inside run_model.
5
6The reason for this is that doing all fits in one sweep saves computation time
7"""
8
9from RandomArray import normal, seed
10from Numeric import concatenate, NewAxis
11import cPickle
12from load_mesh.loadASCII import import_points_file, import_mesh_file
13from anuga.pyvolution.least_squares import fit_to_mesh
14import project
15
16bathymetry_filename = project.bathymetry_filename
17mesh_filename = project.working_dir + project.mesh_filename
18
19
20# Read mesh
21mesh_dict = import_mesh_file(mesh_filename)
22vertex_coordinates = mesh_dict['vertices']
23triangles = mesh_dict['triangles']
24
25# Read basic bathymetry from NetCDF .pts file
26print 'Reading bathymetry'
27base_bathymetry = import_points_file(bathymetry_filename)
28
29
30# Extract relevant Numeric arrays
31print 'Extract bathymetry'
32point_coordinates = base_bathymetry['pointlist']
33elevation = base_bathymetry['attributelist']['elevation']
34N = len(elevation)
35print 'read', N, 'points:' + 'z in [%f, %f]' %(min(elevation),
36                                               max(elevation))
37
38# Create different bathymetries based on sampling (Suresh to fill this in)
39print 'Create different realisations'
40Z = None
41block_number = 0
42seed(13,17)
43for n in range(project.number_of_realisations):
44
45    print 'Creating sample #%d' %n
46    # nth realisation of bathymetry
47
48    if n == 0:
49        z = elevation # Keep the first realisation as original
50    else:   
51        z = elevation + normal(project.mean,
52                               project.std_dev,
53                               N)
54
55    #Concatenate into array with each realisation a column
56    if Z is None:
57        Z = z[:, NewAxis] # The first realisation of this block
58    else:   
59        Z = concatenate((Z, z[:, NewAxis]), axis=1)
60
61    if (n+1)%project.blocksize == 0 or n == project.number_of_realisations-1:
62        #Fit all of them to computational mesh
63        print 'fit block %d to the mesh' %block_number
64        V = fit_to_mesh(vertex_coordinates,
65                        triangles,
66                        point_coordinates,
67                        Z,
68                        alpha = 0.0001,
69                        verbose = True,
70                        expand_search = True,
71                        precrop = True,
72                        use_cache = False)
73
74        Z = None
75
76        #Store this block for use with model
77        print 'Store fitted quantities in block %d for use with model'\
78              %block_number
79
80        filename = project.working_dir + project.basename +\
81                   '_block%d.pck' %block_number
82
83        block_number += 1   
84       
85       
86        fid = open(filename, 'wb')
87        cPickle.dump(V, fid)
88        fid.close()
89       
Note: See TracBrowser for help on using the repository browser.