# Changeset 1393

Ignore:
Timestamp:
May 15, 2005, 10:25:56 PM (19 years ago)
Message:

Location:
inundation/ga/storm_surge
Files:
10 edited

Unmodified
Removed
• ## inundation/ga/storm_surge/Merimbula/run_merimbula_lake.py

 r1363 domain.store = True     #Store for visualisation purposes domain.format = 'sww'   #Native netcdf visualisation format domain.filename = 'Merimbula_2003_4days_dry' from normalDate import ND domain.filename = 'Merimbula_2003_4days_dry_%s'%ND() print domain.filename #----------------------

 r1387 def rectangular_with_ghosts(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)): """Setup a rectangular grid of triangles with m+1 by n+1 grid points len2: y direction (bottom to top) Also returns a list of Return to lists: points and elements suitable for creating a Mesh or FVMesh object, e.g. Mesh(points, elements) from config import epsilon #E = m*n*2        #Number of triangular elements #P = (m+1)*(n+1)  #Number of initial vertices from Numeric import zeros, Float, Int delta1 = float(len1)/m delta2 = float(len2)/n #Dictionary of vertex objects vertices = {} points = [] #Calculate number of points Np = (m+1)*(n+1) class VIndex: def __init__(self, n,m): self.n = n self.m = m def __call__(self, i,j): return j+i*(self.n+1) class EIndex: def __init__(self, n,m): self.n = n self.m = m def __call__(self, i,j): return 2*(j+i*self.n) I = VIndex(n,m) E = EIndex(n,m) points = zeros( (Np,2), Float) for i in range(m+1): for j in range(n+1): vertices[i,j] = len(points) points.append([i*delta1 + origin[0], j*delta2 + origin[1]]) points[I(i,j),:] = [i*delta1 + origin[0], j*delta2 + origin[1]] #Construct 2 triangles per rectangular element and assign tags to boundary #Calculate number of triangles Nt = 2*m*n #Construct 2 triangles per rectangular element and assign tags to boundary elements = [] elements = zeros( (Nt,3), Int) boundary = {} ghosts = {} nt = -1 for i in range(m): for j in range(n): v1 = vertices[i,j+1] v2 = vertices[i,j] v3 = vertices[i+1,j+1] v4 = vertices[i+1,j] #Update boundary dictionary and create elements i1 = I(i,j+1) i2 = I(i,j) i3 = I(i+1,j+1) i4 = I(i+1,j) #Lower Element nt = E(i,j) if i == m-1: boundary[(len(elements), 2)] = 'right' ghosts[nt] = E(1,j) if i == 0: ghosts[nt] = E(m-2,j) if i == m-1: boundary[nt, 2] = 'right' if j == 0: boundary[(len(elements), 1)] = 'bottom' elements.append([v4,v3,v2]) #Lower element boundary[nt, 1] = 'bottom' elements[nt,:] = [i4,i3,i2] #Upper Element nt = E(i,j)+1 if i == m-1: ghosts[nt] = E(1,j)+1 if i == 0: ghosts[nt] = E(m-2,j)+1 if i == 0: boundary[(len(elements), 2)] = 'left' boundary[nt, 2] = 'left' if j == n-1: boundary[(len(elements), 1)] = 'top' elements.append([v1,v2,v3]) #Upper element boundary[nt, 1] = 'top' elements[nt,:] = [i1,i2,i3] ghosts = {} i=0 for j in range(n): v1 = vertices[i,j+1] v2 = vertices[i,j] v3 = vertices[i+1,j+1] v4 = vertices[i+1,j] ghosts.append(elements.index([ return points, elements, boundary, ghosts return points, elements, boundary

 r1387 from mesh_factory import rectangular points, vertices, boundary = rectangular(30, 30) N = 4 M = 4 points, vertices, boundary, ghosts = rectangular_with_ghosts(N, M) print ghosts #Create advection domain with direction (1,-1) #Check that the boundary value gets propagated to all elements for t in domain.evolve(yieldstep = 0.01, finaltime = 1.5): for t in domain.evolve(yieldstep = 0.01, finaltime = 0.01): domain.write_time()
• ## inundation/ga/storm_surge/pyvolution/data_manager.py

 r1363 from Scientific.IO.NetCDF import NetCDFFile import types from time import sleep from time import sleep from os import stat #Get NetCDF retries = 0 file_open = False while not file_open and retries < 10: try: retries = 0 file_open = False while not file_open and retries < 10: try: fid = NetCDFFile(self.filename, 'a') #Open existing file except IOError: #This could happen if someone was reading the file. #In that case, wait a while and try again msg = 'Warning (store_timestep): File %s could not be opened'\ %self.filename msg += ' - trying again' print msg retries += 1 sleep(1) else: file_open = True if not file_open: msg = 'File %s could not be opened for append' %self.filename raise msg except IOError: #This could happen if someone was reading the file. #In that case, wait a while and try again msg = 'Warning (store_timestep): File %s could not be opened'\ %self.filename msg += ' - trying again' print msg retries += 1 sleep(1) else: file_open = True if not file_open: msg = 'File %s could not be opened for append' %self.filename raise msg q = (1-ratio)*Q[index]+ ratio*Q[index+1] else: q = Q[index] q = Q[index] #Return vector of interpolated values return q
• ## inundation/ga/storm_surge/pyvolution/netherlands.py

 r1392 #Wall #if x[i] > 0.995: #    z[i] = -x[i]/20+0.3 if x[i] > 0.995: z[i] = -x[i]/20+0.3 return z/2 #Set bed-slope and friction inflow_stage = 0.2 manning = 0.00 inflow_stage = 0.1 manning = 0.03 Z = Weir(inflow_stage) t0 = time.time() for t in domain.evolve(yieldstep = 1.0, finaltime = 100.0): from visual import * for t in domain.evolve(yieldstep = 0.1, finaltime = 15.0): domain.write_time() #domain.visualiser.scale_z = 1.0 domain.visualiser.update_quantity_color('xmomentum',scale_z = 4.0) #domain.visualiser.update_quantity_color('xmomentum',scale_z = 4.0) #integral_label.text='Integral=%10.5e'%domain.quantities['stage'].get_integral()
• ## inundation/ga/storm_surge/pyvolution/quantity.py

 r1290 # Create a V like [[0 1 2], [3 4 5]....[3*m-2 3*m-1 3*m]] # These vert_id's will relate to the verts created bellow # These vert_id's will relate to the verts created below m = len(self.domain)  #Number of volumes M = 3*m        #Total number of unique vertices
• ## inundation/ga/storm_surge/pyvolution/test_quantity.py

 r1018 quantity = Quantity(self.mesh4) #Try constants first const = 5 #Try constants first const = 5 quantity.set_values(const, location = 'vertices') #print 'Q', quantity.get_integral() assert allclose(quantity.get_integral(), self.mesh4.get_area() * const) #Try with a linear function #print 'Q', quantity.get_integral() assert allclose(quantity.get_integral(), self.mesh4.get_area() * const) #Try with a linear function def f(x, y): return x+y ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2 assert allclose (quantity.get_integral(), ref_integral) ref_integral = (4.0/3 + 8.0/3 + 10.0/3 + 10.0/3) * 2 assert allclose (quantity.get_integral(), ref_integral) runner = unittest.TextTestRunner() runner.run(suite)
• ## inundation/ga/storm_surge/zeus/anuga-workspace.zwi

 r1387 Debug steve pyvolution Merimbula.zpi parallel.zpi
• ## inundation/ga/storm_surge/zeus/pyvolution.zpi

 r1271 ..\pyvolution\netherlands.py..\pyvolution\netherlands-chris.py..\pyvolution\normalDate.py..\pyvolution\pmesh2domain.py..\pyvolution\pressure_force.py
Note: See TracChangeset for help on using the changeset viewer.