[2424] | 1 | """Stochastic study of the ANUGA implementation of the |
---|
| 2 | shallow water wave equation. |
---|
| 3 | |
---|
| 4 | This script runs the model for one realisation of bathymetry as |
---|
| 5 | given in the file bathymetry.txt and outputs a full simulation is \ |
---|
| 6 | sww NetCDF format. |
---|
| 7 | |
---|
| 8 | The left boundary condition is a timeseries defined in |
---|
| 9 | NetCDF file: input_wave.tms |
---|
| 10 | |
---|
| 11 | Note: This scripts needs create_mesh.py to have been run |
---|
| 12 | |
---|
| 13 | Suresh Kumar and Ole Nielsen 2006 |
---|
| 14 | """ |
---|
| 15 | |
---|
| 16 | |
---|
| 17 | #------------------------------------------------------------------------------ |
---|
| 18 | # Import necessary modules |
---|
| 19 | #------------------------------------------------------------------------------ |
---|
| 20 | |
---|
| 21 | # Standard modules |
---|
| 22 | import os |
---|
| 23 | import time |
---|
[2448] | 24 | import cPickle |
---|
[2424] | 25 | |
---|
| 26 | # Related major packages |
---|
| 27 | from pyvolution.shallow_water import Domain |
---|
| 28 | from pyvolution.shallow_water import Reflective_boundary |
---|
| 29 | from pyvolution.shallow_water import Transmissive_Momentum_Set_Stage_boundary |
---|
| 30 | from pyvolution.pmesh2domain import pmesh_to_domain_instance |
---|
| 31 | from pyvolution.data_manager import xya2pts |
---|
| 32 | from pyvolution.util import file_function |
---|
[2484] | 33 | from caching.caching import cache |
---|
[2424] | 34 | |
---|
| 35 | # Application specific imports |
---|
| 36 | import project # Definition of file names and polygons |
---|
| 37 | |
---|
| 38 | |
---|
[2984] | 39 | #----------------------------------------------------------------------------- |
---|
| 40 | # Read in processor information |
---|
| 41 | #----------------------------------------------------------------------------- |
---|
[2448] | 42 | |
---|
[3008] | 43 | try: |
---|
| 44 | import pypar |
---|
| 45 | except: |
---|
| 46 | print 'Could not import pypar' |
---|
| 47 | myid = 0 |
---|
| 48 | numprocs = 1 |
---|
[3009] | 49 | processor_name = 'local host' |
---|
[3008] | 50 | else: |
---|
| 51 | myid = pypar.rank() |
---|
| 52 | numprocs = pypar.size() |
---|
| 53 | processor_name = pypar.Get_processor_name() |
---|
[2448] | 54 | |
---|
[2984] | 55 | print 'I am process %d of %d running on %s' %(myid, numprocs, processor_name) |
---|
| 56 | |
---|
| 57 | |
---|
[2424] | 58 | #----------------------------------------------------------------------------- |
---|
| 59 | # Setup computational domain |
---|
| 60 | #----------------------------------------------------------------------------- |
---|
[2984] | 61 | #print 'Creating domain from', project.mesh_filename |
---|
[2424] | 62 | |
---|
[2979] | 63 | domain = Domain(project.working_dir + project.mesh_filename, |
---|
[2966] | 64 | use_cache=False, |
---|
[2984] | 65 | verbose=False) |
---|
[2866] | 66 | |
---|
| 67 | |
---|
[2984] | 68 | #print 'Number of triangles = ', len(domain) |
---|
| 69 | #print domain.statistics() |
---|
[2424] | 70 | |
---|
[2535] | 71 | |
---|
[2979] | 72 | domain.set_datadir(project.working_dir) |
---|
[2424] | 73 | domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum']) |
---|
| 74 | |
---|
| 75 | |
---|
| 76 | #------------------------------------------------------------------------------ |
---|
| 77 | # Setup boundary conditions |
---|
| 78 | #------------------------------------------------------------------------------ |
---|
| 79 | |
---|
[2984] | 80 | function = file_function(project.boundary_filename, domain, verbose = False) |
---|
[2424] | 81 | Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function) #Input wave |
---|
| 82 | Br = Reflective_boundary(domain) #Wall |
---|
| 83 | |
---|
| 84 | # Bind boundary objects to tags |
---|
| 85 | domain.set_boundary({'wave': Bts, 'wall': Br}) |
---|
| 86 | |
---|
| 87 | |
---|
| 88 | #------------------------------------------------------------------------------ |
---|
[2448] | 89 | # Setup initial conditions |
---|
[2424] | 90 | #------------------------------------------------------------------------------ |
---|
[2448] | 91 | domain.set_quantity('friction', 0.0) |
---|
| 92 | domain.set_quantity('stage', 0.0) |
---|
[2424] | 93 | |
---|
[2448] | 94 | # Get prefitted realisations |
---|
[2484] | 95 | |
---|
| 96 | finaltime = 22.5 |
---|
| 97 | timestep = 0.05 |
---|
| 98 | |
---|
| 99 | |
---|
[3009] | 100 | realisation = 0 |
---|
[2979] | 101 | for filename in os.listdir(project.working_dir): |
---|
[2448] | 102 | if filename.startswith(project.basename) and filename.endswith('.pck'): |
---|
[3008] | 103 | print 'P%d: Reading %s' %(myid, filename) |
---|
| 104 | fid = open(project.working_dir + filename) |
---|
[2448] | 105 | V = cPickle.load(fid) |
---|
| 106 | fid.close() |
---|
| 107 | |
---|
[2989] | 108 | #if myid == 0: |
---|
| 109 | # print 'V', V[6:7,:] |
---|
| 110 | |
---|
[2448] | 111 | # For each column (each realisation) |
---|
| 112 | for i in range(V.shape[1]): |
---|
| 113 | |
---|
[2984] | 114 | # Distribute work in round-robin fashion |
---|
| 115 | if i%numprocs == myid: |
---|
| 116 | |
---|
| 117 | name = project.basename + '_P%d' %myid |
---|
| 118 | domain.set_name(name) #Output name |
---|
[3008] | 119 | print 'V', V.shape |
---|
[2984] | 120 | domain.set_quantity('elevation', V[:,i]) #Assign bathymetry |
---|
[2989] | 121 | |
---|
| 122 | print 'P%d: Setting quantity %d: %s' %(myid, i, str(V[:4,i])) |
---|
| 123 | |
---|
[2984] | 124 | domain.set_time(0.0) #Reset time |
---|
[2966] | 125 | |
---|
[2984] | 126 | #--------------------------------------------------- |
---|
| 127 | # Evolve system through time |
---|
| 128 | #--------------------------------------------------- |
---|
| 129 | print 'P%d: Running realisation %d of %d in block %s'\ |
---|
| 130 | %(myid, realisation, V.shape[1], filename) |
---|
| 131 | |
---|
| 132 | t0 = time.time() |
---|
| 133 | for t in domain.evolve(yieldstep = timestep, |
---|
| 134 | finaltime = finaltime): |
---|
| 135 | pass |
---|
[3035] | 136 | domain.write_time() |
---|
[2984] | 137 | |
---|
[2448] | 138 | |
---|
[2984] | 139 | print 'P%d: Simulation of realisation %d took %.2f seconds'\ |
---|
| 140 | %(myid, realisation, time.time()-t0) |
---|
[2448] | 141 | |
---|
[2484] | 142 | |
---|
| 143 | |
---|
| 144 | |
---|
[2984] | 145 | #--------------------------------------------------- |
---|
| 146 | # Now extract the 3 timeseries (Ch 5-7-9) and store them |
---|
| 147 | # in three files for this realisation |
---|
[2989] | 148 | |
---|
| 149 | print 'P%d: Extracting time series for realisation %d from file %s'\ |
---|
| 150 | %(myid, realisation, project.working_dir + domain.filename + '.sww') |
---|
[2984] | 151 | f = file_function(project.working_dir + domain.filename + '.sww', |
---|
| 152 | quantities='stage', |
---|
| 153 | interpolation_points=project.gauges, |
---|
| 154 | verbose=False) |
---|
[2484] | 155 | |
---|
| 156 | |
---|
[2984] | 157 | simulation_name = project.working_dir + \ |
---|
| 158 | project.basename + '_realisation_%d' %realisation |
---|
[2484] | 159 | |
---|
[2989] | 160 | print 'P%d: Writing to file %s'\ |
---|
| 161 | %(myid, simulation_name + '_' + name + '.txt') |
---|
| 162 | |
---|
[2984] | 163 | for k, name in enumerate(project.gauge_names): |
---|
| 164 | fid = open(simulation_name + '_' + name + '.txt', 'w') |
---|
| 165 | for t in f.get_time(): |
---|
| 166 | #For all precomputed timesteps |
---|
| 167 | val = f(t, point_id = k)[0] |
---|
| 168 | fid.write('%f %f\n' %(t, val)) |
---|
[2484] | 169 | |
---|
[2984] | 170 | fid.close() |
---|
[2484] | 171 | |
---|
| 172 | |
---|
| 173 | |
---|
| 174 | realisation += 1 |
---|
[2984] | 175 | |
---|
| 176 | |
---|
| 177 | pypar.finalize() |
---|