Changeset 7452
 Timestamp:
 Sep 1, 2009, 10:29:19 AM (11 years ago)
 Location:
 anuga_core/source
 Files:

 8 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/interface.py
r7326 r7452 18 18 19 19 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Time_boundary 20 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Transmissive_boundary 21 20 22 from anuga.abstract_2d_finite_volumes.util import file_function 23 24 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 21 25 22 26 from anuga.shallow_water.data_manager import export_grid, create_sts_boundary … … 26 30 from anuga.utilities.polygon import Polygon_function 27 31 32 from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance 33 34 35 # 36 # rectangular domains 37 # 38 def rectangular_cross_domain(*args, **kwargs): 39 points, vertices, boundary = rectangular_cross(*args, **kwargs) 40 return Domain(points, vertices, boundary) 41 42 # 43 # Create domain from file 44 # 45 def create_domain_from_file(file): 46 return pmesh_to_domain_instance(file,Domain) 28 47 29 48 # 
anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r7352 r7452 557 557 # @note 'polyline' may contain multiple sections allowing complex shapes. 558 558 # @note Assume absolute UTM coordinates. 559 def new_get_energy_through_cross_section(self, polyline,559 def get_energy_through_cross_section(self, polyline, 560 560 kind='total', 561 561 verbose=False): … … 588 588 cross_section = Cross_section(self, polyline, verbose) 589 589 590 return cross_section.get_energy_through_cross_section( )590 return cross_section.get_energy_through_cross_section(kind) 591 591 592 592 … … 598 598 # @note 'polyline' may contain multiple sections allowing complex shapes. 599 599 # @note Assume absolute UTM coordinates. 600 def get_energy_through_cross_section(self, polyline,600 def old_get_energy_through_cross_section(self, polyline, 601 601 kind='total', 602 602 verbose=False): … … 2786 2786 self.midpoints = ensure_geospatial(self.midpoints, self.domain.geo_reference) 2787 2787 2788 ## 2789 # @brief set verbose mode 2790 def set_verbose(self,verbose=True): 2791 """Set verbose mode true or flase 2792 """ 2793 2794 self.verbose=verbose 2788 2795 2789 2796 ## … … 2820 2827 ## 2821 2828 # @brief calculate current energy flow through cross section 2822 def get_energy_through_cross_section(self ):2829 def get_energy_through_cross_section(self, kind='total'): 2823 2830 """Obtain average energy head [m] across specified cross section. 2824 2831 … … 2837 2844 """ 2838 2845 2846 from anuga.config import g, epsilon, velocity_protection as h0 2847 2839 2848 # Get interpolated values 2840 2849 stage = self.domain.get_quantity('stage') … … 2843 2852 ymomentum = self.domain.get_quantity('ymomentum') 2844 2853 2845 w = stage.get_values(interpolation_points= midpoints, use_cache=True)2846 z = elevation.get_values(interpolation_points= midpoints, use_cache=True)2847 uh = xmomentum.get_values(interpolation_points= midpoints,2854 w = stage.get_values(interpolation_points=self.midpoints, use_cache=True) 2855 z = elevation.get_values(interpolation_points=self.midpoints, use_cache=True) 2856 uh = xmomentum.get_values(interpolation_points=self.midpoints, 2848 2857 use_cache=True) 2849 vh = ymomentum.get_values(interpolation_points= midpoints,2858 vh = ymomentum.get_values(interpolation_points=self.midpoints, 2850 2859 use_cache=True) 2851 2860 h = wz # Depth … … 2853 2862 # Compute total length of polyline for use with weighted averages 2854 2863 total_line_length = 0.0 2855 for segment in se gments:2864 for segment in self.segments: 2856 2865 total_line_length += segment.length 2857 2866 … … 2879 2888 2880 2889 # Add to weighted average 2881 weigth = se gments[i].length/total_line_length2890 weigth = self.segments[i].length/total_line_length 2882 2891 average_energy += segment_energy*weigth 2883 2892 
anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py
r7398 r7452 1644 1644 verbose=False) 1645 1645 assert num.allclose(Et, w + 0.5*u*u/g) 1646 1647 1648 def test_cross_section_class(self): 1649 """test_cross_section_class(self): 1650 1651 Test that the total and specific energy through a cross section can be 1652 correctly obtained at runtime from the ANUGA cross section class. 1653 1654 This test creates a flat bed with a known flow through it, creates a cross 1655 section and tests that the correct flow and energies are calculated 1656 1657 The specifics are 1658 e = 1 m 1659 u = 2 m/s 1660 h = 2 m 1661 w = 3 m (width of channel) 1662 1663 q = u*h*w = 12 m^3/s 1664 1665 This run tries it with georeferencing and with elevation = 1 1666 """ 1667 1668 import time, os 1669 from Scientific.IO.NetCDF import NetCDFFile 1670 from mesh_factory import rectangular 1671 1672 # Create basic mesh (20m x 3m) 1673 width = 3 1674 length = 20 1675 t_end = 1 1676 points, vertices, boundary = rectangular(length, width, length, width) 1677 1678 # Create shallow water domain 1679 domain = Domain(points, vertices, boundary, 1680 geo_reference=Geo_reference(56, 308500, 6189000)) 1681 1682 domain.default_order = 2 1683 domain.set_quantities_to_be_stored(None) 1684 1685 e = 1.0 1686 w = 1.0 1687 h = we 1688 u = 2.0 1689 uh = u*h 1690 1691 Br = Reflective_boundary(domain) # Side walls 1692 Bd = Dirichlet_boundary([w, uh, 0]) # 2 m/s across the 3 m inlet: 1693 1694 # Initial conditions 1695 domain.set_quantity('elevation', e) 1696 domain.set_quantity('stage', w) 1697 domain.set_quantity('xmomentum', uh) 1698 domain.set_boundary({'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br}) 1699 1700 # Interpolation points down the middle 1701 I = [[0, width/2.], 1702 [length/2., width/2.], 1703 [length, width/2.]] 1704 interpolation_points = domain.geo_reference.get_absolute(I) 1705 1706 # Shortcuts to quantites 1707 stage = domain.get_quantity('stage') 1708 xmomentum = domain.get_quantity('xmomentum') 1709 ymomentum = domain.get_quantity('ymomentum') 1710 1711 1712 # Create some cross sections 1713 cross_sections = [] 1714 for i in range(5): 1715 x = length/2. + i*0.23674563 # Arbitrary 1716 polyline = [[x, 0], [x, width]] 1717 1718 polyline = domain.geo_reference.get_absolute(polyline) 1719 1720 cross_sections.append(Cross_section(domain,polyline)) 1721 1722 1723 1724 for t in domain.evolve(yieldstep=0.1, finaltime=t_end): 1725 # Check that quantities are they should be in the interior 1726 w_t = stage.get_values(interpolation_points) 1727 uh_t = xmomentum.get_values(interpolation_points) 1728 vh_t = ymomentum.get_values(interpolation_points) 1729 1730 assert num.allclose(w_t, w) 1731 assert num.allclose(uh_t, uh) 1732 assert num.allclose(vh_t, 0.0, atol=1.0e6) 1733 1734 1735 # Check flows and energies through the middle 1736 for cross_section in cross_sections: 1737 1738 Q = cross_section.get_flow_through_cross_section() 1739 1740 assert num.allclose(Q, uh*width) 1741 1742 Es = cross_section.get_energy_through_cross_section(kind='specific') 1743 1744 assert num.allclose(Es, h + 0.5*u*u/g) 1745 1746 Et = cross_section.get_energy_through_cross_section(kind='total') 1747 1748 assert num.allclose(Et, w + 0.5*u*u/g) 1749 1750 1751 1752 1753 1646 1754 1647 1755 def test_another_runup_example(self): 
anuga_core/source/anuga/visualiser/offline.py
r4325 r7452 1 from Numeric import array, Float, ravel, zeros 1 #from Numeric import array, Float, ravel, zeros 2 import numpy as num 2 3 from Scientific.IO.NetCDF import NetCDFFile 3 4 from Tkinter import Button, E, Tk, W, Label, StringVar, Scale, HORIZONTAL … … 80 81 else: 81 82 N_vert = len(fin.variables[quantityName]) 82 x = ravel(array(fin.variables['x'], Float))83 y = ravel(array(fin.variables['y'], Float))83 x = num.ravel(num.array(fin.variables['x'], num.float)) 84 y = num.ravel(num.array(fin.variables['y'], num.float)) 84 85 if dynamic is True: 85 q = array(fin.variables[quantityName][frameNumber], Float)86 else: 87 q = ravel(array(fin.variables[quantityName], Float))86 q = num.array(fin.variables[quantityName][frameNumber], num.float) 87 else: 88 q = num.ravel(num.array(fin.variables[quantityName], num.float)) 88 89 89 90 q *= self.height_zScales[quantityName] … … 123 124 for q in filter(lambda n:n != 'x' and n != 'y' and n != 'z' and n != 'time' and n != 'volumes', fin.variables.keys()): 124 125 if len(fin.variables[q].shape) == 1: # Not a timevarying quantity 125 quantities[q] = ravel(array(fin.variables[q], Float))126 quantities[q] = num.ravel(num.array(fin.variables[q], num.float)) 126 127 else: # Timevarying, get the current timestep data 127 quantities[q] = array(fin.variables[q][self.frameNumber], Float)128 quantities[q] = num.array(fin.variables[q][self.frameNumber], num.float) 128 129 fin.close() 129 130 return quantities 
anuga_core/source/anuga/visualiser/realtime.py
r6113 r7452 1 from Numeric import Float, zeros, shape 1 #from import Float, zeros, shape 2 import numpy as num 2 3 from Tkinter import Button, E, Tk, W 3 4 from threading import Event … … 47 48 48 49 # Also build vert_index  a list of the x & y values of each vertex 49 self.vert_index = zeros((N_vert,2), Float)50 self.vert_index = num.zeros((N_vert,2), num.float) 50 51 for n in range(N_tri): 51 52 self.vtk_cells.InsertNextCell(3) … … 56 57 def update_height_quantity(self, quantityName, dynamic=True): 57 58 N_vert = len(self.source.get_vertex_coordinates()) 58 qty_index = zeros(N_vert, Float)59 qty_index = num.zeros(N_vert, num.float) 59 60 triangles = self.source.get_triangles() 60 61 vertex_values, _ = self.source.get_quantity(quantityName).get_vertex_values(xy=False, smooth=False) 
anuga_core/source/anuga_parallel/test_parallel_sw_flow.py
r7449 r7452 15 15 import os 16 16 import sys 17 import pypar17 #import pypar 18 18 import numpy as num 19 20 from anuga.utilities.numerical_tools import ensure_numeric21 from anuga.utilities.polygon import is_inside_polygon22 19 23 20 from anuga.interface import Domain … … 30 27 31 28 32 from anuga_parallel.interface import distribute, myid, numprocs, send, receive 29 from anuga_parallel.interface import distribute, myid, numprocs, send, receive, barrier 33 30 34 31 # … … 38 35 finaltime = 6.0 39 36 nprocs = 4 40 N = 2 941 M = 2 942 verbose = True37 N = 25 38 M = 25 39 verbose = True 43 40 44 41 # … … 51 48 # Setup Test 52 49 ########################################################################## 53 def evolution_test(parallel=False, G = None ):50 def evolution_test(parallel=False, G = None, seq_interpolation_points=None): 54 51 55 52 # … … 65 62 # 66 63 if parallel: 67 if myid == 0 : print 'DISTRIBUTING PARALLEL DOMAIN'64 if myid == 0 and verbose : print 'DISTRIBUTING PARALLEL DOMAIN' 68 65 domain = distribute(domain, verbose=False) 69 66 … … 91 88 # 92 89 # Find which sub_domain in which the interpolation points are located 93 # 94 interpolation_points = [[0.4,0.51], [0.6,0.51], [0.8,0.51], [0.9,0.51]] 90 # 91 # Sometimes the interpolation points sit exactly 92 # between to centroids, so in the parallel run we 93 # reset the interpolation points to the centroids 94 # found in the sequential run 95 # 96 interpolation_points = [[0.4,0.5], [0.6,0.5], [0.8,0.5], [0.9,0.5]] 97 98 if parallel: 99 interpolation_points = seq_interpolation_points 100 101 95 102 gauge_values = [] 96 103 tri_ids = [] … … 108 115 tri_ids.append(2) 109 116 117 110 118 if verbose: print 'P%d has points = %s' %(myid, tri_ids) 111 119 112 113 120 if not parallel: 121 c_coord = domain.get_centroid_coordinates() 122 interpolation_points = [] 123 for id in tri_ids: 124 if id<1: 125 print 'ERROR: All interpolation points be within the sequential domain!' 126 interpolation_points.append(c_coord[id,:]) 127 114 128 # 115 129 # Evolve system through time … … 152 166 assert_(success) 153 167 154 return G 168 return G, interpolation_points 155 169 156 170 # Test an nprocsway run of the shallow water equations … … 182 196 # array G 183 197 # 184 pypar.barrier()198 barrier() 185 199 if myid == 0 and verbose: print 'SEQUENTIAL START' 186 200 187 G = evolution_test(parallel=False)201 G , interpolation_points = evolution_test(parallel=False) 188 202 G = num.array(G,num.float) 189 203 190 pypar.barrier()204 barrier() 191 205 192 206 # … … 196 210 if myid ==0 and verbose: print 'PARALLEL START' 197 211 198 evolution_test(parallel=True, G=G )212 evolution_test(parallel=True, G=G, seq_interpolation_points = interpolation_points) 199 213 200 214 
anuga_core/source/pymetis/pymetis/metis.c
r3414 r7452 1 1 #include <Python.h> 2 #include < Numeric/arrayobject.h>2 #include <numpy/arrayobject.h> 3 3 4 4 /* This must be the same as the metis idxtype */ … … 50 50 int edgecut; 51 51 int numflag = 0; // The metis routine requires an int * for numflag. 52 int dims[1]; // PyArray_FromDimsAndData needs an int[] of array sizes.52 npy_intp dims[1]; // PyArray_SimpleNewFromData needs an npy_intp[] of array sizes. 53 53 54 54 PyObject * elements; … … 104 104 105 105 dims[0] = ne; 106 epart_pyarr = (PyArrayObject *)PyArray_FromDimsAndData(1, dims, PyArray_INT, (char *)epart); 106 epart_pyarr = (PyArrayObject *)PyArray_SimpleNewFromData(1, dims, PyArray_INT, (void *)epart); 107 //epart_pyarr = (PyArrayObject *)PyArray_FromDimsAndData(1, dims, PyArray_INT, (char *)epart); 107 108 dims[0] = nn; 108 npart_pyarr = (PyArrayObject *)PyArray_FromDimsAndData(1, dims, PyArray_INT, (char *)npart); 109 npart_pyarr = (PyArrayObject *)PyArray_SimpleNewFromData(1, dims, PyArray_INT, (void *)npart); 110 //npart_pyarr = (PyArrayObject *)PyArray_FromDimsAndData(1, dims, PyArray_INT, (char *)npart); 109 111 110 112 
anuga_core/source/pymetis/test_metis.py
r5846 r7452 3 3 from sys import path 4 4 5 from Numericimport array, allclose5 from numpy import array, allclose 6 6 7 7 import unittest
Note: See TracChangeset
for help on using the changeset viewer.