# -*- coding: utf-8 -*- #------------------------------------------------------------------------------ # Import necessary modules #------------------------------------------------------------------------------ # Standard modules import os import time import sys import numpy # Related major packages import anuga from anuga.file_conversion.file_conversion import timefile2netcdf from anuga import Inlet_operator # Application specific imports import project # Definition of file names and polygons time00 = time.time() #------------------------------------------------------------------------------ # Preparation of topographic data # Convert ASC 2 DEM 2 PTS using source data and store result in source data #------------------------------------------------------------------------------ # Create DEM from asc data anuga.asc2dem(project.name_stem+'.asc', use_cache=project.cache, verbose=project.verbose) # Create pts file for DEM anuga.dem2pts(project.name_stem+'.dem', use_cache=project.cache, verbose=project.verbose) #------------------------------------------------------------------------------ # Create the triangular mesh and domain based on # overall clipping polygon with a tagged # boundary as defined in project.py #------------------------------------------------------------------------------ domain = anuga.create_domain_from_regions(project.bounding_polygon, boundary_tags={'e0': [0], 'e1': [1], 'e2': [2], 'e3': [3], 'e4': [4], 'e5': [5], 'e6': [6], 'e7': [7], 'e8': [8], 'e9': [9], 'e10': [10], 'e11': [11], 'e12': [12], 'e13': [13], 'e14': [14], 'e15': [15], 'e16': [16], 'e17': [17], 'e18': [18], 'e19': [19], 'e20': [20], 'e21': [21]}, maximum_triangle_area=project.default_res, mesh_filename=project.meshname, use_cache=project.cache, verbose=project.verbose) # Print some stats about mesh and domain print 'Number of triangles = ', len(domain) print 'The extent is ', domain.get_extent() print domain.statistics() #------------------------------------------------------------------------------ # Setup parameters of computational domain #------------------------------------------------------------------------------ domain.set_name('spring_' + project.scenario) # Name of sww file domain.set_datadir('.') # Store sww output here domain.set_minimum_storable_height(0.01) # Store only depth > 1cm domain.set_flow_algorithm('DE0') #------------------------------------------------------------------------------ # Setup initial conditions #------------------------------------------------------------------------------ tide = project.tide domain.set_quantity('stage', tide) #domain.set_quantity('friction', 0.0) domain.set_quantity('friction', 0.05) domain.set_quantity('elevation', filename=project.name_stem + '.pts', use_cache=project.cache, verbose=project.verbose, alpha=0.1) time01 = time.time() print 'That took %.2f seconds to fit data' %(time01-time00) if project.just_fitting: import sys sys.exit() #------------------------------------------------------------------------------ # Setup boundary conditions #------------------------------------------------------------------------------ print 'Available boundary tags', domain.get_boundary_tags() Bd_outlet = anuga.Dirichlet_boundary([14.0, 0, 0]) # Use a fixed water level at the basin outlet Bd = anuga.Dirichlet_boundary([tide, 0, 0]) # A water level lower than elevation domain.set_boundary({'e0': Bd_outlet, 'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd, 'e5': Bd, 'e6': Bd, 'e7': Bd, 'e8': Bd, 'e9': Bd, 'e10': Bd, 'e11': Bd, 'e12': Bd, 'e13': Bd, 'e14': Bd, 'e15': Bd, 'e16': Bd, 'e17': Bd, 'e18': Bd, 'e19': Bd, 'e20': Bd, 'e21': Bd}) #------------------------------------------------------------------------------ # Rate operator #------------------------------------------------------------------------------ #from anuga.operators.rate_operators import Rate_operator #surR = anuga.Quantity(domain) #factor = 1.0 / 1000.0 / 3600.0 # from mm/hr to m/sec #surR_op = anuga.Rate_operator(domain, rate = surR, factor=factor) # inlet cross-sections: #line0 = [[244326,3335740], [245821,3333810]] #line1 = [[239732,3324250], [239766,3323140]] #line2 = [[237899,3317860], [237416,3316920]] # inlet polygons: poly0 = [[244326,3335740], [245026,3335740], [245821,3334510], [245821,3333810]] poly1 = [[239732,3324250], [239965,3324250], [239999,3323140], [239766,3323140]] poly2 = [[237899,3317860], [238199,3317860], [237716,3316920], [237416,3316920]] # obs streamflow: strmflowFile0 = 'USGS 08068275 Spring Ck nr Tomball_0826-0830_datetime.txt' strmflowFile1 = 'USGS 08068780 Little Cypress Ck nr Cypress_0826-0830_datetime.txt' strmflowFile2 = 'USGS 08068740 Cypress Ck at House-Hahl Rd nr Cypress_0826-0830_datetime.txt' # convert to netCDF: timefile2netcdf(strmflowFile0, quantity_names=["strmflow"], time_as_seconds=False) timefile2netcdf(strmflowFile1, quantity_names=["strmflow"], time_as_seconds=False) timefile2netcdf(strmflowFile2, quantity_names=["strmflow"], time_as_seconds=False) # netCDF file name: tmsFile0 = strmflowFile0[:-4] + ".tms" tmsFile1 = strmflowFile1[:-4] + ".tms" tmsFile2 = strmflowFile2[:-4] + ".tms" Q0 = anuga.file_function(tmsFile0, quantities=['strmflow']) Q1 = anuga.file_function(tmsFile1, quantities=['strmflow']) Q2 = anuga.file_function(tmsFile2, quantities=['strmflow']) inlet0 = None inlet1 = None inlet2 = None # inlet operator: #inlet0 = Inlet_operator(domain, line0, Q0, logging=True, description='inlet0', verbose = False) #inlet1 = Inlet_operator(domain, line1, Q1, logging=True, description='inlet1', verbose = False) #inlet2 = Inlet_operator(domain, line2, Q2, logging=True, description='inlet2', verbose = False) inlet0 = Inlet_operator(domain, poly0, Q0, logging=True, description='inlet0', verbose = False) inlet1 = Inlet_operator(domain, poly1, Q1, logging=True, description='inlet1', verbose = False) inlet2 = Inlet_operator(domain, poly2, Q2, logging=True, description='inlet2', verbose = False) #if inlet0 is not None and verbose: inlet0.print_statistics() if inlet0 is not None: inlet0.print_statistics() #if inlet1 is not None and verbose: inlet1.print_statistics() if inlet1 is not None: inlet1.print_statistics() #if inlet2 is not None and verbose: inlet2.print_statistics() if inlet2 is not None: inlet2.print_statistics() #------------------------------------------------------------------------------ # Evolve system through time #------------------------------------------------------------------------------ time02 = time.time() #dir_dem = '/scratch/luoxy/Rsch_Flooding2/EF5_2/output_asc_utm/' # directory of surface-runoff "dem" files #listing = os.listdir(dir_dem) #listing.sort() #surRfileList = [] # list for storing surface-runoff "dem" files #for file in listing: # if file.endswith('.dem'): # dirFile = os.path.join(dir_dem, file) # surRfileList.append(dirFile) yieldstep = 900.0 # interval between outputs to screen (unit: second) finaltime = 432000.0 # end of simulation (= 5 days) (unit: second) for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime): # t starts from 0 #stepSize = 600.0 # unit: second #its = int( (t + 0.001) // stepSize ) # index of yieldstep (starts from 0) ('//' means floor division, e.g.: 4.0//3.0 = 1.0) #current_surR_file = surRfileList[ its ] # current surface-runoff "dem" file (list index starts from 0) #print 'current_surR_file:', current_surR_file #surR.set_values_from_utm_grid_file(current_surR_file, location='centroids') #print 'surR.centroid_values=', surR.centroid_values domain.print_timestepping_statistics() #domain.print_operator_timestepping_statistics() stage = domain.get_quantity('stage') elev = domain.get_quantity('elevation') height = stage - elev # water depth #print ' Integral water depth (m) = ', height.get_integral() print ' Integral water volume = ', height.get_integral() # --------------------------------- time03 = time.time() print 'That took %.2f seconds to simulate.' %(time03 - time02)