source: anuga_work/production/tweed_heads/run_mesh_export_ascii.py @ 7479

Last change on this file since 7479 was 7479, checked in by miriam, 15 years ago

Scripts copied from Hobart, rewritten for tweed valley

File size: 5.1 KB
Line 
1"""This scripts recreates fitting problem where values seem to go towards zero
2in areas with missing data.
3"""
4from os.path import join
5from anuga.utilities.polygon import read_polygon
6from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf
7from anuga.shallow_water.data_manager import dem2pts
8from anuga.geospatial_data.geospatial_data import Geospatial_data
9from anuga.interface import create_domain_from_regions
10from anuga.interface import Reflective_boundary
11
12import project, os
13import sys
14from anuga.shallow_water.data_manager import start_screen_catcher
15from anuga.shallow_water.data_manager import copy_code_files
16from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
17from anuga.shallow_water.data_manager import sww2dem
18from os import sep
19
20# Application specific imports
21from setup_model import project   # Definition of file names and polygons
22
23
24#-------------------------------------------------------------------------------
25# Copy scripts to time stamped output directory and capture screen
26# output to file. Copy script must be before screen_catcher
27#-------------------------------------------------------------------------------
28copy_code_files(project.output_run, __file__, 
29                os.path.join(os.path.dirname(project.__file__),
30                             project.__name__+'.py'))
31start_screen_catcher(project.output_run, 0, 1)
32
33#-----------------------------------------------------------------------------
34# Preparation of topographic data
35#
36# Convert ASC 2 DEM 2 PTS using source data and store result in source data
37# Do for coarse and fine data
38# Fine pts file to be clipped to area of interest
39#-----------------------------------------------------------------------------
40print 'combined_elevation_basename', project.combined_elevation_basename
41
42#------------------------------------------------------------------------------
43# Create ANUGA mesh and fit
44#------------------------------------------------------------------------------
45tags = {'x': range(len(project.bounding_polygon))} # All segments tagged the same
46domain = create_domain_from_regions(project.bounding_polygon,
47                                    boundary_tags=tags,
48                                    maximum_triangle_area=project.bounding_polygon_maxarea,
49                                    interior_regions=project.interior_regions,
50                                    mesh_filename=project.meshes, 
51                                    use_cache=True,
52                                    verbose=True)
53print domain.statistics()
54
55
56domain.set_name(project.scenario_name)
57domain.set_datadir(project.output_run) 
58domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm
59
60#------------------------------------------------------------------------------
61# Setup initial conditions
62#------------------------------------------------------------------------------
63
64print 'Setup initial conditions'
65domain.set_quantity('stage', 0.0)
66domain.set_quantity('friction', 0.01)
67domain.set_quantity('elevation', 
68                    filename=project.combined_elevation+'.pts',
69                    use_cache=True,
70                    verbose=True,
71                    alpha=0.1)
72
73#------------------------------------------------------------------------------
74# Setup boundary conditions
75#------------------------------------------------------------------------------
76
77print 'Set boundary - available tags:', domain.get_boundary_tags()
78Br = Reflective_boundary(domain)
79domain.set_boundary({'x': Br}) 
80
81#-----------------------------------------------------------------------------
82# Take a few steps
83#-----------------------------------------------------------------------------
84
85
86for t in domain.evolve(yieldstep=1, 
87                       finaltime=3): 
88    print domain.timestepping_statistics()
89    print domain.boundary_statistics(tags='x')
90
91print 'Done'
92
93
94#-----------------------------------------------------------------------------
95# Exporting Grid
96# Define allowed variable names and associated equations to generate values.
97# This would not normally change.
98#-----------------------------------------------------------------------------
99var_equations = {'stage':     'stage',
100                 'momentum':  '(xmomentum**2 + ymomentum**2)**0.5',
101                 'depth':     'stage-elevation',
102                 'speed':     '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6)',
103                 'elevation': 'elevation' }
104
105# one or more key strings from var_equations above
106var = ['elevation', 'stage']
107
108######
109# Start running the various conversions we require.
110######
111
112for which_var in var:
113    if which_var not in var_equations:
114        print 'Unrecognized variable name: %s' % which_var
115        break
116   
117    name = join(output_run, scenario_name)
118    outname = name + '_' + which_var
119    quantityname = var_equations[which_var]
120
121    print 'start sww2dem: name=%s' % name
122   
123    sww2dem(name, basename_out = outname,
124                quantity = quantityname,
125                timestep = 0,
126                cellsize = 20,   
127                reduction = max, 
128                verbose = True,
129                format = 'asc')
130
131
132
Note: See TracBrowser for help on using the repository browser.