source: anuga_validation/Hinwood_2008/run_simulations.py @ 7660

"""

Script for running a breaking wave simulation of Jon Hinwood's wave tank.

Duncan Gray, GA - 2007 
"""


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

# Standard modules
import time
from os.path import join

# Related major packages
from anuga.shallow_water import Domain, Reflective_boundary, Time_boundary
from anuga.fit_interpolate.interpolate import interpolate_sww2csv
from anuga.abstract_2d_finite_volumes.util import file_function
from anuga.utilities.interp import interp

# Scenario specific imports
import project
import create_mesh
from prepare_time_boundary import csv2tms


def main(metadata_dic,
         maximum_triangle_area,
         yieldstep,
         boundary_path=None,
         friction=0.,  # planed wood 0.012 http://www.lmnoeng.com/manningn.htm
         outputdir_name=None,
         isTest=False,
         width=1.0,
         use_limits=True,
         end_tag = ''):
    """
    Run a simulation.
    """

    id = metadata_dic['scenario_id']
    
    if isTest is True:
        yieldstep = 1.0
        finaltime = 15.
        maximum_triangle_area=0.1
        outputdir_name += '_test'         
    else:
        finaltime = None
        
    outputdir_name = paras2outputdir_tag(mta=maximum_triangle_area,
                        yieldstep=yieldstep,
                        width=width,
                        use_limits=use_limits,
                        friction=friction,
                        end_tag=end_tag,
                        outputdir_name=outputdir_name)
    basename = outputdir_name       

    if finaltime is None:
        finaltime = metadata_dic['wave_times'][1] + 0.1 
        
    boundary_file_path = id + '_boundary.tsm'
        
    # Convert the boundary file, .csv to .tsm
    csv2tms(boundary_file_path, metadata_dic['xleft'][1])
   
    mesh_filename =  join(project.output_dir,basename + '.msh')

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

    create_mesh.generate(mesh_filename, metadata_dic, width=width,
                         maximum_triangle_area=maximum_triangle_area)
    
    #-------------------------------------------------------------------------
    # Setup computational domain
    #-------------------------------------------------------------------------
    domain = Domain(mesh_filename, use_cache = False, verbose = True)  

    print 'Number of triangles = ', len(domain)
    print 'The extent is ', domain.get_extent()
    print domain.statistics()
    
    domain.set_name(basename)
    domain.set_datadir(project.output_dir)
    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
    domain.set_minimum_storable_height(0.0001)

    if use_limits is True:
        domain.set_default_order(2) # Use second order spatial scheme
        domain.set_timestepping_method('rk2')
        domain.use_edge_limiter = True
        domain.tight_slope_limiters = True
        
        domain.beta_w      = 0.6
        domain.beta_uh     = 0.6
        domain.beta_vh     = 0.6
    

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

    domain.set_quantity('stage', 0.) #the z origin is the still water level
    domain.set_quantity('friction', friction)
    elevation_function = Elevation_function(metadata_dic)
    domain.set_quantity('elevation', elevation_function)

    function = file_function(boundary_file_path, domain, verbose=True)
        
    Br = Reflective_boundary(domain)
    Bts = Time_boundary(domain, function)
    domain.set_boundary( {'wall': Br, 'wave': Bts} )

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

    for t in domain.evolve(yieldstep, finaltime):    
        domain.write_time()
    print 'That took %.2f seconds' %(time.time()-t0)
    print 'finished'


    #-------------------------------------------------------------------------
    # Calculate gauge info
    #-------------------------------------------------------------------------

    if isTest is not True:
        points = []
        for gauge_x in metadata_dic['gauge_x']:
            points.append([gauge_x, 0.0])

        interpolate_sww2csv(join(project.output_dir,basename+".sww"),
                            points,
                            join(project.output_dir,basename + "_depth.csv"),
                            join(project.output_dir,basename + \
                                 "_velocit_x.csv"),
                           join(project.output_dir, basename + \
                                "_velocity_y.csv"),
                            join(project.output_dir,basename + "_stage.csv"))
  

def paras2outputdir_tag(mta,
                        yieldstep,
                        width,
                        use_limits,
                        friction,
                        end_tag,
                        outputdir_name=None):
    """
    Given run parameters return a string that can be used to name files.
    """
    if outputdir_name is None:
        outputdir_name = ''
        
    if use_limits is True:
        outputdir_name += '_lmts'
    else:
        outputdir_name += '_nolmts'
    outputdir_name += '_wdth_' + str(width)
    outputdir_name += '_z_' + str(friction)
    outputdir_name += '_ys_' + str(yieldstep)
    outputdir_name += '_mta_' + str(mta)
    outputdir_name += end_tag

    return outputdir_name


class Elevation_function:
    """
    
    Create a callable instance that returns the bed surface of the
    flume given the scenario metadata.
    
    """
    def __init__(self, metadata_dic):
        self.xslope_position = [metadata_dic['xleft'][0],
                                metadata_dic['xtoe'][0],
                                metadata_dic['xbeach'][0],
                                metadata_dic['xright'][0]]
        self.yslope_height = [metadata_dic['xleft'][1],
                              metadata_dic['xtoe'][1],
                              metadata_dic['xbeach'][1],
                              metadata_dic['xright'][1]]
        
    def __call__(self, x,y):       
        z = interp(self.yslope_height, self.xslope_position, x)
        return z 

#-------------------------------------------------------------
if __name__ == "__main__":

    # Import scenario metadata for the simulations that is to be run
    from scenarios import scenarios

    # Basic parameters of computer flume model (default: w=0.1, max=0.0001, ys=0.01, fr=0.0) 
    width = 0.1
    maximum_triangle_area=0.01
    yieldstep = 0.5
    friction=0.0

    # Select coarse fast run (isTest is True) or detailed slower run (isTest is False)
    #isTest=True
    isTest=False

    # Loop through experiments and run simulations
    for run_data in scenarios:
        main(run_data,
             maximum_triangle_area=maximum_triangle_area,
             yieldstep=yieldstep,
             width=width,
             isTest=isTest,
             outputdir_name=run_data['scenario_id'],
             use_limits=False,
             friction=friction,
             end_tag='_A')