"""Script for running a dam break simulation of UQ's dam break tank.


Ole Nielsen and Duncan Gray, GA - 2006 
"""


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

# Standard modules
import time
import sys
from shutil import copy
from os import path

# Related major packages
import anuga

from anuga.abstract_2d_finite_volumes.region import Set_region 
from anuga.fit_interpolate.interpolate import interpolate_sww2csv, \
     file_function
from anuga.utilities.file_utils import copy_code_files
from anuga.file.csv_file import load_csv_as_dict
from numerical_tools import  err  

# Application specific imports
import create_mesh
import project

def main():     
    slope= 0
    friction = 0.01 
    inital_depth = 0.2
    gate_position = 0.75
   
    return scenario(slope, friction, inital_depth, gate_position)
            

def scenario(slope, friction, inital_depth, gate_position):

    #-------------------------------------------------------------------------
    # Create the triangular mesh
    #-------------------------------------------------------------------------

    create_mesh.generate(project.mesh_filename,
                         gate_position,
                         #is_course=True) # this creates the mesh
                         is_course=False) # this creates the mesh

    head,tail = path.split(project.mesh_filename)
    
    #-------------------------------------------------------------------------
    # Setup computational domain
    #-------------------------------------------------------------------------
    domain = anuga.Domain(project.mesh_filename, use_cache = False)
    domain.set_name(project.basename)
    domain.set_datadir('.')
    domain.set_quantities_to_be_stored({'stage':2,
                                        'xmomentum': 2,
                                        'ymomentum': 2,
                                        'elevation': 1})
    domain.set_minimum_storable_height(0.001)
    domain.set_store_vertices_uniquely(True)  # for writting to sww

    #-------------------------------------------------------------------------
    # Setup initial conditions
    #-------------------------------------------------------------------------


    def elevation_tilt(x, y):
        return x*slope
        
    domain.set_quantity('stage', elevation_tilt)
    domain.set_quantity('friction', friction) 
    domain.set_quantity('elevation',elevation_tilt)

    domain.set_region('dam','stage',inital_depth,
                                 location = 'unique vertices') 

    Br = anuga.Reflective_boundary(domain)
    Bd = anuga.Dirichlet_boundary([0,0,0])  # to drain the water out.
    domain.set_boundary( {'wall': Br, 'edge': Bd} )

    #-------------------------------------------------------------------------
    # Evolve system through time
    #-------------------------------------------------------------------------
    t0 = time.time()

    for t in domain.evolve(yieldstep = 0.1, finaltime = 10):
        pass

    # Load actual experimental results
    actual,title_index_dic = load_csv_as_dict(project.actual_filename)
    
    gauge_locations = [[0.4,0.2]]
    quantities = ['stage', 'elevation']
    file_instance = file_function( project.basename +".sww",
				      quantities = quantities,
				      interpolation_points = gauge_locations,
				      verbose = False,
				      use_cache = False)
    # create a list of the simulated_depths at the actual data times. 
    simulated_depths = []
    for atime in actual['time']:
        quantities_slice = file_instance(float(atime),
                                         point_id=0)
        depth = quantities_slice[0] - quantities_slice[1]
        simulated_depths.append(depth)
    flume_depths = actual["0.4:0.2"]
    flume_depths = [float(i) for i in flume_depths]
    # calc the norm
    #print "map(None, simulated_depths, flume_depths)", \
    #      map(None, simulated_depths, flume_depths)
    norm = err(simulated_depths,
               flume_depths, 2, relative = True)  # 2nd norm (rel. RMS
    return norm
#-------------------------------------------------------------
if __name__ == "__main__":
    main()