source: development/stochastic_study/run_simple_model.py @ 3514

Last change on this file since 3514 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: 5.6 KB
Line 
1"""Stochastic study of the ANUGA implementation of the
2shallow water wave equation.
3
4Very simple flat bed with one spike in the middle where different values are sampled.
5Three sides are reflective and one Dirichlet condition
6
7Output is measured at gauges around the spike.
8
9The system will evolve dynamically towards a steady state.
10
11Suresh Kumar and Ole Nielsen 2006
12"""
13
14
15#------------------------------------------------------------------------------
16# Import necessary modules
17#------------------------------------------------------------------------------
18
19# Standard modules
20import os
21import time
22
23# Related major packages
24from Numeric import allclose, zeros, Float
25from RandomArray import normal, seed
26
27# Application specific imports
28from anuga.pyvolution.mesh_factory import rectangular
29from anuga.pyvolution.shallow_water import Domain
30from anuga.pyvolution.shallow_water import Reflective_boundary
31from anuga.pyvolution.shallow_water import Dirichlet_boundary
32from anuga.pyvolution.util import file_function
33from caching.caching import cache
34import project                 # Definition of file names and polygons
35
36
37#-----------------------------------------------------------------------------
38# Read in processor information
39#-----------------------------------------------------------------------------
40
41try:
42    import pypar
43except:
44    print 'Could not import pypar'
45    myid = 0
46    numprocs = 1
47    processor_name = 'local host'
48else:   
49    myid = pypar.rank()
50    numprocs = pypar.size()
51    processor_name = pypar.Get_processor_name()
52
53print 'I am process %d of %d running on %s' %(myid, numprocs, processor_name)
54
55
56#-----------------------------------------------------------------------------
57# Setup computational domain
58#-----------------------------------------------------------------------------
59
60# Create basic triangular mesh
61points, vertices, boundary = rectangular(10, 10, 10, 10)
62
63#Create shallow water domain
64domain = Domain(points, vertices, boundary)
65domain.set_name(project.basename)
66domain.set_datadir(project.working_dir) 
67print domain.statistics()
68
69#------------------------------------------------------------------------------
70# Setup boundary conditions
71#------------------------------------------------------------------------------
72
73Br = Reflective_boundary(domain) #Wall
74Bd = Dirichlet_boundary([1, 0, 0])
75
76# Bind boundary objects to tags
77domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
78
79
80#------------------------------------------------------------------------------
81# Setup initial conditions
82#------------------------------------------------------------------------------
83domain.set_quantity('friction', 0.0)
84domain.set_quantity('stage', 0.0)
85
86# Initial conditions
87def f(x,y):
88    z = zeros( x.shape, Float )
89    for i in range(len(x)):
90        if allclose([x[i], y[i]], 5):
91            print i
92            z[i] = 1
93    return z
94
95domain.set_quantity('elevation', f)
96
97
98z = domain.get_quantity('elevation').get_values()
99print z
100#print z != 0
101
102import sys; sys.exit() 
103
104# Run model for different realisations of spike
105print 'Create %d different realisations' %project.number_of_realisations
106
107# Array of bed elevation values at vertices
108
109#print z
110N = z.shape[0]
111
112timestep = 0.1
113finaltime = 1
114
115seed(13,17)
116for realisation in range(project.number_of_realisations):
117
118    print 'Creating realisation #%d' %realisation
119    #z[N/2-1,:] = normal(project.mean, project.std_dev, 1)
120
121    # Distribute work in round-robin fashion
122    if realisation%numprocs == myid:
123        name = project.basename + '_P%d' %myid   
124        domain.set_name(name)                # Output name
125        #domain.set_quantity('elevation', z)  # Assign bathymetry
126        domain.set_time(0.0)                 # Reset time
127
128        #---------------------------------------------------
129        # Evolve system through time
130        #---------------------------------------------------
131        print 'P%d: Running realisation %d of %d'\
132              %(myid, realisation, project.number_of_realisations)
133
134        t0 = time.time()
135        for t in domain.evolve(yieldstep = timestep,
136                               finaltime = finaltime):
137            pass
138            domain.write_time()
139                   
140       
141        print 'P%d: Simulation of realisation %d took %.2f seconds'\
142              %(myid, realisation, time.time()-t0)
143
144
145
146       
147        #---------------------------------------------------
148        # Now extract the 3 timeseries (Ch 5-7-9) and store them
149        # in three files for this realisation
150
151        print 'P%d: Extracting time series for realisation %d from file %s'\
152              %(myid, realisation, project.working_dir + domain.filename + '.sww')
153        f = file_function(project.working_dir + domain.filename + '.sww',
154                          quantities='stage',
155                          interpolation_points=project.gauges,
156                          verbose=False)
157
158
159        simulation_name = project.working_dir + \
160                          project.basename + '_realisation_%d' %realisation
161
162        print 'P%d: Writing to file %s'\
163              %(myid, simulation_name + '_' + name + '.txt')
164       
165        for k, name in enumerate(project.gauge_names):
166            fid = open(simulation_name + '_' + name + '.txt', 'w')
167            for t in f.get_time():
168                #For all precomputed timesteps
169                val = f(t, point_id = k)[0]
170                fid.write('%f %f\n' %(t, val))
171               
172            fid.close()
173
174           
175pypar.finalize()
Note: See TracBrowser for help on using the repository browser.