Changeset 5742
- Timestamp:
- Sep 6, 2008, 4:30:16 PM (16 years ago)
- Location:
- anuga_work/development/anuga_1d
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/development/anuga_1d/dam_h_elevation.py
r5741 r5742 1 1 import os 2 2 from math import sqrt, sin, cos, pi, exp 3 from shallow_water_domain _newimport *3 from shallow_water_domain import * 4 4 from Numeric import zeros, Float 5 5 #from analytic_dam_sudi import AnalyticDam … … 31 31 32 32 L=2000.0 33 N= 40033 N=200 34 34 35 35 cell_len=L/N … … 61 61 y=zeros(len(x), Float) 62 62 for i in range(len(x)): 63 if 1100.0 <= x[i] <=1200.0:63 if 800.0 <= x[i] <=1200.0: 64 64 y[i]=10.0 65 65 else: 66 y[i]= 10.066 y[i]=5.0 67 67 return y 68 68 … … 108 108 109 109 domain.set_quantity('stage',stage_perturb) 110 domain.set_quantity('elevation', elevation_ parabol)110 domain.set_quantity('elevation', elevation_box) 111 111 #domain.set_quantity('xmomentum', xmom_sincos) 112 112 domain.order=domain.default_order … … 152 152 153 153 plot(X,ElevationQ,X,HeightQ) 154 plot1.set_ylim([9.99,10.01])154 #plot1.set_ylim([9.99,10.01]) 155 155 xlabel('Position') 156 156 ylabel('Stage') … … 180 180 181 181 raw_input("Press any key") 182 -
anuga_work/development/anuga_1d/domain.py
r5741 r5742 57 57 58 58 self.max_speed_array = zeros(N, Float) 59 59 60 60 61 self.normals = zeros((N, 2), Float) … … 875 876 self.number_of_steps = 0 876 877 self.number_of_first_order_steps = 0 877 self.max_speed = zeros(N, Float)878 self.max_speed_array = 0.0 878 879 879 880 -
anuga_work/development/anuga_1d/quantity.py
r5741 r5742 860 860 # (either from this module or C-extension) 861 861 #backup_centroid_values(self) 862 N = self.domain.number_of_elements 863 for i in range(N): 864 self.centroid_backup_values[i] = self.centroid_values[i] 862 863 self.centroid_backup_values[:] = (self.centroid_values).astype('f') 865 864 866 865 def saxpy_centroid_values(self,a,b): … … 868 867 # (either from this module or C-extension) 869 868 #saxpy_centroid_values(self,a,b) 870 N = self.domain.number_of_elements 871 for i in range(N): 872 self.centroid_backup_values[i] = a*self.centroid_values[i] + b*self.centroid_backup_values[i] 873 869 self.centroid_backup_values[:] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f') 870 874 871 class Conserved_quantity(Quantity): 875 872 """Class conserved quantity adds to Quantity: -
anuga_work/development/anuga_1d/shallow_water_domain_new.py
r5741 r5742 49 49 #Shallow water domain 50 50 class Domain(Generic_Domain): 51 51 52 52 def __init__(self, coordinates, boundary = None, tagged_elements = None): 53 53 … … 60 60 from config import minimum_allowed_height, g 61 61 self.minimum_allowed_height = minimum_allowed_height 62 self.h0 = minimum_allowed_height 62 63 self.g = g 63 64 … … 88 89 self.quantities_to_be_stored = ['stage', 'xmomentum', 'elevation'] 89 90 90 self.__doc__ = 'shallow_water_domain_new' 91 91 self.__doc__ = 'shallow_water_domain_new' 92 92 93 93 def set_quantities_to_be_stored(self, q): -
anuga_work/development/anuga_1d/test_shallow_water.py
r5741 r5742 41 41 #print domain.quantities['stage'].explicit_update 42 42 #print domain.quantities['xmomentum'].explicit_update 43 print stage_ud43 #print stage_ud 44 44 45 45 assert allclose( domain.quantities['stage'].explicit_update, stage_ud ) … … 101 101 Compare shallow_water_domain gravity calculation 102 102 """ 103 104 103 104 def slope_one(x): 105 105 return x 106 106 107 107 domain = Domain(self.points) 108 108 domain.set_quantity('stage',4.0) 109 109 domain.set_quantity('elevation',slope_one) 110 110 domain.set_boundary({'exterior' : Reflective_boundary(domain)}) 111 111 112 112 gravity(domain) 113 113 114 115 114 #print domain.quantities['stage'].vertex_values 115 #print domain.quantities['elevation'].vertex_values 116 116 #print domain.quantities['xmomentum'].explicit_update 117 117 … … 123 123 Compare still lake solution for various versions of shallow_water_domain 124 124 """ 125 126 125 126 def slope_square(x): 127 127 return maximum(4.0-(x-5.0)*(x-5.0), 0.0) 128 128 129 129 domain = Domain(self.points2) 130 130 domain.set_quantity('stage',10.0) 131 131 domain.set_quantity('elevation',slope_square) 132 132 domain.set_boundary({'exterior' : Reflective_boundary(domain)}) 133 133 … … 141 141 pass 142 142 143 ## print domain.quantities['stage'].vertex_values144 ## print domain.quantities['elevation'].vertex_values145 ## print domain.quantities['xmomentum'].vertex_values146 ## 147 ##148 ## print domain.quantities['stage'].centroid_values149 ## print domain.quantities['elevation'].centroid_values150 ## print domain.quantities['xmomentum'].centroid_values143 ## print domain.quantities['stage'].vertex_values 144 ## print domain.quantities['elevation'].vertex_values 145 ## print domain.quantities['xmomentum'].vertex_values 146 ## 147 ## 148 ## print domain.quantities['stage'].centroid_values 149 ## print domain.quantities['elevation'].centroid_values 150 ## print domain.quantities['xmomentum'].centroid_values 151 151 152 152 assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values ) … … 158 158 Compare still lake solution for various versions of shallow_water_domain 159 159 """ 160 161 160 161 def slope_square(x): 162 162 return maximum(4.0-(x-5.0)*(x-5.0), 0.0) 163 163 164 164 domain = Domain(self.points2) 165 165 domain.set_quantity('stage',10.0) 166 166 domain.set_quantity('elevation',slope_square) 167 167 domain.set_boundary({'exterior' : Reflective_boundary(domain)}) 168 168 … … 182 182 Compare still lake solution for various versions of shallow_water_domain 183 183 """ 184 185 184 185 def slope_square(x): 186 186 return maximum(4.0-(x-5.0)*(x-5.0), 0.0) 187 187 188 188 domain = Domain(self.points2) 189 189 domain.set_quantity('stage',10.0) 190 190 domain.set_quantity('elevation',slope_square) 191 191 domain.set_boundary({'exterior' : Reflective_boundary(domain)}) 192 192 … … 238 238 Stage = domain.quantities['stage'] 239 239 Xmom = domain.quantities['xmomentum'] 240 # Ymom = domain.quantities['ymomentum']240 # Ymom = domain.quantities['ymomentum'] 241 241 Bed = domain.quantities['elevation'] 242 242 … … 244 244 #stage = Stage.edge_values 245 245 #xmom = Xmom.edge_values 246 # ymom = Ymom.edge_values246 # ymom = Ymom.edge_values 247 247 #bed = Bed.edge_values 248 248 … … 259 259 #print 'stage_bdry',stage_bdry 260 260 #print 'xmom_bdry', xmom_bdry 261 # ymom_bdry = Ymom.boundary_values262 263 # flux = zeros(3, Float) #Work array for summing up fluxes261 # ymom_bdry = Ymom.boundary_values 262 263 # flux = zeros(3, Float) #Work array for summing up fluxes 264 264 flux = zeros(2, Float) #Work array for summing up fluxes 265 265 ql = zeros(2, Float) -
anuga_work/development/anuga_1d/util_ext.h
r5563 r5742 279 279 } 280 280 281 282 281 283 double get_python_double(PyObject *O, char *name) { 282 284 PyObject *TObject; 285 #define BUFFER_SIZE 80 286 char buf[BUFFER_SIZE]; 283 287 double tmp; 288 int n; 284 289 285 290 … … 287 292 TObject = PyObject_GetAttrString(O, name); 288 293 if (!TObject) { 289 PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_double could not obtain double from object"); 294 n = snprintf(buf, BUFFER_SIZE, "util_ext.h: get_python_double could not obtain double %s.\n", name); 295 //printf("name = %s",name); 296 PyErr_SetString(PyExc_RuntimeError, buf); 297 290 298 return 0.0; 291 299 } … … 298 306 } 299 307 308 309 310 300 311 int get_python_integer(PyObject *O, char *name) { 301 312 PyObject *TObject; 302 int tmp; 313 #define BUFFER_SIZE 80 314 char buf[BUFFER_SIZE]; 315 long tmp; 316 int n; 303 317 304 318 … … 306 320 TObject = PyObject_GetAttrString(O, name); 307 321 if (!TObject) { 308 PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_integer could not obtain double from object"); 322 n = snprintf(buf, BUFFER_SIZE, "util_ext.h: get_python_integer could not obtain double %s.\n", name); 323 //printf("name = %s",name); 324 PyErr_SetString(PyExc_RuntimeError, buf); 309 325 return 0; 310 326 } 311 327 312 tmp = Py Float_AsDouble(TObject);328 tmp = PyInt_AsLong(TObject); 313 329 314 330 Py_DECREF(TObject);
Note: See TracChangeset
for help on using the changeset viewer.