Changeset 5411


Ignore:
Timestamp:
Jun 23, 2008, 1:01:17 PM (17 years ago)
Author:
Leharne
Message:

Updates to okushiri convergence study

Location:
anuga_work/development/convergence_okushiri_2008
Files:
2 added
3 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/development/convergence_okushiri_2008/create_okushiri_truescale.py

    r5343 r5411  
    66
    77from anuga.pmesh.mesh import *
    8 from anuga.pmesh.mesh_interface import create_mesh_from_regions
    98from anuga.coordinate_transforms.geo_reference import Geo_reference
    109from anuga.geospatial_data import Geospatial_data
    … …  
    112111
    113112
    114     #--------------------------------------------------------------------------
    115     # Create the triangular mesh based on overall clipping polygon with a
    116     # tagged
    117     # boundary and interior regions defined in project_truescale.py along with
    118     # resolutions (maximal area of per triangle) for each polygon
    119     #--------------------------------------------------------------------------
    120 
    121 
    122     base_resolution = 1 # Use this to coarsen or refine entire mesh.
    123 
    124     # Basic geometry and bounding polygon
    125     xleft   = 0
    126     xright  = 2179.2
    127     ybottom = 0
    128     ytop    = 1360.8
    129 
    130     point_sw = [xleft, ybottom]
    131     point_se = [xright, ybottom]
    132     point_nw = [xleft, ytop]   
    133     point_ne = [xright, ytop]
    134 
    135     bounding_polygon = [point_se,
    136                         point_ne,
    137                         point_nw,
    138                         point_sw]
    139          
    140     # Localised refined area for gulleys
    141     xl = 1920
    142     xr = 2120
    143     yb = 640
    144     yt = 920
    145     p0 = [xl, yb]
    146     p1 = [xr, yb]
    147     p2 = [xr, yt]
    148     p3 = [xl, yt]
    149113   
    150     gulleys = [p0, p1, p2, p3]
    151    
    152 
    153     # Island area and drawdown region (original)
    154     #island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5]
    155     #island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3]
    156     #island_2 = [xleft + (xright-xleft)/2+0.3, ybottom + 2*(ytop-ybottom)/3-0.3]
    157     #island_3 = [xleft + (xright-xleft)/2+0.3, ybottom + (ytop-ybottom)/3+0.3]
    158     #island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3]
    159     #island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.2]
    160     #island_6 = [xl-.01, yb]  #OK
    161     #island_7 = [xl-.01, yt]  #OK     
    162 
    163 
    164     # Island area and drawdown region
    165     island_0 = [xleft + 2*(xright-xleft)/3+480, ytop-200]
    166     island_1 = [xleft + 2*(xright-xleft)/3+200, ybottom + 2*(ytop-ybottom)/3]
    167     island_2 = [xleft + (xright-xleft)/2+160, ybottom + 2*(ytop-ybottom)/3-120]
    168     island_3 = [xleft + (xright-xleft)/2+160, ybottom + (ytop-ybottom)/3+120]
    169     island_4 = [xleft + 2*(xright-xleft)/3+160, ybottom + (ytop-ybottom)/3-120]
    170     island_5 = [xleft + 2*(xright-xleft)/3+480, ybottom+320]
    171     island_6 = [xl-4, yb]  # Keep right edge just off the gulleys
    172     island_7 = [xl-4, yt]
    173  
    174     island = [island_0, island_1, island_2,
    175               island_3, island_4, island_5,
    176               island_6, island_7]
    177 
    178 
    179     # Region spanning half right hand side of domain just inside boundary (org)
    180     #rhs_nw = [xleft + (xright-xleft)/3+1, ytop-0.02]
    181     #rhs_sw = [xleft + (xright-xleft)/3+1, ybottom+0.02]
    182     #rhs_se = [xright-0.02, ybottom+0.02]
    183     #rhs_ne = [xright-0.02, ytop-0.02]
    184 
    185     # Region spanning half right hand side of domain just inside boundary
    186     rhs_nw = [xleft + (xright-xleft)/3+400, ytop-560]
    187     rhs_sw = [xleft + (xright-xleft)/3+400, ybottom+200]
    188     rhs_se = [xright-4, ybottom+80]
    189     rhs_ne = [xright-4, ytop-80]       
    190 
    191     rhs_region = [rhs_nw, rhs_ne, rhs_se, rhs_sw]
    192 
    193      
    194     # Interior regions and creation of mesh
    195     interior_regions = [[rhs_region, 80],
    196                         [island, 32*base_resolution],
    197                         [gulleys, 3.2*base_resolution]]   
    198 
    199     meshname = project_truescale.mesh_filename + '.msh'
    200     m = create_mesh_from_regions(bounding_polygon,
    201                                  boundary_tags={'wall': [0, 1, 3],
    202                                                 'wave': [2]},     
    203                                  maximum_triangle_area=16000*base_resolution,
    204                                  interior_regions=interior_regions,
    205                                  filename=project_truescale.mesh_filename,
    206                                  verbose=True)
    207 

  • anuga_work/development/convergence_okushiri_2008/project_truescale.py

    r5396 r5411  
    1 """Common filenames for truescale Okushiri Island validation
     1"""Common filenames for truescale Okushiri Island convergence study
    22Formats are given as ANUGA native netCDF where applicable.
    33
    … …  
    1111from anuga.utilities.system_tools import get_user_name, get_host_name
    1212
    13 codename = 'project_grad.py' # FIXME can be obtained automatically as __file__
    14 
    15 home = join(getenv('INUNDATIONHOME'), 'data', 'anuga_validation',
    16             'convergence_okushiri_2008')# Location of Data   
     13home = join(getenv('INUNDATIONHOME'),'data', 'anuga_validation',
     14            'convergence_okushiri_2008') # Location of Data
    1715user = get_user_name()
    1816host = get_host_name()
    … …  
    2018umask(002)
    2119
    22 #time stuff
    23 time = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
    24 run_time = time+'_run'
    25 
    26 #Set anuga directories
    27 anuga_dir = join(home, 'anuga')
    28 
    29 dir_comment='_'+setup+'_'+str(tide)+'_'+\
    30              str(user)+'_'+boundary_event+'_'+which_data
    31 
    32 mesh_dir = join(anuga_dir, 'meshes')
    33 mesh_name = join(mesh_dir, 'okushiri_truescale')
    34 
    35 polygons_dir = join(anuga_dir, 'polygons') # Created with ArcGIS (csv files)
    36 tide_dir = join(anuga_dir, 'tide_data')
    37 
    38 #output locations
    39 output_dir = join(anuga_dir, 'outputs')+sep
    40 output_run_time_dir = output_dir +run_time+dir_comment+sep
    41 output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
    42 
    43 #gauges
    44 gauges_dir = join(anuga_dir,'gauges')
    45 gauge_name = 'gauge_location_onslow.csv'
    46 gauges_dir_name = gauges_dir + gauge_name
     20#-------------------
     21# Input file names
     22#-------------------
    4723
    4824# Given boundary wave
    … …  
    5531bathymetry_filename = 'okushiri_truescale_bathymetry.pts'
    5632
    57 # Triangular mesh
    58 mesh_filename = 'okushiri_truescale.msh'
     33
     34#------------------------------------
     35# Output file names and directories
     36#------------------------------------
    5937
    6038# Model output
    6139output_filename = 'okushiri_truescale.sww'
    6240
     41# Time stuff
     42time = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
     43run_time = time+'_run'
     44
     45# Run parameters
     46finaltime=450
     47setup='original'
     48
     49if setup =='original':
     50    print 'original resolution'
     51    base_resolution=1
     52    yieldstep=1
     53if setup =='double':
     54    print 'double original resolution'
     55    base_resolution=0.5
     56    yieldstep=1
     57if setup =='half':
     58    print 'half original resolution'
     59    base_resolution=2
     60    yieldstep=1
     61if setup =='no polygons':
     62    print 'half original resolution'
     63    base_resolution=2
     64    yieldstep=1
     65 
     66   
     67# Set anuga directories
     68anuga_dir = join(home,'anuga')+sep
     69
     70dir_comment='_'+setup+'_'+str(user)
     71
     72mesh_dir = join(anuga_dir, 'meshes')+sep
     73mesh_name = join(mesh_dir, 'okushiri_truescale')
     74
     75polygons_dir = join(anuga_dir, 'polygons')+sep # Created with ArcGIS (csv files)
     76
     77# Output locations
     78output_dir = join(anuga_dir, 'outputs')+sep
     79output_run_time_dir = output_dir+run_time+dir_comment+sep
     80output_run_time_dir_name = output_run_time_dir + output_filename #Used by post processing
     81
     82# Gauges
     83gauges_dir = join(anuga_dir,'gauges')+sep
     84gauge_name = 'gauge_location_okushiri.csv'
     85gauges_dir_name = gauges_dir+gauge_name
     86
     87# Vertex coordinates
     88vertex_filename = output_run_time_dir+setup+'vertex_coordinates.txt'
     89
     90#------------------------------
     91# Polygon definitions
     92#------------------------------
     93
     94poly_all = read_polygon(polygons_dir+'bounding_polygon.csv')
     95res_poly_all = 16000*base_resolution
     96
     97#print 'Area of bounding polygon', polygon_area(poly_all)/1000000.0
     98
     99poly_gulleys = read_polygon(polygons_dir+'gulleys_polygon.csv')
     100res_gulleys = 3.2*base_resolution
     101
     102poly_island = read_polygon(polygons_dir+'island_polygon.csv')
     103res_island = 32*base_resolution
     104
     105poly_rhs = read_polygon(polygons_dir+'rhs_polygon.csv')
     106res_rhs = 80*base_resolution
     107
     108interior_regions = [[poly_gulleys,res_gulleys],[poly_island,res_island],
     109                    [poly_rhs,res_rhs]]
     110                   
     111trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
     112
     113print 'min number triangles', trigs_min
    63114
    64115
     116
     117
     118
     119

  • anuga_work/development/convergence_okushiri_2008/run_okushiri_truescale.py

    r5396 r5411  
    1 """Validation of the AnuGA implementation of the shallow water wave equation.
     1"""Anuga convergence study using a true-scale version of the Okushiri
     2Island tsunami wavetank experiment
    23
    3 This script sets up a modified version of the Okushiri Island benchmark
     4This script runs a modified version of the Okushiri Island benchmark
    45
    56as published at
    … …  
    1213
    1314This version up-scales the original 1:400 scale wave-tank experiment, to
    14 "true-scale".
     15"true-scale" with mesh defined using interior polygons.
    1516
    16 The original validation data was downloaded and made available in this directory
    17 for convenience but the original data is available at
     17The original validation data is available at
    1818http://www.cee.cornell.edu/longwave/index.cfm?page=benchmark&problem=2
    1919where a detailed description of the problem is also available.
    2020
    21 
    22 Run create_okushiri_truescale.py to process the boundary condition and build a the
    23 mesh before running this script.
     21Run create_okushiri_truescale.py to convert bathymetry and input boundary
     22condition into NetCDF format.
    2423
    2524"""
     25#----------------------------------------
     26# Import necessary modules
     27#----------------------------------------
    2628
    27 # Module imports
     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
    2838from anuga.shallow_water import Domain
    2939from anuga.shallow_water import Reflective_boundary
    3040from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
    3141from anuga.abstract_2d_finite_volumes.util import file_function
    32 
    33 import project_truescale
     42from anuga.shallow_water.data_manager import copy_code_files, start_screen_catcher
     43from anuga.pmesh.mesh_interface import create_mesh_from_regions
    3444
    3545
    36 #-------------------------
     46import project_truescale # Definition of filenames and interior polygons
     47
     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#------------------------------
    3768# Create Domain from mesh
    38 #-------------------------
    39 domain = Domain(project_truescale.mesh_filename, use_cache=True, verbose=True)
     69#------------------------------
     70domain = Domain(project_truescale.mesh_name+'.msh', use_cache=True, verbose=True)
    4071print domain.statistics()
    4172
    … …  
    4374
    4475# Write vertex coordinates to file
    45 filename='okushiri_mesh_vertex_coordinates_truescale.txt'
     76filename=project.vertex_filename
    4677fid=open(filename,'w')
    4778fid.write('x (m), y (m)\n')
    … …  
    6899#-------------------------
    69100domain.set_name(project_truescale.output_filename)  # Name of output sww file
    70 domain.set_datadir(project_truescale.output_dir)
     101domain.set_datadir(project_truescale.output_run_time_dir) # Name of output directory
    71102domain.set_default_order(2)               # Apply second order scheme
    72103domain.set_all_limiters(0.9)              # Max second order scheme (old lim)
    73 domain.set_minimum_storable_height(0.4)   # Don't store h < 0.4m
     104domain.set_minimum_storable_height(0.1)   # Don't store h < 0.1m
    74105domain.tight_slope_limiters = 1
    75106domain.beta_h = 0.0
    76107
    77 #Timings on AMD64-242 (beta_h=0)
    78 #  tight_slope_limiters = 0:
    79 #    3035s - 3110s
    80 #  tight_slope_limiters = 1:
    81 #    3000s - 3008s
    82 #
    83 # beta_h>0: In the order of 3200s
    84108
    85109#-------------------------
    … …  
    109133t0 = time.time()
    110134
    111 for t in domain.evolve(yieldstep = 1, finaltime = 450):
     135for t in domain.evolve(yieldstep=project_truescale.yieldstep, finaltime = 450):
    112136    print domain.timestepping_statistics(track_speeds=False,
    113137                                         triangle_id=triangle_id)
Note: See TracChangeset for help on using the changeset viewer.