Changeset 7579


Ignore:
Timestamp:
Dec 7, 2009, 6:18:23 PM (9 years ago)
Author:
ole
Message:

Cleanup of school exercise. Optimisation to follow

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/classroom/WORKING-RipCurrent.py

    r7578 r7579  
    55# Import necessary modules
    66#------------------------------------------------------------------------------
    7 from anuga.abstract_2d_finite_volumes.mesh_factory import *
    8 from anuga.shallow_water import *
    9 from anuga.shallow_water.shallow_water_domain import *
    10 from math import *
    11 from numpy import *
     7from anuga.interface import create_domain_from_regions
     8from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
     9from anuga.shallow_water.shallow_water_domain import Reflective_boundary
     10from anuga.shallow_water.shallow_water_domain import Time_boundary
     11from pylab import figure, quiver, show, cos, sin, pi
    1212import numpy
    13 from matplotlib import *
    14 from pylab import *
    1513import csv
    16 
    17 #Parameters
    18 
     14import time
     15
     16
     17#------------------------------------------------------------------------------
     18# Parameters
     19#------------------------------------------------------------------------------
    1920filename = "WORKING-RIP-LAB_Expt-Geometry_Triangular_Mesh"
    20 location_of_shore = 140    #the position along the y axis of the shorefront
    21 sandbar = 1.2 #height of sandbar
    22 sealevel = 0 #height of coast above sea level
    23 steepness = 8000 #period of sandbar - larger number gives smoother slope - longer period
     21location_of_shore = 140 # The position along the y axis of the shorefront
     22sandbar = 1.2           # Height of sandbar
     23sealevel = 0            # Height of coast above sea level
     24steepness = 8000        # Period of sandbar -
     25                        # larger number gives smoother slope - longer period
    2426halfchannelwidth = 5
    2527bank_slope = 0.1
     
    2729timestep = 1
    2830
     31
    2932#------------------------------------------------------------------------------
    3033# Setup computational domain
     
    3235length = 120
    3336width = 170
    34 seafloor_resolution = 60.0 # Resolution: Max area of triangles in the mesh for seafloor
     37seafloor_resolution = 60.0 # Resolution: Max area of triangles in the mesh
    3538feature_resolution = 1.0
    3639beach_resolution = 10.0
     
    4245meshname = str(filename)+'.msh'
    4346
    44 #interior regions
     47# Interior regions
    4548feature_regions = [[feature_boundary_polygon, feature_resolution],
    4649                   [beach_interior_polygon, beach_resolution]]
    4750
    48 create_mesh_from_regions(sea_boundary_polygon,
    49                          boundary_tags={'bottom': [0],
    50                                         'right' : [1],
    51                                         'top' :   [2],
    52                                         'left':   [3]},
    53                          maximum_triangle_area = seafloor_resolution,
    54                          filename = meshname,
    55                          interior_regions = feature_regions,
    56                          use_cache = False,
    57                          verbose = False)
    58 domain = Domain(meshname, use_cache=True, verbose=True)
     51domain = create_domain_from_regions(sea_boundary_polygon,
     52                                    boundary_tags={'bottom': [0],
     53                                                   'right' : [1],
     54                                                   'top' :   [2],
     55                                                   'left':   [3]},
     56                                    maximum_triangle_area = seafloor_resolution,
     57                                    mesh_filename = meshname,
     58                                    interior_regions = feature_regions,
     59                                    use_cache = True,
     60                                    verbose = True)
     61                                   
    5962domain.set_name(filename)                  # Output name
    6063print domain.statistics()
    6164
     65
    6266#------------------------------------------------------------------------------
    6367# Setup initial conditions
    6468#------------------------------------------------------------------------------
    65 
    6669def topography(x,y):
    6770    """Complex topography defined by a function of vectors x and y."""
    6871
    69     #general slope and buildings
     72    # General slope and buildings
    7073    z=0.05*(y-(location_of_shore))
    7174   
     
    7679    for i in range(N):
    7780        if y[i]>150:
    78             z[i] = (0.1*(y[i]-150)) + 0.05*(y[i] - (location_of_shore))
     81            z[i] = (0.1*(y[i]-150)) + 0.05*(y[i]-(location_of_shore))
    7982
    8083    return z
     
    8588   
    8689    N = len(x)
    87     for i in range(N):
    88         if -1*(bank_slope)*x[i] + 112 < y[i] < -1*(bank_slope)*x[i] + 124 and 0<x[i]<((length/2)-halfchannelwidth):
     90   
     91    # It would be great with a comment about what this does
     92    for i in range(N):
     93        ymin = -1*(bank_slope)*x[i] + 112
     94        ymax = -1*(bank_slope)*x[i] + 124
     95        xmin = 0
     96        xmax = (length/2)-halfchannelwidth
     97        if ymin < y[i] < ymax and xmin < x[i]< xmax:
    8998            z[i] += sandbar*((cos((y[i]-118)/steepness)))
    90     for i in range(N):
    91         if -1*(bank_slope)*(x[i]-(length/2)) + (-1*(bank_slope)*(length/2)+112) < y[i] < -1*(bank_slope)*(x[i]-(length/2)) + (-1*(bank_slope)*(length/2)+124) and ((length/2)+halfchannelwidth)<x[i]<183:
     99   
     100    # It would be great with a comment about what this does       
     101    for i in range(N):
     102        ymin = -1*(bank_slope)*(x[i]-(length/2)) + (-1*(bank_slope)*(length/2)+112)
     103        ymax = -1*(bank_slope)*(x[i]-(length/2)) + (-1*(bank_slope)*(length/2)+124)
     104        xmin = (length/2)+halfchannelwidth
     105        xmax = 183
     106        if ymin < y[i] < ymax and xmin < x[i] < xmax:
    92107            z[i] += sandbar*(cos((y[i]-118)/steepness))
     108           
    93109    return z
    94110
    95 domain.set_quantity('elevation', topography)           # elevation is a function
    96 domain.set_quantity('friction', 0.01)                  # Constant friction
    97 
    98 domain.add_quantity('elevation', topography3)
    99 
    100 domain.set_quantity('stage', 0)   # Dry initial condition
     111domain.set_quantity('elevation', topography)  # Apply base elevation function
     112domain.add_quantity('elevation', topography3) # Add elevation modification
     113domain.set_quantity('friction', 0.01)         # Constant friction
     114domain.set_quantity('stage', 0)               # Constant initial condition at mean sea level
     115
    101116
    102117#------------------------------------------------------------------------------
     
    107122Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
    108123
    109 
    110 amplitude = 0.4 #amplitude of wave (wave height m)
    111 period =  5     #wave period (sec)
    112 
    113124def wave(t):
    114 
    115     A = amplitude # Amplitude [m] (Wave height)
    116     T = period  # Wave period [s]
     125    """Define wave driving the system
     126    """
     127   
     128    A = 0.4  # Amplitude of wave [m] (wave height)
     129    T = 5    # Wave period [s]
    117130
    118131    if t < 30000000000:
     
    126139domain.set_boundary({'left': Br, 'right': Br, 'top': Bo, 'bottom': Bt})
    127140
     141
    128142#------------------------------------------------------------------------------
    129143# Evolve system through time
    130144#------------------------------------------------------------------------------
    131145
    132 list1 = [x for x in csv.reader(open('New_gauges.csv','r'),dialect='excel',delimiter=",") ]
    133 
    134 print list1
    135 
    136 u = numpy.zeros(len(list1), 'd')
    137 v = numpy.zeros(len(list1), 'd')
    138 
    139 
    140 for t in domain.evolve(yieldstep = (timestep), finaltime = (simulation_length)):
     146# Read in gauge locations for interpolation and convert them to floats
     147gauge_file = open('New_gauges.csv', 'r')
     148G = [(float(x[0]), float(x[1])) for x in csv.reader(gauge_file,
     149                                                    dialect='excel',
     150                                                    delimiter=',')]
     151gauges = numpy.array(G)
     152number_of_gauges = len(gauges)
     153print gauges                               
     154
     155# Allocate space for velocity values
     156u = numpy.zeros(number_of_gauges)
     157v = numpy.zeros(number_of_gauges)
     158
     159t0 = time.time()
     160for t in domain.evolve(yieldstep = timestep, finaltime = simulation_length):
    141161    print domain.timestepping_statistics()
    142     S = domain.get_quantity('stage').get_values(interpolation_points=list1)
    143     E = domain.get_quantity('elevation').get_values(interpolation_points=list1)
     162    S = domain.get_quantity('stage').get_values(interpolation_points=gauges)
     163    E = domain.get_quantity('elevation').get_values(interpolation_points=gauges)
    144164    depth = S-E
    145165
    146     uh = domain.get_quantity('xmomentum').get_values(interpolation_points=list1)
    147     vh = domain.get_quantity('ymomentum').get_values(interpolation_points=list1)
     166    uh = domain.get_quantity('xmomentum').get_values(interpolation_points=gauges)
     167    vh = domain.get_quantity('ymomentum').get_values(interpolation_points=gauges)
    148168    u += uh/depth
    149169    v += vh/depth
    150 
    151 n_time_intervals = (simulation_length)/(timestep)
    152 u_average = u / (n_time_intervals)
    153 v_average = v / (n_time_intervals)
    154 
    155 print "there were", n_time_intervals, "time steps"
    156 
    157 print "sum y velocity", v
    158 print "average y velocity", v_average
    159 print "sum x velocity", u
    160 print "average x velocity", u_average
    161 
    162 x_output = file("x_velocity.txt", 'w')
    163 y_output = file("y_velocity.txt", 'w')
    164 
    165 print >> x_output, " "
    166 print >> y_output, " "
    167 
    168 print >> x_output, u_average
    169 print >> y_output, v_average
    170 
    171 X = [x[0] for x in csv.reader(open('New_gauges.csv','r'),dialect='excel',delimiter=",") ]
    172 Y = [x[1] for x in csv.reader(open('New_gauges.csv','r'),dialect='excel',delimiter=",") ]
     170print 'Evolution took %.2f seconds' % (time.time()-t0)
     171
     172   
     173#------------------------------------------------------------------------------
     174# Post processing
     175#------------------------------------------------------------------------------   
     176
     177# Use this
     178
     179#     from anuga.fit_interpolate.interpolate import Interpolation_function
     180
     181#     # Get mesh and quantities from sww file
     182#     X = get_mesh_and_quantities_from_file(filename,
     183#                                           quantities=quantity_names,
     184#                                           verbose=verbose)
     185#     mesh, quantities, time = X
     186
     187#     # Find all intersections and associated triangles.
     188#     segments = mesh.get_intersecting_segments(polyline, verbose=verbose)
     189
     190#     # Get midpoints
     191#     interpolation_points = segment_midpoints(segments)
     192
     193#     # Interpolate
     194#     if verbose:
     195#         log.critical('Interpolating - total number of interpolation points = %d'
     196#                      % len(interpolation_points))
     197
     198#     I = Interpolation_function(time,
     199#                                quantities,
     200#                                quantity_names=quantity_names,
     201#                                vertex_coordinates=mesh.nodes,
     202#                                triangles=mesh.triangles,
     203#                                interpolation_points=interpolation_points,
     204#                                verbose=verbose)
     205
     206#     return segments, I
     207
     208
     209
     210
     211
     212
     213n_time_intervals = simulation_length/timestep
     214u_average = u/n_time_intervals
     215v_average = v/n_time_intervals
     216
     217#print "there were", n_time_intervals, "time steps"
     218
     219#print "sum y velocity", v
     220#print "average y velocity", v_average
     221#print "sum x velocity", u
     222#print "average x velocity", u_average
     223
     224x_output = file('x_velocity.txt', 'w')
     225y_output = file('y_velocity.txt', 'w')
     226
     227#print >> x_output, " "
     228#print >> y_output, " "
     229
     230#print >> x_output, u_average
     231#print >> y_output, v_average
     232
     233
     234X = gauges[:,0]
     235Y = gauges[:,1]
     236                             
    173237U = u_average.tolist()
    174238V = v_average.tolist()
    175239
    176 print "U = ", U
    177 print "U has type", type(U)
    178 
    179 
    180 from matplotlib import *
    181 from pylab import *
     240#print "U = ", U
     241#print "U has type", type(U)
     242
    182243
    183244figure()
Note: See TracChangeset for help on using the changeset viewer.