Changeset 2312 for development/cairns_2006

Timestamp:

Feb 1, 2006, 2:03:12 PM (19 years ago)

Author:

jakeman

Message:

create_mesh generates an irregular grid with five regions and run_cairns accepts large bathymetry
files for simulation

Location:

development/cairns_2006

Files:

• create_mesh.py (modified) (8 diffs)
• run_cairns.py (modified) (8 diffs)

development/cairns_2006/create_mesh.py

@@ -22 +22 @@
     coord_sw = convert_latlon_to_xycoords(-157,-24,latlong_origin)
     coord_nw = convert_latlon_to_xycoords(-157,-8,latlong_origin)
-    
-    """
-    xleft = 0
-    ybottom = 0
-    xright = 100
-    ytop = 100
-
-    point_sw = [xleft, ybottom]
-    point_se = [xright, ybottom]
-    point_nw = [xleft, ytop]    
-    point_ne = [xright, ytop]
-
-    """
-    
-    """
-    #Outline
-    point_sw = [coord_sw[0], coord_sw[1]]
-    point_se = [coord_se[0], coord_se[1]]
-    point_nw = [coord_nw[0], coord_nw[1]]    
-    point_ne = [coord_ne[0], coord_ne[1]]
-    """
 
     # Sets origin which is needed by an instance of a mesh object
     geo = Geo_reference(xllcorner = coord_se[0],
                         yllcorner = coord_se[1])
-    
-    #    geo = Geo_reference(xllcorner = xleft,
-    #                        yllcorner = ybottom)
-
-    """
-    coord_nw = [3.6e6,1e6]
-    coord_ne = [0.0,1e6]
-    coord_sw = [3.6e6,0.0]
-    coord_se = [0.0,0.0]
-
-    geo = Geo_reference(xllcorner = 0.0,
-                        yllcorner = 0.0)
-    """
     
     print "***********************"
@@ -69 +35 @@
     #Boundary
     dict = {}
-
-    """
-    dict['points'] = [point_se,
-                      point_ne,
-                      point_nw,
-                      point_sw]
-    
-    """
     
     dict['points'] = [coord_se,   
@@ -82 +40 @@
                       coord_nw,
                       coord_sw]
-    
-#    print dict['points']
     
     dict['segments'] = [[0,1], [1,2], [2,3], [3,0]]  #The outer border
@@ -91 +47 @@
                             'border',
                             'border']
-
         
     m.addVertsSegs(dict)
     
+
     dict = {}
     # rupture zone
     # Tonga-Kermadec Subduction zone
     # located between 180E and -173W and 15S and 35S
-    p0 = convert_latlon_to_xycoords(-173,-20,latlong_origin)
-    p1 = convert_latlon_to_xycoords(-173,-15,latlong_origin)
-    p2 = convert_latlon_to_xycoords(-171,-15,latlong_origin)
-    p3 = convert_latlon_to_xycoords(-171,-20,latlong_origin)
+    r0 = convert_latlon_to_xycoords(-173,-20,latlong_origin)
+    r1 = convert_latlon_to_xycoords(-173,-15,latlong_origin)
+    r2 = convert_latlon_to_xycoords(-171,-15,latlong_origin)
+    r3 = convert_latlon_to_xycoords(-171,-20,latlong_origin)
 
-    """
-    rupture_south = degminsec2decimal_degrees(-15,0,0)
-    rupture_north = degminsec2decimal_degrees(-17,0,0)
-    rupture_west = degminsec2decimal_degrees(-165,0,0)
-    rupture_east = degminsec2decimal_degrees(-166,0,0)
-
-    p0 = [rupture_south, rupture_west]
-    p1 = [rupture_south, rupture_east]
-    p2 = [rupture_north, rupture_east]
-    p3 = [rupture_north, rupture_west]
-    """
-    
-
-    dict['points']  = [p0, p1, p2, p3]
+    dict['points']  = [r0, r1, r2, r3]
     dict['segments'] = [[0,1], [1,2], [2,3], [3,0]]
     dict['segment_tags'] = ['', '', '', '']
@@ -124 +67 @@
     #Interior regions
     # Cairns lat and long (16.575S, 145.45E)
-    c0 = convert_latlon_to_xycoords(145,-24,latlong_origin)
-    c1 = convert_latlon_to_xycoords(145,-8,latlong_origin)
-    c2 = convert_latlon_to_xycoords(146,-8,latlong_origin)
-    c3 = convert_latlon_to_xycoords(153,-24,latlong_origin)
+    c0 = convert_latlon_to_xycoords(149,-22,latlong_origin)
+    c1 = convert_latlon_to_xycoords(146,-19,latlong_origin)
+    c2 = convert_latlon_to_xycoords(145,-15,latlong_origin)
+    c3 = convert_latlon_to_xycoords(145,-13,latlong_origin)
+    c4 = convert_latlon_to_xycoords(146,-13,latlong_origin)   
+    c5 = convert_latlon_to_xycoords(147,-17,latlong_origin)
+    c6 = convert_latlon_to_xycoords(153,-21,latlong_origin)
+    c7 = convert_latlon_to_xycoords(153,-22,latlong_origin)
 
-    """
-    coast_south = degminsec2decimal_degrees(-16,0,0)
-    coast_north = degminsec2decimal_degrees(-17,0,0)
-    coast_west = degminsec2decimal_degrees(145,0,0)
-    coast_east = degminsec2decimal_degrees(145,75,0)
-    
-    d0 = [coast_south, coast_west]
-    d1 = [coast_south, coast_east]
-    d2 = [coast_north, coast_east]
-    d3 = [coast_north, coast_west]
-    """
-
-    dict['points']  = [c0, c1, c2, c3]
-    dict['segments'] = [[0,1], [1,2], [2,3], [3,0]]
-    dict['segment_tags'] = ['', '', '', '']
+    dict['points']  = [c0, c1, c2, c3, c4, c5, c6, c7]
+    dict['segments'] = [[0,1], [1,2], [2,3], [3,4], [4,5], [5,6], [6,7], [7,0]]
+    dict['segment_tags'] = ['', '', '', '', '', '', '', '']
     m.addVertsSegs(dict)
 
@@ -179 +114 @@
     m.addVertsSegs(dict)
 
-    """
-    Generates a rectangular region in the mesh.  Could reduce the number of
-    traingles needed for this refined area by creating a polygon that better
-    resmebles the shape of the reef so areas of deep sea floor are not covered
-    by un-necessary small traiangels
 
-
-    #Interior regions
-    reef1_south = degminsec2decimal_degrees(-20,44,0)
-    reef1_north = degminsec2decimal_degrees(-20,42,0)
-    reef1_west = degminsec2decimal_degrees(116,48,0)
-    reef1_east = degminsec2decimal_degrees(116,53,30)
+    #Exclude PNG because land mass has no effect on the wave that hits Cairns
+    # Papua New Guinea
+    p0 = convert_latlon_to_xycoords(145.5,-12.0,latlong_origin)
+    p1 = convert_latlon_to_xycoords(145.5,-8.0,latlong_origin)
+    p2 = convert_latlon_to_xycoords(154.5,-8.0,latlong_origin)
+    p3 = convert_latlon_to_xycoords(154.5,-12.0,latlong_origin)
     
-    k0 = [reef1_south, reef1_west]
-    k1 = [reef1_south, reef1_east]
-    k2 = [reef1_north, reef1_east]
-    k3 = [reef1_north, reef1_west]    
-    
-    dict['points']  = reef1_polygon = [k0, k1, k2, k3] 
+    dict['points']  = [p0, p1, p2, p3]
     dict['segments'] = [[0,1], [1,2], [2,3], [3,0]]
     dict['segment_tags'] = ['', '', '', '']
     m.addVertsSegs(dict)
-    """
-    
-    base_resolution = 5e11
+
 
     """
@@ -212 +135 @@
     m^2
     """
-    print coord_nw[0]-0.1
-    print coord_ne[1]-0.1
+
+    # base resolution = 5e11 gives good size test grid
+    #base_resolution = 5e11
+    base_resolution = 5e9
+    
     ocean = m.addRegionEN(coord_nw[0]-0.1, coord_nw[1]-0.1)
     ocean.setMaxArea(base_resolution)
+    
 
-    #    ocean = m.addRegionEN(xright-0.1, ytop-0.1)
-    #    ocean.setMaxArea(base_resolution)
-
-    rupture = m.addRegionEN(p0[0]+0.1, p0[1]+0.1)
-    rupture.setMaxArea(base_resolution)
+    rupture = m.addRegionEN(r0[0]+0.1, r0[1]+0.1)
+    rupture.setMaxArea(100*base_resolution)
 
     coast = m.addRegionEN(c0[0]+0.1, c0[1]+0.1)
-    coast.setMaxArea(0.001*base_resolution)
+    coast.setMaxArea(0.003*base_resolution)
 
     fiji = m.addRegionEN(f0[0]+0.1, f0[1]+0.1)
@@ -233 +157 @@
 
     melanesia = m.addRegionEN(m0[0]+0.1, m0[1]+0.1)
-    melanesia.setMaxArea(0.001*base_resolution)
+    melanesia.setMaxArea(0.01*base_resolution)
+    
+    png = m.addRegionEN(p0[0]+0.1, p0[1]+0.1)
+    png.setMaxArea(0.01*base_resolution)
     
     #a10000a refers to the max area passed to triangle

development/cairns_2006/run_cairns.py

@@ -10 +10 @@
 import sys, os
 from pyvolution.shallow_water import Domain, Reflective_boundary,\
-     File_boundary, Transmissive_Momentum_Set_Stage_boundary, Transmissive_boundary
+     File_boundary, Transmissive_Momentum_Set_Stage_boundary,\
+     Transmissive_boundary, Dirichlet_boundary
 from pyvolution.mesh_factory import rectangular_cross
 from pyvolution.pmesh2domain import pmesh_to_domain_instance
@@ -19 +20 @@
 from utilities.polygon import read_polygon, Polygon_function
 
+from conversion import convert_latlon_to_xycoords
+
 #-------------------------------
 # Domain
 #-------------------------------
 print 'Creating domain'
+
 """
 M = 36 #float(len1)/m
@@ -36 +40 @@
 
 """
+
 domain = cache(pmesh_to_domain_instance,
                (project.mesh_filename, Domain),
@@ -53 +58 @@
 print 'Initial values'
 
-p0 = read_polygon('cairns_fault.xya')
+latlong_origin = 145, -24
+r0 = convert_latlon_to_xycoords(-173,-20,latlong_origin)
+r1 = convert_latlon_to_xycoords(-173,-15,latlong_origin)
+r2 = convert_latlon_to_xycoords(-171,-15,latlong_origin)
+r3 = convert_latlon_to_xycoords(-171,-20,latlong_origin)
+
+#p0 = read_polygon('cairns_fault.xya')
 #p0 = read_polygon('cairns_fault_degrees.xya')
-
-domain.set_quantity('stage',Polygon_function([(p0,1000.0)]))
+p0 = [r0, r1, r2, r3]
+domain.set_quantity('stage',Polygon_function([(p0,1.0)]))
 
 #print domain.quantities['stage'].vertex_values
@@ -77 +88 @@
 print 'Boundaries'
 
-
 #Bts = Transmissive_Momentum_Set_Stage_boundary(domain, tide_function)
-
+#Bd = Dirichlet_boundary([0,0,0])
 Br = Reflective_boundary(domain)
 domain.set_boundary({'border': Br, 'border': Br, 'border': Br, 'border': Br})
 
 #Bt = Transmissive_boundary(domain)
-#domain.set_boundary({'left': Bt, 'right': Bt, 'top': Bt, 'bottom': Bt})
+#domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
 
 #-------------------------------
@@ -90 +100 @@
 #-------------------------------
 
-domain.visualise = True
-domain.visualise_color_stage = True
+#domain.visualise = True
+#domain.visualise_color_stage = True
 
-#domain.visualise = False
-#domain.visualise_color_stage = False
+#Set domain.visualise = False for large problem sizes
+domain.visualise = False
+domain.visualise_color_stage = False
 
 base = os.path.basename(sys.argv[0])
@@ -107 +118 @@
 t0 = time.time()
 yieldstep = 60.0
-finaltime = 480.0
+finaltime = 21600.0
 
 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
@@ -113 +124 @@
 
 print 'That took %.2f seconds' %(time.time()-t0)
-