source: anuga_work/development/convergence_okushiri_2008/run_okushiri_truescale.py @ 5442

Last change on this file since 5442 was 5442, checked in by ole, 16 years ago

Retired h-limiter and beta_h as per ticket:194.
All unit tests and validation tests pass.

File size: 4.9 KB
RevLine 
[5411]1"""Anuga convergence study using a true-scale version of the Okushiri
2Island tsunami wavetank experiment
[5343]3
[5411]4This script runs a modified version of the Okushiri Island benchmark
[5343]5
6as published at
7
8THE THIRD INTERNATIONAL WORKSHOP ON LONG-WAVE RUNUP MODELS
9June 17-18 2004
10Wrigley Marine Science Center
11Catalina Island, California
12http://www.cee.cornell.edu/longwave/
13
14This version up-scales the original 1:400 scale wave-tank experiment, to
[5411]15"true-scale" with mesh defined using interior polygons.
[5343]16
[5411]17The original validation data is available at
[5343]18http://www.cee.cornell.edu/longwave/index.cfm?page=benchmark&problem=2
19where a detailed description of the problem is also available.
20
[5411]21Run create_okushiri_truescale.py to convert bathymetry and input boundary
22condition into NetCDF format.
[5343]23
24"""
[5411]25#----------------------------------------
26# Import necessary modules
27#----------------------------------------
[5343]28
[5411]29# Standard modules
30from os import sep,umask
31from os.path import dirname, basename
32from os import mkdir, access, F_OK
33from shutil import copy
34import time
35import sys
36
37# Related major packages
[5343]38from anuga.shallow_water import Domain
39from anuga.shallow_water import Reflective_boundary
40from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
41from anuga.abstract_2d_finite_volumes.util import file_function
[5411]42from anuga.shallow_water.data_manager import copy_code_files, start_screen_catcher
43from anuga.pmesh.mesh_interface import create_mesh_from_regions
[5343]44
45
[5411]46import project_truescale # Definition of filenames and interior polygons
[5343]47
[5411]48copy_code_files(project_truescale.output_run_time_dir,'run_okushiri_truescale.py',
49                'project_truescale.py' )
50myid = 0
51numprocs = 1
52start_screen_catcher(project_truescale.output_run_time_dir, myid, numprocs)
53
54#--------------------------------------------------------------------------
55# Create the triangular mesh based on overall bounding polygon with a
56# tagged boundary and interior regions defined in project_truescale.py
57# along with resolutions (maximal area of per triangle) for each polygon
58#--------------------------------------------------------------------------
59
60create_mesh_from_regions(project_truescale.poly_all,
61                         boundary_tags={'wall': [0, 1, 3],'wave': [2]},
62                         maximum_triangle_area=project_truescale.res_poly_all,
63                         #interior_regions=project_truescale.interior_regions,
64                         filename=project_truescale.mesh_name+'.msh',
65                         verbose=True)
66
67#------------------------------
[5343]68# Create Domain from mesh
[5411]69#------------------------------
70domain = Domain(project_truescale.mesh_name+'.msh', use_cache=True, verbose=True)
[5343]71print domain.statistics()
72
[5396]73z = domain.get_vertex_coordinates()
[5343]74
[5396]75# Write vertex coordinates to file
[5411]76filename=project.vertex_filename
[5396]77fid=open(filename,'w')
78fid.write('x (m), y (m)\n')
79for i in range(len(z)):
80    pt=z[i]
81    x=pt[0]
82    y=pt[1]
83    fid.write('%.6f, %.6f\n' %(x, y))
84       
[5343]85#-------------------------
86# Initial Conditions
87#-------------------------
88domain.set_quantity('friction', 0.0)
89domain.set_quantity('stage', 0.0)
90domain.set_quantity('elevation',
91                    filename=project_truescale.bathymetry_filename,
92                    alpha=0.02,                   
93                    verbose=True,
94                    use_cache=True)
95
96
97#-------------------------
98# Set simulation parameters
99#-------------------------
100domain.set_name(project_truescale.output_filename)  # Name of output sww file
[5411]101domain.set_datadir(project_truescale.output_run_time_dir) # Name of output directory
[5343]102domain.set_default_order(2)               # Apply second order scheme
103domain.set_all_limiters(0.9)              # Max second order scheme (old lim)
[5411]104domain.set_minimum_storable_height(0.1)   # Don't store h < 0.1m
[5343]105
106
107#-------------------------
108# Boundary Conditions
109#-------------------------
110
111# Create boundary function from timeseries provided in file
112function = file_function(project_truescale.boundary_filename,
113                         domain, verbose=True)
114
115# Create and assign boundary objects
116Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
117Br = Reflective_boundary(domain)
118domain.set_boundary({'wave': Bts, 'wall': Br})
119
120
121# Select triangle containing ch5 for diagnostic output
122# around known gauge
123triangle_id = domain.get_triangle_containing_point([1808.4, 478.4])
124# This should get triangle id 32833 with centroid (4.5244, 1.1972)
125
126
127#-------------------------
128# Evolve through time
129#-------------------------
130import time
131t0 = time.time()
132
[5411]133for t in domain.evolve(yieldstep=project_truescale.yieldstep, finaltime = 450):
[5343]134    print domain.timestepping_statistics(track_speeds=False,
135                                         triangle_id=triangle_id)
136    print domain.boundary_statistics()
137
138print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.