""" -----------------------------------------------------------------model_ini.py ANUGA MODELLING TEMPLATE (Model Initialisation) This initialisation script sets model run time flags and model-scenario-event specific parameters and file references for use in the subsequent model_xx.py scripts To run in parallel you must use the folowing command line mpirun -np n python model_run.py where n is the total number of processors being used for the parallel run. This script is imported by the model_data/build/checks/run/results scripts <<<< THIS SCRIPT CONTAINS MOST OF THE CODE THAT WILL NORMALLY NEED TO >>>> <<<< BE EDITED TO CREATE A PARTICULAR SIMULATION >>>> Version : 1.00 August 2015 flood template modified for tsunami modelling Author : E Rigby (0437 250 500) ted.rigby@rienco.com.au -------------------------------------------------------------------------------- """ ################################################################################# # # # IMPORT THE REQUIRED PYTHON MODULES # # # ################################################################################# # Standard python modules import os import sys from os import sep, mkdir from os.path import join, exists from time import time,localtime, strftime # Related ANUGA modules import anuga ################################################################################# # # # OPTIONALLY IDENTIFY THE USER & ORGANISATIONAL DETAILS # # An optional block included in case more corporate ID desired # ################################################################################# org_name = 'Rienco Consulting' org_add1 = '5 Pekham St' org_add2 = 'Sorell' org_add3 = 'Tas 7172' org_add4 = 'Australia' org_phone = '0437 250 500' org_web = 'www.rienco.com.au' user_name = 'Ted' user_email = 'ted.rigby@rienco.com.au' ################################################################################# # # # SET THE CORE SIMULATION DESCRIPTORS FOR THIS RUN # # # ################################################################################# model = 'SETAS' # Set model name #scenario = 'GA2009' # Set scenario scenario = 'MRT2016_Scenario4' # Set scenario event_number = '58280P' # 58260, 64477 the event number - See GA 2009 event_chosen.xls runtype = 'FINAL' # (TRIAL, BASIC, STABLE, FINAL) prefix = model+'_'+scenario+'_'+event_number+'_'+runtype # run ID model_verbose = True # set to True to enable additional model specific output to screen model_checks = True # set to True to enable additional model checks & plots to the checks dir screen_log = False # creates a logfile of screen IO if true # Note the anuga code uses final_time in evolve as duration rather than final_time # This will eventually be corrected and will require mods to the following once implemented starttime = 5922.5 # start of simulation = sts start time for this event (approx 1.6 hrs)) duration = 14400.0 # simulation duration of 14400 secs (4 hrs) (was 14hrs for MRT/GA2009) # note end simulation after a duration of 14400 secs (at approx 5.6 hrs) # note yieldstep is set later with the run resolution reflecting coarse or fine simulations ################################################################################# # # # OPTIONALLY COPY SCREEN I-O TO LOGFILE & IDENTIFY BUILD # # # ################################################################################# # create the build identification of this run user = anuga.get_user_name() # system names used in file stamping host = anuga.get_host_name() # ditto t0 = time() # record when the run started in secs build_time = strftime('%Y%m%d_%H%M%S',localtime()) build_ID = build_time+'_build'+'_'+str(user)+'_on_'+str(host) anuga_revision = anuga.get_revision_number() anuga_ID = str(anuga_revision) print ' ' print prefix+ ' >>>> #########################################################################' print prefix+ ' >>>>' print prefix+ ' >>>> COMMENCING THIS MODEL SIMULATION USING ANUGA REVISION NUMBER %s ' %(anuga_ID) print prefix+ ' >>>>' print prefix+ ' >>>> for the simulation named ' print prefix+ ' >>>> ', prefix print prefix+ ' >>>> and model Build ID of ' print prefix+ ' >>>> ', build_ID print prefix+ ' >>>> #########################################################################\n' print ' ' ################################################################################# # # # SET THE REQUIRED ANUGA SPECIFIC RUN TIME FLAGS & VARIABLES # # # ################################################################################# minimum_allowed_height = 0.001 # sets min depth (m) recognised in calcs minimum_storable_height = 0.010 # writes d=0 (m)to sww if d < 0.010 - improves graphics maximum_allowed_speed = 5.00 # limits speed to 5m/sec max if d < minimum_allowed_height anuga_verbose = True # set to True to enable additional anuga internal output flow_algorithm = 'DE0' # set flow algorithm - refer anuga..config.py for options # NOTES RE Flow_algorithm selection # DE0 used Euler timestepping (and so is first order accurate in time). # DE1 uses the average of two Euler steps to get a second order time stepping method. # DE2 uses the average of three Euler steps to get a third order time stepping method. # DE0 is very diffusive, so if you have a periodic forcing, the height of the wave will # decrease quite markedly. So not so good for complicated tsunami forcing problems. # DE0 works quite well for flood modelling. # Testing at GA suugests DE2 currently provides the best replication of tsunami behaviour # DE2 will typically take 3 times as long as DE0 and DE1 takes 2 times as long as DE0. # Reccomend use of DE0 generally for speeed and switch to DE2 for final production runs # when modelling a tsunami. ################################################################################# # # # SET THE SIMULATION SPECIFIC RUN TIME FLAGS & VARIABLES # # # ################################################################################# # assign names for the mesh and domain for this simulation mesh_filename = prefix+'.tsh' # msh for binary out and tsh for ascii output domain_name = prefix ##### Assign the different internal mesh resolutions/yieldsteps to be used in this simulation # Set the various model resolutions to suit the run type if runtype == 'TRIAL' : # only run with a extremely coarse mesh mesh_factor = 64.0 yieldstep = 600 # 10min elif runtype == 'BASIC' : # run with a relatively coarse mesh mesh_factor = 16.0 yieldstep = 240 # 6min elif runtype == 'STABLE' : # run with reasonably fine mesh mesh_factor = 4.0 yieldstep = 120 # 2min else : # default is FINAL with unfactored res/yieldstep as below mesh_factor = 1.0 yieldstep = 60 # 1min ##### Assign the resolution class to factored resolution values for use in model_data # As used for priority areas in GA2009 # vfine_res = mesh_factor * 500.0 # m2 vrfine limit 32x32m tri (typ high detail built up areas) # As used for aos1 and aos2 areas in GA2009 # fine_res = mesh_factor * 1500.0 # m2 fine limit 55x55m tri (typ areas of interest generally) # As used for sw area in GA2009 # medium_res = mesh_factor * 30000.0 # m2 medium limit 245x245m tri (typ coastline generally and large estuaries) # As used for default bounding poly res in GA2009 # coarse_res = mesh_factor * 100000.0 # m2 coarse limit 450x450m tri (typ open deepwater areas) # As used for resolution areas in MRT2016 xfine_res = 50 * mesh_factor # m2 10x10m triangles in high res plot area (Blackmans Bay) vfine_res = 200 * mesh_factor # m2 20x20m triangles in std res plot areas of interest fine_res = 1250 * mesh_factor # m2 50x50m triangles in coastal strip -10 to 10m AHD medium_res = 20000 * mesh_factor # m2 200x200m triangles in moderate depth ocean -10 to -30m AMD coarse_res = 80000 * mesh_factor # m2 400x400m triangles in deeper water -30m and deeper default_res = coarse_res # applies to all bounding_poly pre internal res being added ##### Assign the initial default n to apply everywhere prior to applying spatial ns default_n = 0.010 # reflects estuary bed n - changes for land etc will be stamped over ##### Assign initial water level for this simulation HAT = 0.8 # SETAS MSL = 0.0 # SETAS LAT = -0.0 # SETAS # tide = xxx # could be a timeseries initial_water_level = HAT ################################################################################ # # # ASSIGN AND CHECK THAT THE REQUIRED SIMULATION DIRECTORIES EXIST # # # ################################################################################ #------------------------------------------------------------------------------- # Checks the current model directories required to run the simulation are present # This assumes a directory structure in which the model run scripts, model data files, # model check files and model result files are in a directory called event - as outlined below. # # model\scenario\event\data ( eg elevation, roughness, polygons, ICs,) # \ Event \ boundary timeseries # \checks ( eg screen logs, check files and plots) # \scripts ( py ini, data, build, check, run, results scripts) # \results (sww result file) # \ timeseries (timeseries files - (csv and plot files of eg stage H -v- time) # \ spatial (spatial eg asc max surface grid files) # # # eg # \SETAS2015\GA2009\58260\data \ 58260 \ various timeseries files # \ various data files # \ checks \ various check files # \ scripts \ python files # \ results \ output.sww file # \ timeseries\ various timeseries files # \ spatial \ various spatial files # #NOTE : \ scripts is the current directory when model_run.py is run # all directories are assigned relative to the scripts directory #------------------------------------------------------------------------------- ##### Assign the required input directories # ANUGA data data (input) directory. All input data is read from this directory data_dir = join('..' ,'data') resolution_dir = join(data_dir, 'ResolutionPolygons2016_5revised') gauges_dir = data_dir+sep+'gauges' # ANUGA results (output sww) directory and parent of timeseries ansd spatial output dirs results_dir = join('..','results') # ANUGA results/timeseries (output) directory timeseries_dir = results_dir+sep+'timeseries' # ANUGA output spatial output files directory spatial_dir = results_dir+sep+'spatial' # ANUGA output check files directory checks_dir = join('..','checks') ##### Check the required output directories are present and create if not if not exists(results_dir) : # the required results directory does not exist mkdir(results_dir) print prefix+ ' >>>> NOTE - Results directory %s was missing and has been created ' %(results_dir) if not exists(checks_dir) : # the required checks directory does not exist mkdir(checks_dir) print prefix+ ' >>>> NOTE - checks directory %s was missing and has been created ' %(checks_dir) if not exists(timeseries_dir) : # the required time series results directory does not exist mkdir(timeseries_dir) print prefix+ ' >>>> NOTE - Timeseries directory %s was missing and has been created ' %(timeseries_dir) if not exists(spatial_dir) : # the required spatial results directory does not exist mkdir(spatial_dir) print prefix+ ' >>>> NOTE - Spatial directory %s was missing and has been created ' %(spatial_dir) ################################################################################ # # # ASSIGN THE VARIOUS INPUT FILES FOR THIS SIMULATION # # # ################################################################################ ##### Assign the backdrop coastline polygon filename for model check plots # Note - can not use shp for bounding poly as order in sts bdry points is # counterclockwise contray to GIS rules. coastal_poly_filename = 'setas_coastal_polygon.shp' ##### Assign the bounding polygon for this simulation # Note - can not use shp for bounding poly as order in sts bdry points is # counterclockwise contrary to GIS rules. Need to re-order sts if want to # use shp to input bounding poly in future! bounding_poly_filename = 'bounding_polygon_GA2009_sts.csv' ##### Assign the internal polygon filenames for the MRT 2016 simulation # Internal resolution polygons for MRT2016 - coarse is default set in bounding poly int1_res_poly_filename = 'OceanApproachesNested.shp' # medium res - ocean ring 1 30m deep int2_res_poly_filename = 'CoastalZoneNested.shp' # fine res - coastal ring 2 10m deep - nested in ocean zone # the CoastalZoneAbove10m is added in at this point to allow larger triangles in the higher ground of islands that # are otherwise all at coastal res from Anuga's need to have only simple polys (no holes) int3_res_poly_filename = 'CoastalZoneAbove10mNested.shp' # coarse res - land ring 3a > 10m AHD - nested in coastal zone int4_res_poly_filename = 'PlotAreas.shp' # vfine res - plot areas merged ting 3b - nested in coastal zone int5_res_poly_filename = 'StudyArea.shp' # High res - study area ring 4 - nested in plot areas # same order as above ie int1_res_poly_filename polys will be modelled at int1_res int1_res = medium_res int2_res = fine_res int3_res = coarse_res int4_res = vfine_res int5_res = xfine_res ##### Assign the ANUGA input combined elevation spatial dataset filename for this simulation (xyz,shp or asc) elevation_points_filename = 'combined_elevation_MRT2016.txt' # assign the combined land/sea elevation points name ##### Assign the ANUGA input surface roughness (n) spatial dataset filename for this simulation mannings_grid_filename = 'manningsv2_20160302.asc' # assign the combined land/sea roughness shapefile name ##### Assign the ANUGA input boundary conditions timeseries for this simulation event_sts_filename = 'hobart.sts' # a preprepared netcdf of stage timeseries at various locations ##### Assign the ANUGA gauge locations dataset filename at which time-series output from this simulation is required gauges_filename = 'time_of_arrival_all.csv' # assign the gauges (points) shp filename ################################################################################# # # # CREATE AND CHECK THE REQUIRED PATHNAMES ARE PRESENT # # # ################################################################################# # create the full pathname for the output msh file mesh_pathname = join(checks_dir,mesh_filename) # dump the mesh into the checks directory # create the full pathname for coastal_filename and check exists coastal_poly_pathname = join(data_dir,coastal_poly_filename) if not exists(coastal_poly_pathname) : # check and stop run if absent print prefix+ ' >>>> ERROR - specified coastal poly file %s missing ' %(coastal_poly_pathname) sys.exit(0) # terminate the program # create the full pathname for bounding_poly_filename and check exists bounding_poly_pathname = join(data_dir,bounding_poly_filename) if not exists(bounding_poly_pathname) : # check and stop run if absent print prefix+ ' >>>> ERROR - specified bounding poly file %s missing ' %(bounding_poly_pathname) sys.exit(0) # terminate the program # create the full pathname for the 5int res polys and check exists int1_res_poly_pathname = join(data_dir,int1_res_poly_filename) # Coarsest of the interior polys first int2_res_poly_pathname = join(data_dir,int2_res_poly_filename) int3_res_poly_pathname = join(data_dir,int3_res_poly_filename) int4_res_poly_pathname = join(data_dir,int4_res_poly_filename) int5_res_poly_pathname = join(data_dir,int5_res_poly_filename) if not exists(int1_res_poly_pathname) : # check and stop run if absent print prefix+ ' >>>> ERROR - specified int res poly file %s missing ' %(int1_res_poly_pathname) sys.exit(0) if not exists(int2_res_poly_pathname) : # check and stop run if absent print prefix+ ' >>>> ERROR - specified int res poly file %s missing ' %(int2_res_poly_pathname) sys.exit(0) if not exists(int3_res_poly_pathname) : # check and stop run if absent print prefix+ ' >>>> ERROR - specified int res poly file %s missing ' %(int3_res_poly_pathname) sys.exit(0) if not exists(int4_res_poly_pathname) : # check and stop run if absent print prefix+ ' >>>> ERROR - specified int res poly file %s missing ' %(int4_res_poly_pathname) sys.exit(0) if not exists(int5_res_poly_pathname) : # check and stop run if absent print prefix+ ' >>>> ERROR - specified int res poly file %s missing ' %(int5_res_poly_pathname) sys.exit(0) # Create the full pathname for the elevation points dataset and check exists elevation_points_pathname = join (data_dir, elevation_points_filename) if not exists(elevation_points_pathname) : # check and stop run if absent print prefix+ ' >>>> ERROR - specified elevation file %s missing ' %(elevation_points_pathname) sys.exit(0) # Create the full pathname for the mannings roughness polys and check exists mannings_grid_pathname = join(data_dir,mannings_grid_filename) if not exists(mannings_grid_pathname) : # check and stop run if absent print prefix+ ' >>>> ERROR - specified mannings roughness file %s missing ' %(mannings_grid_pathname) sys.exit(0) # Create the full pathname for the netcdf sts boundary conditions file and check exists event_sts_pathname = join(data_dir,event_sts_filename) if not exists(event_sts_pathname) : # check and stop run if absent print prefix+ ' >>>> ERROR - specified boundary conditions file %s missing ' %(event_sts_pathname) sys.exit(0) # Create the full pathname for the file defining the output gauge lo0cations and check exists gauges_pathname = join(gauges_dir,gauges_filename) if not exists(gauges_pathname) : # check and stop run if absent print prefix+ ' >>>> ERROR - specified Gauges file %s missing ' %(gauges_pathname) sys.exit(0)