Changeset 4797


Ignore:
Timestamp:
Nov 7, 2007, 2:19:42 PM (17 years ago)
Author:
nick
Message:

still setting up validation

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_validation/conical_island/run_circular.py

    r4777 r4797  
    1010from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
    1111from anuga.shallow_water import Dirichlet_boundary, Time_boundary
     12from anuga.shallow_water.data_manager import get_maximum_inundation_data
    1213from anuga.abstract_2d_finite_volumes.util import file_function
    1314from anuga.pmesh.mesh_interface import create_mesh_from_regions
     15from anuga.utilities.polygon import read_polygon, plot_polygons
    1416#from cmath import cos,sin,cosh,pi
    1517from math import cos,pi,sin,tan,sqrt
     
    2224#-------------------------
    2325
     26def get_xy(x=0,y=0,r=1,angle=1):
     27    x1 = x + cos(angle)*r
     28    y1 = y + sin(-angle)*r
     29#    print "get",x,y,r,x1,y1,"angle",angle, cos(angle),sin(angle)
     30    return x1,y1,
    2431
    2532def prepare_timeboundary(filename):
     
    5057    for i, line in enumerate(lines):
    5158        fields = line.split()
    52         print 'fields',i,line
     59#        print 'fields',i,line
    5360
    5461        T[i] = float(fields[0])
     
    8289
    8390
    84 
    85 poly_all= [[0,0],[0,25],[30,25],[30,0]]
    86 poly=[[9,10],[9,17.6],[16.6,17.6],[16.6,10]]
    87 internal_poly=[[poly,1]]
     91angles = (0.0, 22.5, 45.0, 67.5, 75.0, 80.0, 85.0, 87.5, 90.0, 92.5,
     92          95.0, 100.0, 105.0, 112.5, 135.0, 157.5, 180.0, 202.5, 225.0,
     93          247.5, 270.0, 292.5, 315.0, 337.5)
     94angles1 = (0.0, 22.5, 45.0, 67.5, 90.0, 112.5, 135.0, 157.5, 180.0, 202.5, 225.0,
     95          247.5, 270.0, 292.5, 315.0, 337.5)
     96
     97x = 12.96
     98y = 13.80
     99r1 = 1.5
     100r2 = 2.5
     101d=0.75
     102res=0.0005
     103
     104poly = []
     105internal_poly=[]
     106for i,angle in enumerate(angles):
     107    angle = ((angle-180)/180.0)*pi
     108#    p1 = get_xy(x,y,r1,angle-0.01745*d)
     109#    p2 = get_xy(x,y,r2,angle-0.01745*d)
     110#    p3 = get_xy(x,y,r2,angle+0.01745*d)
     111#    p4 = get_xy(x,y,r1,angle+0.01745*d)
     112#    poly[i] = [[p1[0],p1[1]],[p2[0],p2[1]],[p3[0],p3[1]],[p4[0],p4[1]]]
     113
     114    p = get_xy(x,y,4,angle)
     115    poly.insert(0,[p[0],p[1]])
     116internal_poly.insert(0,(poly,1))
     117
     118poly1 = []
     119for i,angle in enumerate(angles):
     120    angle = ((angle-180)/180.0)*pi
     121    p = get_xy(x,y,2.,angle)
     122    poly1.insert(0,[p[0],p[1]])
     123   
     124internal_poly.insert(1,(poly1,.1))
     125
     126poly2 = []
     127for i,angle in enumerate(angles):
     128    angle = ((angle-180)/180.0)*pi
     129    p = get_xy(x,y,1.,angle)
     130    poly2.insert(0,[p[0],p[1]])
     131   
     132internal_poly.insert(2,(poly2,1))
     133
     134print internal_poly
     135
     136plot_polygons(internal_poly,'poly.png',)
     137
     138poly_all= [[0,7],[0,25],[30,25],[30,7]]
     139#poly1=[[9,10],[9,17.6],[16.6,17.6],[16.6,10]]
     140#internal_poly.insert(0,(poly1,0.025))
    88141create_mesh_from_regions(poly_all,
    89                          boundary_tags={'wall': [0,1,2],'wave': [3]},
    90                          maximum_triangle_area=10,
     142                         boundary_tags={'wall': [0,2],'wave': [3], 'buffer':[1]},
     143                         maximum_triangle_area=1,
    91144                         interior_regions=internal_poly,
    92145                         filename='temp.msh',
     
    101154#-------------------------
    102155domain.set_quantity('friction', 0.0)
    103 water_height = 0.32
     156water_height = 0.322
    104157domain.set_quantity('stage', water_height)
    105158
     
    121174                z=0.625
    122175
    123             print 'x,y,f',i,j,a,b,z
     176#            print 'x,y,f',i,j,a,b,z
    124177            file.write("%s, %s, %s\n" %(i, j, z))
    125178    file.close()
     
    133186# Set simulation parameters
    134187#-------------------------
    135 domain.set_name('test')  # Name of output sww file
     188domain.set_name('test1')  # Name of output sww file
    136189domain.set_default_order(1)               # Apply second order scheme
    137190domain.set_all_limiters(0.9)              # Max second order scheme (old lim)
     
    160213
    161214# Create and assign boundary objects
    162 boundary_filename="ts2a_edit_input_32"
     215boundary_filename="ts2cnew1_input"
    163216prepare_timeboundary(boundary_filename+'.txt')
    164217
     
    174227Br = Reflective_boundary(domain)
    175228Bd = Dirichlet_boundary([water_height,0,0])
    176 domain.set_boundary({'wave': Bts, 'wall': Bd})
     229domain.set_boundary({'wave': Bts, 'wall': Br, 'buffer':Bd})
    177230
    178231#-------------------------
     
    182235t0 = time.time()
    183236
    184 for t in domain.evolve(yieldstep = 5, finaltime = 75):
     237for t in domain.evolve(yieldstep = 1, finaltime = 32):
    185238    domain.write_time()
    186239#    domain.write_time(track_speeds=False)
    187240
     241for t in domain.evolve(yieldstep = 0.1, finaltime = 36
     242                       ,skip_initial_step = True):
     243    domain.write_time()
     244   
     245for t in domain.evolve(yieldstep = 1, finaltime = 75
     246                       ,skip_initial_step = True):
     247    domain.write_time()
     248
    188249print 'That took %.2f seconds' %(time.time()-t0)
     250
     251
     252
     253#angles = (0.0, 22.5, 45.0, 67.5, 75.0, 80.0, 85.0, 87.5, 90.0, 92.5,
     254#          95.0, 100.0, 105.0, 112.5, 135.0, 157.5, 180.0, 202.5, 225.0,
     255#          247.5, 270.0, 292.5, 315.0, 337.5)
     256#x = 12.96
     257#y = 13.80
     258#r1 = 1.1
     259#r2 = 3.6
     260#d=2
     261#poly = {}
     262for i,angle in enumerate(angles):
     263#    angle = ((angle-180)/180.0)*pi
     264##    angle = angle1*3.14157
     265#    p1 = get_xy(x,y,r1,angle-0.01745*d)
     266#    p2 = get_xy(x,y,r2,angle-0.01745*d)
     267#    p3 = get_xy(x,y,r2,angle+0.01745*d)
     268#    p4 = get_xy(x,y,r1,angle+0.01745*d)
     269#    poly[i] = [[p1[0],p1[1]],[p2[0],p2[1]],[p3[0],p3[1]],[p4[0],p4[1]]]
     270##    print i,poly[i]
     271#   
     272   
     273    run_up, x_y = get_maximum_inundation_data(filename='test.sww',polygon=poly[i], verbose=False)
     274   
     275    print 'maximum_inundation_data',angle, run_up, x_y
     276   
     277
     278
     279
     280
     281
Note: See TracChangeset for help on using the changeset viewer.