Changeset 7870
- Timestamp:
- Jun 23, 2010, 2:33:54 PM (15 years ago)
- Location:
- trunk/anuga_core/source/anuga
- Files:
-
- 3 edited
-
file/test_csv.py (modified) (1 diff)
-
file/test_sww.py (modified) (1 diff)
-
shallow_water/shallow_water_domain.py (modified) (44 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_core/source/anuga/file/test_csv.py
r7866 r7870 3 3 import tempfile 4 4 import numpy as num 5 6 from anuga.utilities.system_tools import get_pathname_from_package 5 7 6 8 from csv_file import load_csv_as_array, load_csv_as_dict, store_parameters, \ -
trunk/anuga_core/source/anuga/file/test_sww.py
r7866 r7870 9 9 from anuga.shallow_water.shallow_water_domain import Domain 10 10 from sww import load_sww_as_domain, weed, get_mesh_and_quantities_from_file 11 from Scientific.IO.NetCDF import NetCDFFile 11 12 12 13 # boundary functions -
trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r7810 r7870 65 65 TRAC administration of ANUGA (User Manuals etc) at 66 66 https://datamining.anu.edu.au/anuga and Subversion repository at 67 $HeadURL$ 67 $HeadURL: https://datamining.anu.edu.au/svn/anuga/trunk/anuga_core/source/ 68 anuga/shallow_water/shallow_water_domain.py $ 68 69 69 70 Constraints: See GPL license in the user guide … … 74 75 """ 75 76 76 # Subversion keywords:77 #78 # $LastChangedDate$79 # $LastChangedRevision$80 # $LastChangedBy$81 82 77 83 78 import numpy as num … … 87 82 88 83 from anuga.shallow_water.forcing import Cross_section 89 from anuga.pmesh.mesh_interface import create_mesh_from_regions 90 from anuga.utilities.numerical_tools import gradient, mean, ensure_numeric 91 from anuga.geospatial_data.geospatial_data import ensure_geospatial 84 from anuga.utilities.numerical_tools import mean 92 85 from anuga.file.sww import SWW_file 93 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a 94 95 from anuga.fit_interpolate.interpolate import Modeltime_too_late, \ 96 Modeltime_too_early 97 98 from anuga.geometry.polygon import inside_polygon, polygon_area, \ 99 is_inside_polygon 86 100 87 import anuga.utilities.log as log 101 88 102 89 import types 103 from types import IntType, FloatType 104 from warnings import warn 105 106 107 ################################################################################ 108 # Shallow water domain 109 ################################################################################ 110 111 ## 112 # @brief Class for a shallow water domain. 90 113 91 class Domain(Generic_Domain): 114 115 #conserved_quantities = ['stage', 'xmomentum', 'ymomentum'] 116 #other_quantities = ['elevation', 'friction'] 117 118 ## 119 # @brief Instantiate a shallow water domain. 120 # @param coordinates 121 # @param vertices 122 # @param boundary 123 # @param tagged_elements 124 # @param geo_reference 125 # @param use_inscribed_circle 126 # @param mesh_filename 127 # @param use_cache 128 # @param verbose 129 # @param evolved_quantities 130 # @param full_send_dict 131 # @param ghost_recv_dict 132 # @param processor 133 # @param numproc 134 # @param number_of_full_nodes 135 # @param number_of_full_triangles 92 """ Class for a shallow water domain.""" 136 93 def __init__(self, 137 94 coordinates=None, … … 153 110 number_of_full_nodes=None, 154 111 number_of_full_triangles=None): 112 """ 113 Instantiate a shallow water domain. 114 coordinates - vertex locations for the mesh 115 vertices - vertex indices for the mesh 116 boundary - boundaries of the mesh 117 # @param tagged_elements 118 # @param geo_reference 119 # @param use_inscribed_circle 120 # @param mesh_filename 121 # @param use_cache 122 # @param verbose 123 # @param evolved_quantities 124 # @param full_send_dict 125 # @param ghost_recv_dict 126 # @param processor 127 # @param numproc 128 # @param number_of_full_nodes 129 # @param number_of_full_triangles 130 """ 155 131 156 132 # Define quantities for the shallow_water domain … … 165 141 166 142 Generic_Domain.__init__(self, 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 143 coordinates, 144 vertices, 145 boundary, 146 conserved_quantities, 147 evolved_quantities, 148 other_quantities, 149 tagged_elements, 150 geo_reference, 151 use_inscribed_circle, 152 mesh_filename, 153 use_cache, 154 verbose, 155 full_send_dict, 156 ghost_recv_dict, 157 processor, 158 numproc, 159 number_of_full_nodes=number_of_full_nodes, 160 number_of_full_triangles=number_of_full_triangles) 185 161 186 162 self.set_defaults() … … 200 176 201 177 202 203 204 ##205 # @brief Set default values, usually read in from a config file206 # @param flag207 178 def set_defaults(self): 208 179 """Set the default values in this routine. That way we can inherit class … … 212 183 from anuga.config import minimum_storable_height 213 184 from anuga.config import minimum_allowed_height, maximum_allowed_speed 214 from anuga.config import g, epsilon, beta_w, beta_w_dry,\185 from anuga.config import g, beta_w, beta_w_dry, \ 215 186 beta_uh, beta_uh_dry, beta_vh, beta_vh_dry, tight_slope_limiters 216 187 from anuga.config import alpha_balance … … 219 190 from anuga.config import use_edge_limiter 220 191 from anuga.config import use_centroid_velocities 221 222 223 192 224 193 self.set_minimum_allowed_height(minimum_allowed_height) … … 247 216 self.use_centroid_velocities = use_centroid_velocities 248 217 249 ## 250 # @brief 251 # @param flag 218 252 219 def set_new_mannings_function(self, flag=True): 253 220 """Cludge to allow unit test to pass, but to … … 263 230 264 231 265 ##266 # @brief267 # @param flag268 232 def set_use_edge_limiter(self, flag=True): 269 233 """Cludge to allow unit test to pass, but to … … 277 241 278 242 279 280 ##281 # @brief282 # @param beta283 243 def set_all_limiters(self, beta): 284 244 """Shorthand to assign one constant value [0,1] to all limiters. … … 299 259 self.quantities['ymomentum'].beta = beta 300 260 301 ## 302 # @brief 303 # @param flag 304 # @param reduction 261 305 262 def set_store_vertices_uniquely(self, flag, reduction=None): 306 263 """Decide whether vertex values should be stored uniquely as … … 406 363 this function will have no effect. 407 364 408 The format, where q is a list of names is for backwards compatibility only. 365 The format, where q is a list of names is for backwards compatibility 366 only. 409 367 It will take the specified quantities to be time dependent and assume 410 368 'elevation' to be static regardless. … … 416 374 return 417 375 418 # Check correc ness376 # Check correctness 419 377 for quantity_name in q: 420 378 msg = ('Quantity %s is not a valid conserved quantity' … … 422 380 assert quantity_name in self.quantities, msg 423 381 424 if type(q) == types.ListType:425 426 msg = 'List arguments to set_quantities_to_be_stored '427 msg += 'has been deprecated and will be removed in future '428 msg += 'versions of ANUGA.'429 msg += 'Please use dictionary argument instead'430 warn(msg, DeprecationWarning)431 432 433 434 # FIXME: Raise deprecation warning435 tmp = {}436 for x in q:437 tmp[x] = 2438 tmp['elevation'] = 1439 q = tmp440 441 382 assert type(q) == types.DictType 442 383 self.quantities_to_be_stored = q … … 567 508 568 509 569 ##570 # @brief Run integrity checks on shallow water domain.571 510 def check_integrity(self): 511 """ Run integrity checks on shallow water domain. """ 572 512 Generic_Domain.check_integrity(self) 573 513 … … 580 520 assert self.conserved_quantities[2] == 'ymomentum', msg 581 521 582 ##583 # @brief584 522 def extrapolate_second_order_sw(self): 585 #Call correct module function (either from this module or C-extension) 523 """Call correct module function 524 (either from this module or C-extension)""" 586 525 extrapolate_second_order_sw(self) 587 526 588 ##589 # @brief590 527 def compute_fluxes(self): 591 #Call correct module function (either from this module or C-extension) 528 """Call correct module function 529 (either from this module or C-extension)""" 592 530 compute_fluxes(self) 593 531 594 ##595 # @brief596 532 def distribute_to_vertices_and_edges(self): 597 # Call correct module function533 """ Call correct module function """ 598 534 if self.use_edge_limiter: 599 535 distribute_using_edge_limiter(self) … … 603 539 604 540 605 ##606 # @brief Evolve the model by one step.607 # @param yieldstep608 # @param finaltime609 # @param duration610 # @param skip_initial_step611 541 def evolve(self, 612 542 yieldstep=None, … … 614 544 duration=None, 615 545 skip_initial_step=False): 616 """Specialisation of basic evolve method from parent class""" 546 """Specialisation of basic evolve method from parent class. 547 548 Evolve the model by 1 step. 549 """ 617 550 618 551 # Call check integrity here rather than from user scripts … … 642 575 yield(t) 643 576 644 ## 645 # @brief 577 646 578 def initialise_storage(self): 647 579 """Create and initialise self.writer object for storing data. … … 655 587 self.writer.store_connectivity() 656 588 657 ## 658 # @brief 659 # @param name 589 660 590 def store_timestep(self): 661 591 """Store time dependent quantities and time. … … 667 597 self.writer.store_timestep() 668 598 669 ## 670 # @brief Get time stepping statistics string for printing. 671 # @param track_speeds 672 # @param triangle_id 599 673 600 def timestepping_statistics(self, 674 601 track_speeds=False, … … 712 639 Ch = Cw-Cz 713 640 714 s= ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\715 % (name.ljust(qwidth), Vh[0], Vh[1], Vh[2])716 717 s+= ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\718 % (name.ljust(qwidth), Eh[0], Eh[1], Eh[2])719 720 s+= ' %s: centroid_value = %.4f\n'\721 % (name.ljust(qwidth), Ch[0])722 723 msg += s641 message = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\ 642 % (name.ljust(qwidth), Vh[0], Vh[1], Vh[2]) 643 644 message += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\ 645 % (name.ljust(qwidth), Eh[0], Eh[1], Eh[2]) 646 647 message += ' %s: centroid_value = %.4f\n'\ 648 % (name.ljust(qwidth), Ch[0]) 649 650 msg += message 724 651 725 652 uh = self.quantities['xmomentum'] … … 739 666 Cu = Cuh/(Ch + epsilon) 740 667 name = 'U' 741 s = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\742 % (name.ljust(qwidth), Vu[0], Vu[1], Vu[2])743 744 s += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\745 % (name.ljust(qwidth), Eu[0], Eu[1], Eu[2])746 747 s += ' %s: centroid_value = %.4f\n'\748 % (name.ljust(qwidth), Cu[0])749 750 msg += s668 message = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n' \ 669 % (name.ljust(qwidth), Vu[0], Vu[1], Vu[2]) 670 671 message += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n' \ 672 % (name.ljust(qwidth), Eu[0], Eu[1], Eu[2]) 673 674 message += ' %s: centroid_value = %.4f\n' \ 675 % (name.ljust(qwidth), Cu[0]) 676 677 msg += message 751 678 752 679 Vv = Vvh/(Vh + epsilon) … … 754 681 Cv = Cvh/(Ch + epsilon) 755 682 name = 'V' 756 s = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\757 % (name.ljust(qwidth), Vv[0], Vv[1], Vv[2])758 759 s += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\760 % (name.ljust(qwidth), Ev[0], Ev[1], Ev[2])761 762 s+= ' %s: centroid_value = %.4f\n'\683 message = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n' \ 684 % (name.ljust(qwidth), Vv[0], Vv[1], Vv[2]) 685 686 message += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n' \ 687 % (name.ljust(qwidth), Ev[0], Ev[1], Ev[2]) 688 689 message += ' %s: centroid_value = %.4f\n'\ 763 690 %(name.ljust(qwidth), Cv[0]) 764 691 765 msg += s692 msg += message 766 693 767 694 # Froude number in each direction … … 771 698 Cfx = Cu/(num.sqrt(g*Ch) + epsilon) 772 699 773 s= ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\774 % (name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2])775 776 s+= ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\777 % (name.ljust(qwidth), Efx[0], Efx[1], Efx[2])778 779 s+= ' %s: centroid_value = %.4f\n'\780 % (name.ljust(qwidth), Cfx[0])781 782 msg += s700 message = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\ 701 % (name.ljust(qwidth), Vfx[0], Vfx[1], Vfx[2]) 702 703 message += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\ 704 % (name.ljust(qwidth), Efx[0], Efx[1], Efx[2]) 705 706 message += ' %s: centroid_value = %.4f\n'\ 707 % (name.ljust(qwidth), Cfx[0]) 708 709 msg += message 783 710 784 711 name = 'Froude (y)' … … 787 714 Cfy = Cv/(num.sqrt(g*Ch) + epsilon) 788 715 789 s= ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\790 % (name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2])791 792 s+= ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\793 % (name.ljust(qwidth), Efy[0], Efy[1], Efy[2])794 795 s+= ' %s: centroid_value = %.4f\n'\796 % (name.ljust(qwidth), Cfy[0])797 798 msg += s716 message = ' %s: vertex_values = %.4f,\t %.4f,\t %.4f\n'\ 717 % (name.ljust(qwidth), Vfy[0], Vfy[1], Vfy[2]) 718 719 message += ' %s: edge_values = %.4f,\t %.4f,\t %.4f\n'\ 720 % (name.ljust(qwidth), Efy[0], Efy[1], Efy[2]) 721 722 message += ' %s: centroid_value = %.4f\n'\ 723 % (name.ljust(qwidth), Cfy[0]) 724 725 msg += message 799 726 800 727 return msg … … 813 740 # Run through boundary array and compute for each segment 814 741 # the normal momentum ((uh, vh) dot normal) times segment length. 815 # Based on sign accumulate this into boundary_inflow and boundary_outflow. 742 # Based on sign accumulate this into boundary_inflow and 743 # boundary_outflow. 816 744 817 745 # Compute flows along boundary … … 842 770 if edge_flow > 0: 843 771 # Flow is inflow 844 total_boundary_inflow += edge_flow 772 total_boundary_inflow += edge_flow 845 773 else: 846 774 # Flow is outflow … … 881 809 882 810 area = self.mesh.get_areas() 883 volume = 0.0884 811 885 812 stage = self.get_quantity('stage').get_values(location='centroids') 886 elevation = self.get_quantity('elevation').get_values(location='centroids') 813 elevation = \ 814 self.get_quantity('elevation').get_values(location='centroids') 887 815 depth = stage-elevation 888 816 … … 897 825 total_boundary_outflow) = self.compute_boundary_flows() 898 826 899 s= '---------------------------\n'900 s+= 'Volumetric balance report:\n'901 s+= '--------------------------\n'902 s+= 'Total boundary inflow [m^3/s]: %.2f\n' % total_boundary_inflow903 s += 'Total boundary outflow [m^3/s]: %.2f\n' % total_boundary_outflow904 s+= 'Net boundary flow by tags [m^3/s]\n'827 message = '---------------------------\n' 828 message += 'Volumetric balance report:\n' 829 message += '--------------------------\n' 830 message += 'Total boundary inflow [m^3/s]: %.2f\n' % total_boundary_inflow 831 message += 'Total boundary outflow [m^3/s]: %.2f\n' % total_boundary_outflow 832 message += 'Net boundary flow by tags [m^3/s]\n' 905 833 for tag in boundary_flows: 906 s += ' %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag]) 907 908 s += 'Total net boundary flow [m^3/s]: %.2f\n' % (total_boundary_inflow + total_boundary_outflow) 909 s += 'Total volume in domain [m^3]: %.2f\n' % self.compute_total_volume() 910 911 # The go through explicit forcing update and record the rate of change for stage and 912 # record into forcing_inflow and forcing_outflow. Finally compute integral 913 # of depth to obtain total volume of domain. 834 message += ' %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag]) 835 836 message += 'Total net boundary flow [m^3/s]: %.2f\n' % \ 837 (total_boundary_inflow + total_boundary_outflow) 838 message += 'Total volume in domain [m^3]: %.2f\n' % \ 839 self.compute_total_volume() 840 841 # The go through explicit forcing update and record the rate of change 842 # for stage and 843 # record into forcing_inflow and forcing_outflow. Finally compute 844 # integral of depth to obtain total volume of domain. 914 845 915 846 # FIXME(Ole): This part is not yet done. 916 847 917 return s848 return message 918 849 919 850 ################################################################################ … … 951 882 from shallow_water_ext import compute_fluxes_ext_central \ 952 883 as compute_fluxes_ext 953 954 N = len(domain) # number_of_triangles955 884 956 885 # Shortcuts … … 1001 930 """Wrapper calling C version of extrapolate_second_order_sw""" 1002 931 1003 import sys1004 932 from shallow_water_ext import extrapolate_second_order_sw as extrapol2 1005 1006 N = len(domain) # number_of_triangles1007 933 1008 934 # Shortcuts … … 1066 992 domain.extrapolate_second_order_sw() 1067 993 else: 1068 raise 'Unknown order'994 raise Exception('Unknown order') 1069 995 else: 1070 996 # Old code: … … 1077 1003 Q.extrapolate_second_order_and_limit_by_vertex() 1078 1004 else: 1079 raise 'Unknown order'1005 raise Exception('Unknown order') 1080 1006 1081 1007 # Take bed elevation into account when water heights are small … … 1119 1045 Q.extrapolate_second_order_and_limit_by_edge() 1120 1046 else: 1121 raise 'Unknown order'1047 raise Exception('Unknown order') 1122 1048 1123 1049 balance_deep_and_shallow(domain) … … 1206 1132 elevation = domain.quantities['elevation'] 1207 1133 1208 h = stage.centroid_values - elevation.centroid_values 1209 z = elevation.vertex_values 1210 1211 x = domain.get_vertex_coordinates() 1212 g = domain.g 1213 1214 gravity_c(g, h, z, x, xmom_update, ymom_update) #, 1.0e-6) 1134 height = stage.centroid_values - elevation.centroid_values 1135 elevation = elevation.vertex_values 1136 1137 point = domain.get_vertex_coordinates() 1138 1139 gravity_c(domain.g, height, elevation, point, xmom_update, ymom_update) 1215 1140 1216 1141 ## … … 1241 1166 ymom_update = ymom.semi_implicit_update 1242 1167 1243 N = len(domain)1244 1168 eps = domain.minimum_allowed_height 1245 1169 g = domain.g 1246 1170 1247 1171 if domain.use_new_mannings: 1248 manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, ymom_update) 1172 manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, \ 1173 ymom_update) 1249 1174 else: 1250 manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, ymom_update) 1175 manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, \ 1176 ymom_update) 1251 1177 1252 1178 … … 1278 1204 ymom_update = ymom.explicit_update 1279 1205 1280 N = len(domain)1281 1206 eps = domain.minimum_allowed_height 1282 g = domain.g1283 1284 1207 1285 1208 if domain.use_new_mannings: 1286 manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, ymom_update) 1209 manning_friction_new(domain.g, eps, x, w, uh, vh, z, eta, xmom_update, \ 1210 ymom_update) 1287 1211 else: 1288 manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, ymom_update) 1289 1290 1291 1292 # FIXME (Ole): This was implemented for use with one of the analytical solutions (Sampson?) 1212 manning_friction_old(domain.g, eps, w, uh, vh, z, eta, xmom_update, \ 1213 ymom_update) 1214 1215 1216 1217 # FIXME (Ole): This was implemented for use with one of the analytical solutions 1293 1218 ## 1294 1219 # @brief Apply linear friction to a surface. … … 1301 1226 """ 1302 1227 1303 from math import sqrt1304 1305 1228 w = domain.quantities['stage'].centroid_values 1306 1229 z = domain.quantities['elevation'].centroid_values … … 1314 1237 ymom_update = domain.quantities['ymomentum'].semi_implicit_update 1315 1238 1316 N = len(domain) # number_of_triangles1239 num_tris = len(domain) 1317 1240 eps = domain.minimum_allowed_height 1318 g = domain.g #Not necessary? Why was this added? 1319 1320 for k in range(N): 1241 1242 for k in range(num_tris): 1321 1243 if tau[k] >= eps: 1322 1244 if h[k] >= eps: … … 1330 1252 surface_roughness_data, 1331 1253 verbose=False): 1332 """Returns an array of friction values for each wet element adjusted for depth. 1254 """Returns an array of friction values for each wet element adjusted for 1255 depth. 1333 1256 1334 1257 Inputs: 1335 1258 domain - computational domain object 1336 1259 default_friction - depth independent bottom friction 1337 surface_roughness_data - N x 5 array of n0, d1, n1, d2, n2 values for each1338 f riction region.1260 surface_roughness_data - N x 5 array of n0, d1, n1, d2, n2 values 1261 for each friction region. 1339 1262 1340 1263 Outputs: 1341 wet_friction - Array that can be used directly to update friction as follows: 1264 wet_friction - Array that can be used directly to update friction as 1265 follows: 1342 1266 domain.set_quantity('friction', wet_friction) 1343 1267 … … 1345 1269 1346 1270 """ 1347 1348 import numpy as num1349 1271 1350 # Create a temp array to store updated depth dependent friction for wet elements 1351 # EHR this is outwardly inneficient but not obvious how to avoid recreating each call?????? 1352 N=len(domain) 1353 wet_friction = num.zeros(N, num.float) 1354 wet_friction[:] = default_n0 # Initially assign default_n0 to all array so sure have no zeros values 1272 default_n0 = 0 # James - this was missing, don't know what it should be 1355 1273 1274 # Create a temp array to store updated depth dependent 1275 # friction for wet elements 1276 # EHR this is outwardly inneficient but not obvious how to avoid 1277 # recreating each call?????? 1278 1279 wet_friction = num.zeros(len(domain), num.float) 1280 wet_friction[:] = default_n0 # Initially assign default_n0 to all array so 1281 # sure have no zeros values 1356 1282 1357 depth = domain.create_quantity_from_expression('stage - elevation') # create depth instance for this timestep 1283 # create depth instance for this timestep 1284 depth = domain.create_quantity_from_expression('stage - elevation') 1358 1285 # Recompute depth as vector 1359 d = depth.get_values(location='centroids')1286 d_vals = depth.get_values(location='centroids') 1360 1287 1361 1288 # rebuild the 'friction' values adjusted for depth at this instant 1362 for i in domain.get_wet_elements(): # loop for each wet element in domain 1363 1289 # loop for each wet element in domain 1290 1291 for i in domain.get_wet_elements(): 1364 1292 # Get roughness data for each element 1365 n0 = float(surface_roughness_data[i,0]) 1366 d1 = float(surface_roughness_data[i,1]) 1367 n1 = float(surface_roughness_data[i,2]) 1368 d2 = float(surface_roughness_data[i,3]) 1369 n2 = float(surface_roughness_data[i,4]) 1293 d1 = float(surface_roughness_data[i, 1]) 1294 n1 = float(surface_roughness_data[i, 2]) 1295 d2 = float(surface_roughness_data[i, 3]) 1296 n2 = float(surface_roughness_data[i, 4]) 1370 1297 1371 1298 1372 1299 # Recompute friction values from depth for this element 1373 1300 1374 if d [i] <= d1:1375 d epth_dependent_friction= n11376 elif d [i] >= d2:1377 d epth_dependent_friction= n21301 if d_vals[i] <= d1: 1302 ddf = n1 1303 elif d_vals[i] >= d2: 1304 ddf = n2 1378 1305 else: 1379 d epth_dependent_friction = n1+((n2-n1)/(d2-d1))*(d[i]-d1)1306 ddf = n1 + ((n2-n1)/(d2-d1))*(d_vals[i]-d1) 1380 1307 1381 1308 # check sanity of result 1382 if (depth_dependent_friction < 0.010 or depth_dependent_friction > 9999.0) : 1383 log.critical('%s >>>> WARNING: computed depth_dependent friction ' 1309 if (ddf < 0.010 or \ 1310 ddf > 9999.0) : 1311 log.critical('>>>> WARNING: computed depth_dependent friction ' 1384 1312 'out of range, ddf%f, n1=%f, n2=%f' 1385 % (model_data.basename, 1386 depth_dependent_friction, n1, n2)) 1313 % (ddf, n1, n2)) 1387 1314 1388 1315 # update depth dependent friction for that wet element 1389 wet_friction[i] = depth_dependent_friction 1390 1391 # EHR add code to show range of 'friction across domain at this instant as sanity check????????? 1316 wet_friction[i] = ddf 1317 1318 # EHR add code to show range of 'friction across domain at this instant as 1319 # sanity check????????? 1392 1320 1393 1321 if verbose : 1394 nvals=domain.get_quantity('friction').get_values(location='centroids') # return array of domain nvals 1395 n_min=min(nvals) 1396 n_max=max(nvals) 1322 # return array of domain nvals 1323 nvals = domain.get_quantity('friction').get_values(location='centroids') 1324 n_min = min(nvals) 1325 n_max = max(nvals) 1397 1326 1398 1327 log.critical(' ++++ calculate_depth_dependent_friction - ' … … 1408 1337 ################################################################################ 1409 1338 1410 from anuga.utilities import compile 1411 if compile.can_use_C_extension('shallow_water_ext.c'): 1412 # Underlying C implementations can be accessed 1413 from shallow_water_ext import assign_windfield_values 1414 else: 1339 def _raise_compile_exception(): 1340 """ Raise exception if compiler not available. """ 1415 1341 msg = 'C implementations could not be accessed by %s.\n ' % __file__ 1416 1342 msg += 'Make sure compile_all.py has been run as described in ' 1417 1343 msg += 'the ANUGA installation guide.' 1418 raise Exception, msg 1419 1420 # Optimisation with psyco 1421 #from anuga.config import use_psyco 1422 #if use_psyco: 1423 #try: 1424 #import psyco 1425 #except: 1426 #import os 1427 #if os.name == 'posix' and os.uname()[4] in ['x86_64', 'ia64']: 1428 #pass 1429 ##Psyco isn't supported on 64 bit systems, but it doesn't matter 1430 #else: 1431 #msg = ('WARNING: psyco (speedup) could not be imported, ' 1432 #'you may want to consider installing it') 1433 #log.critical(msg) 1434 #else: 1435 #psyco.bind(Domain.distribute_to_vertices_and_edges) 1436 #psyco.bind(Domain.compute_fluxes) 1437 1344 raise Exception(msg) 1345 1346 from anuga.utilities import compile 1347 if not compile.can_use_C_extension('shallow_water_ext.c'): 1348 _raise_compile_exception() 1438 1349 1439 1350 if __name__ == "__main__":
Note: See TracChangeset
for help on using the changeset viewer.
