source: development/stochastic_study/run_model.py @ 2878

"""Stochastic study of the ANUGA implementation of the
shallow water wave equation.

This script runs the model for one realisation of bathymetry as
given in the file bathymetry.txt and outputs a full simulation is \
sww NetCDF format.

The left boundary condition is a timeseries defined in
NetCDF file: input_wave.tms

Note: This scripts needs create_mesh.py to have been run

Suresh Kumar and Ole Nielsen 2006
"""


#------------------------------------------------------------------------------
# Import necessary modules
#------------------------------------------------------------------------------

# Standard modules
import os
import time
import cPickle

# Related major packages
from pyvolution.shallow_water import Domain
from pyvolution.shallow_water import Reflective_boundary
from pyvolution.shallow_water import Transmissive_Momentum_Set_Stage_boundary
from pyvolution.pmesh2domain import pmesh_to_domain_instance
from pyvolution.data_manager import xya2pts
from pyvolution.util import file_function
from caching.caching import cache

# Application specific imports
import project                 # Definition of file names and polygons




#-----------------------------------------------------------------------------
# Setup computational domain
#----------------------------------------------------------------------------- 
print 'Creating domain from', project.mesh_filename

#domain = pmesh_to_domain_instance(project.mesh_filename, Domain,
#                                  use_cache=True,
#                                  verbose=True)

domain = Domain(project.mesh_filename,
                use_cache=True,
                verbose=True)                
                

print 'Number of triangles = ', len(domain)
print 'The extent is ', domain.get_extent()
print domain.statistics()

import sys; sys.exit()


domain.set_datadir('.')
domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])

#domain.check_integrity()


#------------------------------------------------------------------------------
# Setup boundary conditions
#------------------------------------------------------------------------------

function = file_function(project.boundary_filename, domain, verbose = True)
Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function) #Input wave
Br = Reflective_boundary(domain) #Wall

# Bind boundary objects to tags
domain.set_boundary({'wave': Bts, 'wall': Br})


#------------------------------------------------------------------------------
# Setup initial conditions
#------------------------------------------------------------------------------
domain.set_quantity('friction', 0.0)
domain.set_quantity('stage', 0.0)

# Get prefitted realisations

finaltime = 22.5
timestep = 0.05



realisation = 0
for filename in os.listdir('.'):
    if filename.startswith(project.basename) and filename.endswith('.pck'):
        print 'Reading %s' %filename
        fid = open(filename)
        V = cPickle.load(fid)
        fid.close()

        # For each column (each realisation)
        for i in range(V.shape[1]):
            domain.set_name(project.basename)         #Output name
            domain.set_quantity('elevation', V[:,i])  #Assign bathymetry
            domain.starttime = 0.0                    #Reset time

            #---------------------------------------------------
            # Evolve system through time
            #---------------------------------------------------
            print 'Running realisation %d of %d in block %s'\
                  %(i, V.shape[1], filename)
            t0 = time.time()
            for t in domain.evolve(yieldstep = timestep, finaltime = finaltime):
                domain.write_time()
                
        
            print 'Realisation %d took %.2f seconds'\
                  %(realisation, time.time()-t0)




            #---------------------------------------------------
            # Now extract the 3 timeseries (Ch 5-7-9) and store them
            # in three files for this realisation
            gauges = [[4.521, 1.196],  [4.521, 1.696],  [4.521, 2.196]] 
            gauge_names = ['ch5', 'ch7', 'ch9']

            
            f = file_function(domain.filename + '.sww',
                              quantities='stage',
                              interpolation_points=gauges,
                              verbose = True)


            simulation_name = domain.filename + '_realisation_%d' %realisation

            for k, name in enumerate(gauge_names):
                fid = open(simulation_name + '_' + name + '.txt', 'w')
                for t in f.T:
                    #For all precomputed timesteps
                    val = f(t, point_id = k)[0]
                    fid.write('%f %f\n' %(t, val))

                fid.close()

                    

            realisation += 1