source: anuga_validation/okushiri_2005/run_okushiri_parallel.py @ 5847

Last change on this file since 5847 was 5847, checked in by steve, 15 years ago

Changed parallel_api so that global mesh only needs to
be constructed on processor 0

File size: 3.6 KB
Line 
1"""Validation of the AnuGA implementation of the shallow water wave equation.
2
3This script sets up Okushiri Island benchmark as published at the
4
5THE THIRD INTERNATIONAL WORKSHOP ON LONG-WAVE RUNUP MODELS
6June 17-18 2004
7Wrigley Marine Science Center
8Catalina Island, California
9http://www.cee.cornell.edu/longwave/
10
11
12The validation data was downloaded and made available in this directory
13for convenience but the original data is available at
14http://www.cee.cornell.edu/longwave/index.cfm?page=benchmark&problem=2
15where a detailed description of the problem is also available.
16
17
18Run create_okushiri.py to process the boundary condition and build a the
19mesh before running this script.
20
21"""
22
23# To know quickly if pypar can be imported
24import imp
25print " imp.find_module('pypar')", imp.find_module('pypar')
26import  pypar
27
28# Module imports
29from anuga.shallow_water import Domain
30from anuga.shallow_water import Reflective_boundary
31from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
32from anuga.abstract_2d_finite_volumes.util import file_function
33
34from anuga_parallel.parallel_api import myid, numprocs, distribute   
35
36import project
37
38use_cache = True
39#-------------------------
40# Create Domain from mesh
41# on processor 0
42#-------------------------
43if myid == 0 :
44    domain = Domain(project.mesh_filename, use_cache=use_cache, verbose=True)
45    print domain.statistics()
46
47    #-------------------------
48    # Initial Conditions
49    # At present need to do the
50    # fitting of the bathymetry
51    # on a global domain, ie before
52    # distributing the domain
53    #-------------------------
54    domain.set_quantity('friction', 0.0)
55    domain.set_quantity('stage', 0.0)
56    domain.set_quantity('elevation',
57                        filename=project.bathymetry_filename,
58                        alpha=0.02,                   
59                        verbose=True,
60                        use_cache=use_cache)
61else:
62    domain = None
63
64#-------------------------
65# Distribute domain if run in parallel
66#-------------------------
67if numprocs > 1:
68    domain = distribute(domain)
69
70
71#-------------------------
72# Set simulation parameters
73#-------------------------
74domain.set_name(project.output_filename)  # Name of output sww file
75domain.set_default_order(2)               # Apply second order scheme
76domain.set_all_limiters(0.9)              # Max second order scheme (old lim)
77domain.set_minimum_storable_height(0.001) # Don't store w < 0.001m
78domain.set_maximum_allowed_speed(0.1)     # Allow a little runoff (0.1 is OK)
79
80
81#-------------------------
82# Boundary Conditions
83#-------------------------
84
85# Create boundary function from timeseries provided in file
86function = file_function(project.boundary_filename,
87                         domain, verbose=True)
88
89# Create and assign boundary objects
90Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
91Br = Reflective_boundary(domain)
92domain.set_boundary({'wave': Bts, 'wall': Br})
93
94# Select triangle containing ch5 for diagnostic output
95# around known gauge
96# This should get triangle id 32833 with centroid (4.5244, 1.1972) on one processor
97try: 
98    triangle_id = domain.get_triangle_containing_point([4.521, 1.196])
99except:
100    triangle_id = None
101
102
103#-------------------------
104# Evolve through time
105#-------------------------
106import time
107t0 = time.time()
108
109for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5):
110    if myid == 0 :
111        domain.write_time()
112
113    if triangle_id is not None :
114        print domain.timestepping_statistics(track_speeds=False,
115                                             triangle_id=triangle_id)
116
117print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.