source: anuga_work/development/ted_rigby/model_run.py @ 7839

"""
----------------------------------------------------------------- model_run.py
ANUGA RURAL FLOOD MODELLING TEMPLATE (Domain Construction & Execution)
Script for building the model domain and setting initial conditions and 
(temporally varying)  boundary conditions appropriate to the particular 
scenario-event before running (evolving) the simulation.

Most model data for the run is read in or otherwise referenced in 
model_data.py and imported into model_run.py to create the model mesh and 
domain before running the simulation. 

Code Version : 1.00 June 2008 Initial release
               1.01 December 2008 Vble_n recoded to speed execution
Author       : E Rigby (0437 250 500) ted.rigby@rienco.com.au

-------------------------------------------------------------------------------
"""

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

# Standard python modules
import os
import time
import sys
from Numeric import zeros, Float

# Related ANUGA modules
from anuga.shallow_water import Domain
from anuga.shallow_water import Reflective_boundary
from anuga.shallow_water import Dirichlet_boundary
from anuga.shallow_water import Time_boundary
from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
from anuga.shallow_water.shallow_water_domain import Inflow
from anuga.abstract_2d_finite_volumes.util import file_function
from anuga.abstract_2d_finite_volumes.quantity import Quantity
from anuga.pmesh.mesh_interface import create_mesh_from_regions
from anuga.utilities.polygon import read_polygon, plot_polygons, Polygon_function, inside_polygon
from anuga.alpha_shape.alpha_shape import alpha_shape_via_files


# Model specific imports
import model_data                # Defines and provides the scenario-event specific data 
import get_tuflow_data           # Provides the routines to read in tuflow specific data

#####################################################################################
#                                                                                   #
#                     S I M U L A T I O N    M A I N L I N E                        #
#                                                                                   #
#####################################################################################

t0 = time.time()                 # record start time of this run (secs)
print model_data.basename+' >>>> Comencing model_run.py code execution at t= %.2f hours'  %(float(time.time()-t0)/3600.0)

#####################################################################################
#                                                                                   #
#                          CREATE THE TRIANGULAR MESH                               #
#                                                                                   #
#####################################################################################

# Note: The mesh is created based on overall clipping polygon with a tagged  
# boundary and interior regions as defined in model_data.py along with       
# resolutions (maximal area per triangle) for each interior polygon       

print model_data.basename+' >>>> Creating variable mesh within bounding_poly reflecting regional resolutions'
# Create the mesh reflecting the bounding (default) and internal region mesh resolutions set in model_data.py
# with boundary tags as defined in model_data.py

create_mesh_from_regions(model_data.bounding_polygon,
                         boundary_tags = model_data.boundary_strings,
                         maximum_triangle_area = model_data.default_res,
                         minimum_triangle_angle=30.0,
                         interior_regions = model_data.interior_resregions,
                         filename = model_data.mesh_filename,
                         use_cache=True,
                         verbose=model_data.anuga_verbose)


print model_data.basename+' >>>> Completed mesh creation at t= %.2f hours ' %(float(time.time()-t0)/3600.0)

######################################################################################
#                                                                                    #
#                       CREATE THE COMPUTATIONAL DOMAIN                              #
#                                                                                    #
######################################################################################

print model_data.basename+' >>>> Creating domain from mesh and setting domain parameters'
# create domain from mesh meshname
domain = Domain(model_data.mesh_filename, use_cache=True, verbose=model_data.anuga_verbose) 
print domain.statistics()                                # confirm what has been done

# set base domain parameters
domain.set_name(model_data.basename)                     # already created in model_data.py
domain.set_datadir(model_data.results_dir)                 # already created and/or checked in model_data.py
domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
domain.set_minimum_storable_height(0.01)                 # Remove very shallow water depths from sww
domain.set_maximum_allowed_speed(8.0)                    # Allow Maximum velocity of 8m/s....   

######################################################################################
#                                                                                    #
#                       ASSIGN THE DOMAIN INITIAL CONDITIONS                         #
#                                                                                    #
######################################################################################

#-------------------------------------------------------------------------------------
# This section assigns the initial (t=0) conditions for the domain
# Domain initial conditions can be fixed, read from a file, or from a function
# Note also that anuga permits almost all initial conditions to subsequently vary 
# with time, so while the following conditions apply to the whole domain at t=0
# all can be modified within the evolve time loop. In this way scour, deposition, 
# blockage, changes in roughness etc with time can be explicitly included in a simulation
#-------------------------------------------------------------------------------------

print model_data.basename+' >>>> Assigning initial conditions to %s domain' %(model_data.basename)

# Elevation data may be assigned to the domain from a csv or (NetCDF) pts file. 

# For a smaller csv that will be read in repeatedly, best to use geospatial_data to readin and 
# export_point_file to save the points permanently in a NetCDF *.pts file. This runs faster
# than the block read of a csv - but - is limited to about 6 million points.
# Alternatively - by directly assigning elevation from a csv file can take advantage of the block 
# read in domain.set_quantity so the elevation dataset can be virtually unlimited in size. 

# In this simulation the size of the elevation dataset is to large for conversion to a (NetCDF) 
# pts format so data is read in from the raw csv file. 

domain.set_quantity('elevation', 
                    filename = model_data.elev_filename,
                    use_cache=True,
                    verbose=model_data.anuga_verbose,
                    alpha=0.01)


# Set the intial water level (stage) across the domain              
InitialWaterLevel = model_data.InitialWaterLevel             # At t=0 water surface is in this model at initial lake level (0.3m)
if domain.get_quantity('elevation') < InitialWaterLevel:     # This can create ponds if land is below initial tailwater level
    domain.set_quantity('stage', InitialWaterLevel)          # set initial stage at initial tailwater level (wet)
else:        
    domain.set_quantity('stage', expression='elevation')     # else set initial stage at land level (dry)


# Assign uniform initial friction value to domain mesh centroids as initial condition
# Note: Friction will be re-assessed  based on roughness data read in from tuflow
# 2d_mat files and computed depth during evolve time loop!!!!
domain.set_quantity('friction', model_data.default_n, location = 'centroids')




#######################################################################################
#                                                                                     #
#             ASSIGN THE (TEMPORAL) SUBAREAL INFLOWS  TO THE DOMAIN                   #
#                                                                                     #
#######################################################################################

#--------------------------------------------------------------------------------------
# This section reads the local or total temporal inflows read in in model_data.py from
# the Tuflow ts1 files and stored as NETCDF tms file in the tms_files subdirectory,
# into the domain.
#
# Note: The current ANUGA inflow function 'pours' flow onto the domain within a circle 
# or polygonal area specified by the user at the specified temporal rate.
# EHR -- the above inflows should eventually be augmented by  inflows at the boundary 
# (total as distinct from local hydrographs) to better correspond with conventional 
# flood modelling.
#--------------------------------------------------------------------------------------

print model_data.basename+" >>>> Assigning the inflow hydrographs from the (%s) tms files directory" %(model_data.tms_dir)


inflow_fields=[]
# read in each line of inflow_hydrographs[tms_filename,xcoord,ycoord] and assign as inflow hydrograph to domain
for i in range(len(model_data.inflow_hydrographs)) :
    
    tms_filename = model_data.inflow_hydrographs[i][0]                   # extract the inflow hydro tms_filename
    xcoord = model_data.inflow_hydrographs[i][1]                         # extract the inflow hydro xcoord
    ycoord = model_data.inflow_hydrographs[i][2]                         # extract the inflow hydro ycoord
    print model_data.basename+' >>>> Reading inflow hydrograph from %s at %7.2f %7.2f and appending to forcing terms ' %(tms_filename,xcoord,ycoord)
        
    # convert time series to temporal function
    flowrate = file_function(tms_filename, 
                             quantities='value',
                             boundary_polygon=None)
    hydrograph = Inflow(domain,
                        center=(xcoord,ycoord),
                        radius=10, 
                        rate=flowrate)

    # append hydrograph to domain
    domain.forcing_terms.append(hydrograph)

print model_data.basename+' >>>> Completed assignment of %i inflow  hydrographs at t= %.2f hours ' %(len(model_data.inflow_hydrographs),float(time.time()-t0)/3600.0)


#######################################################################################
#                                                                                     #
#        ASSIGN THE (TEMPORAL) DOWNSTREAM BOUNDARY CONDITIONS TO THE DOMAIN           #
#                                                                                     #
#######################################################################################          


print model_data.basename+' >>>> Assigning the DSBC from the (%s) tms files directory ' %(model_data.tms_dir)
print model_data.basename+' >>>> Available boundary tags are ', domain.get_boundary_tags()

# convert dsbc time series to function of time 
tms_filename=model_data.dsbc_hydrographs[0]
lake_stage  = file_function(tms_filename,quantities='value',boundary_polygon=None)
# Note: will need to be modified if more than one dsbc is to be applied
print model_data.basename + ' >>>> Applying dsbc %s to model at the Lake_bdry ' %(tms_filename)

MaMt_bdry   = Reflective_boundary(domain)    # all reflective except lake as using inflow()
Mriv_bdry   = Reflective_boundary(domain)    # Really redundant untill get boundary inflows
Cbck_bdry   = Reflective_boundary(domain) 
WFra_bdry   = Reflective_boundary(domain) 
EFra_bdry   = Reflective_boundary(domain) 

lake_bdry   = Transmissive_Momentum_Set_Stage_boundary(domain=domain,function=lake_stage)
closed_bdry = Reflective_boundary(domain)                          

# boundary conditions for current scenario
domain.set_boundary({'Lake': lake_bdry,
                     'MaMt': MaMt_bdry,
                     'Mriv': Mriv_bdry,
                     'Cbck': Cbck_bdry,
                     'WFra': WFra_bdry,
                     'EFra': EFra_bdry,
                     'exterior': closed_bdry })  # this relates to the residue of the bounding poly not explicitly asigned above
    
print model_data.basename+' >>>> Completed assignment of %i dsbc boundary conditions to domain at t= %.2f hours ' %(len(model_data.dsbc_hydrographs),float(time.time()-t0)/3600.0)


#######################################################################################
#                                                                                     #
#  ASSIGN THE (TEMPORALLY DEPTH DEPENDENT) SURFACE ROUGHNESS (FRICTION) PARAMETERS    #
#                                                                                     #
#######################################################################################

#---------------------------------------------------------------------------------------
# This section creates an array of surface roughness data for each element in the domain 
# ready for later use in the evolve loop where the roughness data is  used with depth to
# compute a depth dependent friction at each yieldstep (output timestep to the swww file).
#
# Note that because the polys read in should but may not cover the entire bounding poly area
# It is prudent to initially declare a default material (as in Tuflow) applying the default 
# to the whole domain before applying (overwriting) with the spatial roughness data from the 
# polys obtained from the 2d_mat mid/mif files!!!!
#--------------------------------------------------------------------------------------

print model_data.basename+' >>>> Creating and assigning surface roughness data for use in evolve'

N = len(domain)
# Create the new roughness array to hold the elemental roughness data
# In this array of Nx5 values
#       n0 -- is the equivalent fixed roughness value  
#       d1 -- is the depth below which n1 applies
#       n1 -- is the roughness applying below d1
#       d2 -- is the depth above which n2 applies
#       n2 -- is the roughness applying above d2  

# Create and assign zeros to the roughness array
material_variables = zeros((N, 5), Float)                

# Set default material variables n0, d1, n1, d2, n2 in the roughness array
default_n0 = 0.040
default_d1 = 0.3
default_n1 = 0.050
default_d2 = 1.5
default_n2 = 0.030

material_variables[:,0] = default_n0   # n0
material_variables[:,1] = default_d1   # d1  
material_variables[:,2] = default_n1   # n1  
material_variables[:,3] = default_d2   # d2  
material_variables[:,4] = default_n2   # n2  

print '         ++++ Initially - material_variables data set to default values across whole domain ' 

# Update material_variables from model_data.mat_poly_list

# Obtain a list of real world centroid coords of model elements
points = domain.get_centroid_coordinates(absolute=True)   

for k, poly in enumerate(model_data.mat_poly_list):
    
    # Get indices of triangles inside polygon k
    indices = inside_polygon(points, poly, closed=True)

    # Index into material data for polygon k
    rid = model_data.mat_RID_list[k]
        
    # Material data for polygon     
    n0, d1, n1, d2, n2 = model_data.mat_roughness_data[rid]    
    
    # Update the data for the elements inside the mat_poly    
    for i in indices:
        material_variables[i, 0] = n0
        material_variables[i, 1] = d1        
        material_variables[i, 2] = n1
        material_variables[i, 3] = d2
        material_variables[i, 4] = n2                        
        
    

print '         ++++ Finally - material_variables data updated spatially with the material values set in the mif/mid files ' 

    
print model_data.basename+' >>>> Completed assignment of material (surface roughness) data to domain at t= %.2f hours ' %(float(time.time()-t0)/3600.0)

print model_data.basename+' >>>> Completed domain construction at t= %.2f hours ' %(float(time.time()-t0)/3600.0)



#########################################################################################
#                                                                                       #
#                                     RUN THE SIMULATION                                #
#                                                                                       #
#########################################################################################

#----------------------------------------------------------------------------------------
# This section starts the analysis of flow through the domain from the specified (realworld)
# starttime  up to the specified (realworld) finaltime
# All time is in seconds.                                                 
#
# Note:
# 1. -- The internal computational timestep is set at each timestep to meet a CFL=1 limit.
#       It is therefore important to examine the mesh statistics to see that there are no 
#       unintended small triangles, particularly in areas of deeper water as they can
#       dramatically reduce the computational timerstep and slow the model down.
# 2. -- The model runs on 'internal' time which will normally differ from realworld or
#       'Design' event time as input by the user in various datafiles. This 'internal'time is
#       normally not seen by the user but users should be aware of its existence.
# 3. -- In a similar fashion the model runs with all coordinates recomputed to a local
#       lower left frame origin (to impove numerical precission). The output files are in these
#       'internal' coordinates but each file also stores the new origin coordintes from which
#       the realworld coordinates can be restored (as presently occurs in animate.exe)
# 4. -- 'friction' is recomputed as a function of depth and location at each steptime for
#       all wet cells. Recomputation at every computational step was considered unnecessay
#       and greatly incresed run times.
#----------------------------------------------------------------------------------------
print model_data.basename+' >>>> Starting the simulation for ',model_data.catchment,'-',model_data.simclass,'-',model_data.scenario,'-',model_data.event
starttime     = 40*3600              # start simulation at t=40h relative to user/data time
endtime       = starttime + 4*3600   # end simulation at t=44h relative to user/data time
steptime      = 300                  # write sww  file results out every 300secs = 5min

print model_data.basename+' >>>> Simulation is for %.2f hours with output to the sww file at %.0f minute intervals ' %(float((endtime-starttime)/3600.0),float(steptime/60.0))

#  Initiate the starttime reset
domain.set_starttime(starttime)
print 'Start time', domain.get_time()


for t in domain.evolve(yieldstep = steptime, finaltime = endtime):       # This is a steptime loop (not computaional timestep)
    
    domain.write_time()                                                  # confirm t each steptime
    domain.write_boundary_statistics(tags='Lake', quantities='stage')    # indicate bdry stats each steptime
    depth = domain.create_quantity_from_expression('stage - elevation')  # create depth instance for this timestep

    print '         ++++ Recomputing friction at t= ', t

    # rebuild the 'friction' values adjusted for depth 
    for i in domain.get_wet_elements():                                  # loop for each wet element in domain
    
        # Get roughness variables 
        n0 = material_variables[i,0]
        d1 = material_variables[i,1]
        n1 = material_variables[i,2]
        d2 = material_variables[i,3]
        n2 = material_variables[i,4]                
                

        # Recompute friction values from depth for this element 
        d = depth.get_values(location='centroids', indices=[i])[0]
         
        if d <= d1: 
            depth_dependent_friction = n1
        elif d >= d2:
            depth_dependent_friction = n2                        
        else:
            depth_dependent_friction = n1+(n2-n1)/(d2-d1)*d                         


        domain.set_quantity('friction',
                            numeric=depth_dependent_friction,
                            location='centroids',
                            indices=[i],
                            verbose=model_data.anuga_verbose)                                
              
              
    if model_data.model_verbose :
        friction = domain.get_quantity('friction')
        
        # Print something
        
        

print model_data.basename+' >>>> ',model_data.catchment,'-',model_data.simclass,'-',model_data.scenario,'-',model_data.event,' -- Simulation completed succesfully'
print model_data.basename+' >>>> Completed ', endtime/3600.0, 'hours of model simulation at t= %.2f hours' %(float(time.time()-t0)/3600.0)