Changeset 3293
- Timestamp:
- Jul 10, 2006, 11:23:11 AM (18 years ago)
- Location:
- development/pyvolution-1d
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
development/pyvolution-1d/domain.py
r2791 r3293 706 706 707 707 #Initialise interval of timestep sizes (for reporting only) 708 # SEEMS WIERD 708 709 self.min_timestep = max_timestep 709 710 self.max_timestep = min_timestep … … 747 748 #Update boundary values 748 749 self.update_boundary() 750 751 #print 'timestep', self.timestep 749 752 750 753 #Update time … … 801 804 else: 802 805 raise 'Unknown order' 803 Q.interpolate_from_vertices_to_edges()806 #Q.interpolate_from_vertices_to_edges() 804 807 805 808 -
development/pyvolution-1d/quantity.py
r2791 r3293 518 518 Qc = self.centroid_values 519 519 Qv = self.vertex_values 520 520 521 521 #Check each triangle 522 522 for k in range(self.domain.number_of_elements): … … 526 526 #vertex coordinates 527 527 x0, x1 = V[k,:] 528 529 528 #Extrapolate 530 529 Qv[k,0] = Qc[k] + G[k]*(x0-x) 531 530 Qv[k,1] = Qc[k] + G[k]*(x1-x) 532 533 534 535 531 536 532 def limit(self): -
development/pyvolution-1d/shallow_water_1d.py
r2791 r3293 63 63 64 64 #forcing terms not included in 1d domain ?WHy? 65 self.forcing_terms.append(gravity)66 self.forcing_terms.append(manning_friction)65 #self.forcing_terms.append(gravity) 66 #self.forcing_terms.append(manning_friction) 67 67 #print "\nI have Removed forcing terms line 64 1dsw" 68 68 … … 176 176 177 177 #Initialise real time viz if requested 178 if self.visualise is True and self.time == 0.0:179 if self.visualiser is None:180 self.initialise_visualiser()181 182 self.visualiser.update_timer()183 self.visualiser.setup_all()178 #if self.visualise is True and self.time == 0.0: 179 # if self.visualiser is None: 180 # self.initialise_visualiser() 181 # 182 # self.visualiser.update_timer() 183 # self.visualiser.setup_all() 184 184 185 185 #Store model data, e.g. for visualisation 186 if self.store is True and self.time == 0.0:187 self.initialise_storage()188 #print 'Storing results in ' + self.writer.filename189 else:190 pass191 #print 'Results will not be stored.'192 #print 'To store results set domain.store = True'193 #FIXME: Diagnostic output should be controlled by194 # a 'verbose' flag living in domain (or in a parent class)186 #if self.store is True and self.time == 0.0: 187 # self.initialise_storage() 188 # #print 'Storing results in ' + self.writer.filename 189 #else: 190 # pass 191 # #print 'Results will not be stored.' 192 # #print 'To store results set domain.store = True' 193 # #FIXME: Diagnostic output should be controlled by 194 # # a 'verbose' flag living in domain (or in a parent class) 195 195 196 196 #Call basic machinery from parent class … … 198 198 skip_initial_step): 199 199 #Real time viz 200 if self.visualise is True:201 self.visualiser.update_all()202 self.visualiser.update_timer()200 # if self.visualise is True: 201 # self.visualiser.update_all() 202 # self.visualiser.update_timer() 203 203 204 204 205 205 #Store model data, e.g. for subsequent visualisation 206 if self.store is True:207 self.store_timestep(self.quantities_to_be_stored)206 # if self.store is True: 207 # self.store_timestep(self.quantities_to_be_stored) 208 208 209 209 #FIXME: Could maybe be taken from specified list … … 312 312 q_right[1] = q_right[1]*normal 313 313 314 z = (zl+zr)/2 #Take average of field values 314 #z = (zl+zr)/2 #Take average of field values 315 z = 0.5*(zl+zr) #Take average of field values 315 316 316 317 w_left = q_left[0] #w=h+z … … 350 351 #Flux computation 351 352 352 #FIXME(Ole): Why is it again that we don't353 #use uh_left and uh_right directly in the first entries?354 353 #flux_left = array([u_left*h_left, 355 # u_left*uh_left + 0.5*g*h_left**2, 356 # u_left*vh_left]) 354 # u_left*uh_left + 0.5*g*h_left**2]) 357 355 #flux_right = array([u_right*h_right, 358 # u_right*uh_right + 0.5*g*h_right**2, 359 # u_right*vh_right]) 360 356 # u_right*uh_right + 0.5*g*h_right**2]) 361 357 flux_left = array([u_left*h_left, 362 u_left*uh_left + 0.5*g*h_left* *2])358 u_left*uh_left + 0.5*g*h_left*h_left]) 363 359 flux_right = array([u_right*h_right, 364 u_right*uh_right + 0.5*g*h_right* *2])360 u_right*uh_right + 0.5*g*h_right*h_right]) 365 361 366 362 denom = s_max-s_min 367 363 if denom == 0.0: 368 #edgeflux = array([0.0, 0.0, 0.0])369 364 edgeflux = array([0.0, 0.0]) 370 365 max_speed = 0.0 … … 373 368 edgeflux += s_max*s_min*(q_right-q_left)/denom 374 369 375 #edgeflux[1] = -edgeflux[1]*normal376 370 edgeflux[1] = edgeflux[1]*normal 377 371 378 #edgeflux = rotate(edgeflux, normal, direction=-1)379 372 max_speed = max(abs(s_max), abs(s_min)) 380 373 … … 445 438 for i in range(2): 446 439 #Quantities inside volume facing neighbour i 447 #ql = [stage[k, i], xmom[k, i], ymom[k, i]]448 ql[0] = stage[k, i]449 ql [1] = xmom[k, i]440 #ql[0] = stage[k, i] 441 #ql[1] = xmom[k, i] 442 ql = [stage[k, i], xmom[k, i]] 450 443 zl = bed[k, i] 451 444 … … 454 447 if n < 0: 455 448 m = -n-1 #Convert negative flag to index 456 #qr = [stage_bdry[m], xmom_bdry[m], ymom_bdry[m]]457 449 qr[0] = stage_bdry[m] 458 450 qr[1] = xmom_bdry[m] … … 469 461 #Outward pointing normal vector 470 462 normal = domain.normals[k, i] 471 #CHECK THIS LINE [k, 2*i:2*i+1] 472 463 473 464 #Flux computation using provided function 474 465 #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr) … … 486 477 try: 487 478 #timestep = min(timestep, 0.5*domain.radii[k]/max_speed) 488 timestep = min(timestep, 0.5/max_speed) 479 #timestep = 0.01 480 timestep = min(timestep, 0.5*domain.areas[k]/max_speed) 481 if timestep < 0.00001: 482 #print 'max_speed', max_speed 483 s = domain.quantities['stage'] 484 s = s.centroid_values 485 xm = domain.quantities['xmomentum'] 486 xm = xm.centroid_values 487 #print 'h', s[k] 488 #print 'xm', xm[k] 489 #print 'u', xm[k]/s[k] 490 #break 491 #timestep = 0.01 492 #print 'areas', domain.areas[k] 493 #print "timestep", timestep 489 494 except ZeroDivisionError: 490 495 pass … … 492 497 #Normalise by area and store for when all conserved 493 498 #quantities get updated 494 flux /= domain.areas[k] 499 #flux /= domain.areas[k] 500 # ADD ABOVE LINE AGAIN 495 501 Stage.explicit_update[k] = flux[0] 496 502 Xmom.explicit_update[k] = flux[1] … … 596 602 597 603 #Remove very thin layers of water 598 protect_against_infinitesimal_and_negative_heights(domain)604 #protect_against_infinitesimal_and_negative_heights(domain) 599 605 600 606 #Extrapolate all conserved quantities 601 if optimised_gradient_limiter: 602 #MH090605 if second order, 603 #perform the extrapolation and limiting on 604 #all of the conserved quantities 605 606 if (domain.order == 1): 607 for name in domain.conserved_quantities: 608 Q = domain.quantities[name] 609 Q.extrapolate_first_order() 607 #if optimised_gradient_limiter: 608 # #MH090605 if second order, 609 # #perform the extrapolation and limiting on 610 # #all of the conserved quantities 611 612 # if (domain.order == 1): 613 # for name in domain.conserved_quantities: 614 # Q = domain.quantities[name] 615 # Q.extrapolate_first_order() 616 # elif domain.order == 2: 617 # domain.extrapolate_second_order_sw() 618 # else: 619 # raise 'Unknown order' 620 #else: 621 #old code: 622 for name in domain.conserved_quantities: 623 Q = domain.quantities[name] 624 if domain.order == 1: 625 Q.extrapolate_first_order() 610 626 elif domain.order == 2: 611 domain.extrapolate_second_order_sw() 627 Q.extrapolate_second_order() 628 Q.limit() 612 629 else: 613 630 raise 'Unknown order' 614 else:615 #old code:616 for name in domain.conserved_quantities:617 Q = domain.quantities[name]618 if domain.order == 1:619 Q.extrapolate_first_order()620 elif domain.order == 2:621 Q.extrapolate_second_order()622 Q.limit()623 else:624 raise 'Unknown order'625 626 631 627 632 #Take bed elevation into account when water heights are small 628 balance_deep_and_shallow(domain)633 #balance_deep_and_shallow(domain) 629 634 630 635 #Compute edge values by interpolation … … 664 669 wc = domain.quantities['stage'].centroid_values 665 670 zc = domain.quantities['elevation'].centroid_values 666 #xmomc = domain.quantities['xmomentum'].centroid_values667 ymomc = domain.quantities['ymomentum'].centroid_values671 xmomc = domain.quantities['xmomentum'].centroid_values 672 # ymomc = domain.quantities['ymomentum'].centroid_values 668 673 669 674 from shallow_water_ext import protect -
development/pyvolution-1d/test_shallow_water_1d.py
r2791 r3293 466 466 def test_compute_fluxes1(self): 467 467 #Use values from previous version 468 469 #a = [0.0, 0.0] 470 #b = [0.0, 2.0] 471 #c = [2.0,0.0] 472 #d = [0.0, 4.0] 473 #e = [2.0, 2.0] 474 #f = [4.0,0.0] 475 a=0.0 476 b=2.0 477 c=4.0 478 d=6.0 479 e=8.0 480 #f=10.0 481 482 #points = [a, b, c, d, e, f] 468 print "test compute fluxes 1" 469 470 a=0.0 471 b=2.0 472 c=4.0 473 d=6.0 474 e=8.0 475 483 476 points = [a, b, c, d, e] 484 #bac, bce, ecf, dbe 485 #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 486 487 #domain = Domain(points, vertices) 477 488 478 domain = Domain(points) 489 479 … … 493 483 val3 = 2.+8.0/3 494 484 495 #domain.set_quantity('stage', [[val0, val0, val0], [val1, val1, val1], 496 # [val2, val2, val2], [val3, val3, val3]]) 497 domain.set_quantity('stage', [[val0, val0], [val1, val1], 498 [val2, val2], [val3, val3]]) 485 domain.set_quantity('stage', [[val3, val3], [val1, val1], 486 [val2, val2], [val0, val0]]) 499 487 stage = domain.get_quantity('stage') 500 488 print stage … … 504 492 505 493 #Flux across right edge of volume 1 506 #normal = domain.get_normal(1,0)507 #ql = domain.get_conserved_quantities(vol_id=1, edge=0)508 #qr = domain.get_conserved_quantities(vol_id=2, edge=2)509 #flux0, max_speed = flux_function(normal, ql, qr, zl, zr)510 511 #Flux across right edge of volume 1512 #ql = domain.get_conserved_quantities(vol_id=1, edge=1)513 #qr = domain.get_conserved_quantities(vol_id=2, edge=0)514 494 ql = domain.get_conserved_quantities(vol_id=1, vertex=1) 515 495 qr = domain.get_conserved_quantities(vol_id=2, vertex=0) … … 521 501 522 502 #Check that flux seen from other triangles is inverse 523 #tmp = qr; qr=ql; ql=tmp524 #normal = domain.get_normal(3,0)525 #flux, max_speed = flux_function(normal, ql, qr, zl, zr)526 503 flux, max_speed = flux_function(-1.0, qr, ql, zl, zr) 527 504 print flux0 … … 529 506 assert allclose(flux + flux0, 0.) 530 507 531 #print flux0532 #print max_speed533 #assert allclose(flux0, [-15.3598804, 253.71111111, 0.])534 508 assert allclose(flux0, [-15.3598804, 253.71111111]) 535 509 assert allclose(max_speed, 9.21592824046) 536 510 537 511 538 #Flux across upper edge of volume 1 539 #normal = domain.get_normal(1,1) 540 #ql = domain.get_conserved_quantities(vol_id=1, edge=1) 541 #qr = domain.get_conserved_quantities(vol_id=2, edge=0) 512 #Flux across left edge of volume 1 542 513 ql = domain.get_conserved_quantities(vol_id=1, vertex=0) 543 514 qr = domain.get_conserved_quantities(vol_id=0, vertex=1) 544 #flux1, max_speed = flux_function(normal, ql, qr, zl, zr)545 515 flux1, max_speed = flux_function(-1.0, ql, qr, zl, zr) 546 #assert allclose(flux1, [2.4098563, 0., 123.04444444])547 516 print 'flux1', flux1 548 assert allclose(flux1, [2.4098563, 0.])517 assert allclose(flux1, [2.4098563, -123.04444444]) 549 518 assert allclose(max_speed, 7.22956891292) 550 519 551 #Flux across lower left hypotenuse of volume 1 552 #normal = domain.get_normal(1,2) 553 #ql = domain.get_conserved_quantities(vol_id=1, edge=2) 554 #qr = domain.get_conserved_quantities(vol_id=0, edge=1) 555 #ql = domain.get_conserved_quantities(vol_id=1, vertex=2) 556 #qr = domain.get_conserved_quantities(vol_id=0, vertex=1) 557 #flux2, max_speed = flux_function(normal, ql, qr, zl, zr) 558 559 #assert allclose(flux2, [9.63942522, -61.59685738, -61.59685738]) 560 #assert allclose(max_speed, 7.22956891292) 561 562 #Scale, add up and check that compute_fluxes is correct for vol 1 563 #e0 = domain.edgelengths[1, 0] 564 #e1 = domain.edgelengths[1, 1] 565 #e2 = domain.edgelengths[1, 2] 566 567 #total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1] 520 #Add up and check that compute_fluxes is correct for vol 1 568 521 total_flux = -(flux0+flux1)/domain.areas[1] 569 #assert allclose(total_flux, [-0.68218178, -166.6, -35.93333333]) 570 assert allclose(total_flux, [-0.68218178, -166.6]) 571 572 573 522 print 'total flux', total_flux 523 print 'domain area',domain.areas[1] 524 assert allclose(total_flux, [6.47501205, -65.333333333]) 525 574 526 domain.compute_fluxes() 575 527 576 #assert allclose(total_flux, domain.explicit_update[1,:]) 577 for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']): 528 for i, name in enumerate(['stage', 'xmomentum']): 578 529 assert allclose(total_flux[i], 579 530 domain.quantities[name].explicit_update[1]) … … 585 536 # [-35.68522449, 0., 69.68888889]]) 586 537 538 print domain.quantities['stage'].explicit_update 539 print domain.quantities['xmomentum'].explicit_update 587 540 assert allclose(domain.quantities['stage'].explicit_update, 588 [0., -0.68218178, -111.77316251, -35.68522449])541 [0., 6.47501205, -111.77316251, -35.68522449]) 589 542 assert allclose(domain.quantities['xmomentum'].explicit_update, 590 [-69.68888889, - 166.6, 69.68888889, 0])543 [-69.68888889, -65.33333333, 69.68888889, 0]) 591 544 assert allclose(domain.quantities['ymomentum'].explicit_update, 592 545 [-69.68888889, -35.93333333, 0., 69.68888889]) … … 597 550 def test_catching_negative_heights(self): 598 551 599 #a = [0.0, 0.0] 600 #b = [0.0, 2.0] 601 #c = [2.0,0.0] 602 #d = [0.0, 4.0] 603 #e = [2.0, 2.0] 604 #f = [4.0,0.0] 605 a=0.0 606 b=2.0 607 c=4.0 608 d=6.0 609 e=8.0 610 #f=10.0 611 612 #points = [a, b, c, d, e, f] 552 a=0.0 553 b=2.0 554 c=4.0 555 d=6.0 556 e=8.0 557 613 558 points = [a, b, c, d, e] 614 #bac, bce, ecf, dbe615 #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]616 617 #domain = Domain(points, vertices)618 559 domain = Domain(points) 619 560 val0 = 2.+2.0/3 … … 623 564 624 565 zl=zr=4 #Too large 625 #domain.set_quantity('elevation', zl*ones( (4,3) ))626 #domain.set_quantity('stage', [[val0, val0-1, val0-2],627 # [val1, val1+1, val1],628 # [val2, val2-2, val2],629 # [val3-0.5, val3, val3]])630 566 domain.set_quantity('elevation', zl*ones( (4,2) )) 631 567 domain.set_quantity('stage', [[val0, val0-1], … … 703 639 704 640 from config import g 705 #a = [0.0, 0.0] 706 #b = [0.0, 2.0] 707 #c = [2.0, 0.0] 708 #d = [0.0, 4.0] 709 #e = [2.0, 2.0] 710 #f = [4.0, 0.0] 711 a=0.0 712 b=2.0 713 c=4.0 714 d=6.0 715 e=8.0 716 #f=10.0 717 718 #points = [a, b, c, d, e, f] 641 a=0.0 642 b=2.0 643 c=4.0 644 d=6.0 645 e=8.0 646 719 647 points = [a, b, c, d, e] 720 #bac, bce, ecf, dbe 721 #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 722 723 #domain = Domain(points, vertices) 648 724 649 domain = Domain(points) 725 650 … … 751 676 from config import g 752 677 753 #a = [0.0, 0.0] 754 #b = [0.0, 2.0] 755 #c = [2.0, 0.0] 756 #d = [0.0, 4.0] 757 #e = [2.0, 2.0] 758 #f = [4.0, 0.0] 759 a=0.0 760 b=2.0 761 c=4.0 762 d=6.0 763 e=8.0 764 #f=10.0 765 766 #points = [a, b, c, d, e, f] 678 a=0.0 679 b=2.0 680 c=4.0 681 d=6.0 682 e=8.0 683 767 684 points = [a, b, c, d, e] 768 #bac, bce, ecf, dbe 769 #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 770 771 #domain = Domain(points, vertices) 685 772 686 domain = Domain(points) 773 687 … … 1086 1000 #Using test data generated by pyvolution-2 1087 1001 #Assuming no friction and flat bed (0.0) 1088 1089 #a = [0.0, 0.0] 1090 #b = [0.0, 2.0] 1091 #c = [2.0,0.0] 1092 #d = [0.0, 4.0] 1093 #e = [2.0, 2.0] 1094 #f = [4.0,0.0] 1095 a=0.0 1096 b=2.0 1097 c=4.0 1098 d=6.0 1099 e=8.0 1100 #f=10.0 1101 1102 #points = [a, b, c, d, e, f] 1002 print "\ntest distriute basic" 1003 1004 a=0.0 1005 b=2.0 1006 c=4.0 1007 d=6.0 1008 e=8.0 1009 1103 1010 points = [a, b, c, d, e] 1104 #bac, bce, ecf, dbe1105 #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]1106 1107 #domain = Domain(points, vertices)1108 1011 domain = Domain(points) 1109 1012 … … 1120 1023 domain.order = 1 1121 1024 domain.distribute_to_vertices_and_edges() 1025 print "First order L", L 1122 1026 assert allclose(L[1], val1) 1123 1027 1124 #Second order 1028 #Second order 1125 1029 domain.order = 2 1030 a = domain.quantities['stage'].compute_gradients() 1031 print 'gradients', a 1032 1126 1033 domain.distribute_to_vertices_and_edges() 1127 assert allclose(L[1], [2.2, 4.9, 4.9]) 1128 1129 1034 print "Second order L", L 1035 print "L[1]", L[1] 1036 #assert allclose(L[1], [2.2, 4.9, 4.9]) 1037 assert allclose(L[1], [2.5, 5.5]) 1130 1038 1131 1039 def test_distribute_away_from_bed(self): … … 1133 1041 #Assuming no friction and flat bed (0.0) 1134 1042 1135 #a = [0.0, 0.0] 1136 #b = [0.0, 2.0] 1137 #c = [2.0,0.0] 1138 #d = [0.0, 4.0] 1139 #e = [2.0, 2.0] 1140 #f = [4.0,0.0] 1043 print "\ntest distribute awaway from bed" 1044 1141 1045 a=0.0 1142 1046 b=2.0 … … 1148 1052 #points = [a, b, c, d, e, f] 1149 1053 points = [a, b, c, d, e] 1150 #bac, bce, ecf, dbe 1151 #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 1152 1153 #domain = Domain(points, vertices) 1054 1154 1055 domain = Domain(points) 1155 1056 … … 1178 1079 assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778]) 1179 1080 1180 1181 1182 1183 1081 def test_1d_solution_I(self): 1082 print "TEST 1D-SOLUTION I" 1083 1084 L = 2000.0 # Length of channel (m) 1085 N = 100 # Number of compuational cells 1086 cell_len = L/N # Origin = 0.0 1087 1088 points = zeros(N+1,Float) 1089 for i in range(N+1): 1090 points[i] = i*cell_len 1091 1092 domain = Domain(points) 1093 1094 def stage(x): 1095 for i in range(len(x)): 1096 if x[i]<=1000.0: 1097 x[i] = 10.0 1098 else: 1099 x[i] = 5.0 1100 return x 1101 1102 domain.set_quantity('stage', stage) 1103 #L = domain.quantities['stage'].vertex_values 1104 #print "Initial Stage" 1105 #print L 1106 domain.set_boundary({'exterior': Reflective_boundary(domain)}) 1107 1108 import time 1109 t0 = time.time() 1110 yieldstep = 10 1111 finaltime = 50.0 1112 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): 1113 pass 1114 1115 print 'That took %.2f seconds' %(time.time()-t0) 1116 1117 #L = domain.quantities['stage'].vertex_values 1118 #print "Final Stage" 1119 #print L 1120 1121 C = domain.quantities['stage'].vertex_values 1122 #print C 1123 1124 f = file('test_solution_I.out', 'w') 1125 for i in range(N): 1126 f.write(str(C[i,1])) 1127 f.write("\n") 1128 f.close 1129 """ 1130 def test_1d_solution_II(self): 1131 print "TEST 1D-SOLUTION II" 1132 1133 L = 2000.0 # Length of channel (m) 1134 N = 100 # Number of compuational cells 1135 cell_len = L/N # Origin = 0.0 1136 1137 points = zeros(N+1,Float) 1138 for i in range(N+1): 1139 points[i] = i*cell_len 1140 1141 domain = Domain(points) 1142 1143 def stage(x): 1144 for i in range(len(x)): 1145 if x[i]<=1000.0: 1146 x[i] = 10.0 1147 else: 1148 x[i] = 0.1 1149 return x 1150 1151 1152 domain.set_quantity('stage', stage) 1153 #L = domain.quantities['stage'].vertex_values 1154 #print "Initial Stage" 1155 #print L 1156 domain.set_boundary({'exterior': Reflective_boundary(domain)}) 1157 1158 import time 1159 t0 = time.time() 1160 yieldstep = 1.0 1161 finaltime = 50.0 1162 1163 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): 1164 a=0 1165 1166 print 'That took %.2f seconds' %(time.time()-t0) 1167 #L = domain.quantities['stage'].vertex_values 1168 #print "Final Stage" 1169 #print L 1170 1171 C = domain.quantities['stage'].vertex_values 1172 #print C 1173 1174 f = file('test_solution_II.out', 'w') 1175 for i in range(N): 1176 f.write(str(C[i,1])) 1177 f.write("\n") 1178 f.close 1179 1180 def test_1d_solution_III(self): 1181 print "TEST 1D-SOLUTION III" 1182 L = 2000.0 # Length of channel (m) 1183 N = 100 # Number of compuational cells 1184 cell_len = L/N # Origin = 0.0 1185 1186 points = zeros(N+1,Float) 1187 for i in range(N+1): 1188 points[i] = i*cell_len 1189 1190 domain = Domain(points) 1191 1192 def stage(x): 1193 for i in range(len(x)): 1194 if x[i]<=1000.0: 1195 x[i] = 10.0 1196 else: 1197 x[i] = 0.0 1198 return x 1199 1200 domain.set_quantity('stage', stage) 1201 #L = domain.quantities['stage'].vertex_values 1202 #print "Initial Stage" 1203 #print L 1204 domain.set_boundary({'exterior': Reflective_boundary(domain)}) 1205 1206 import time 1207 t0 = time.time() 1208 yieldstep = 1.0 1209 finaltime = 30.0 1210 1211 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): 1212 a=0 1213 1214 print 'That took %.2f seconds' %(time.time()-t0) 1215 1216 #L = domain.quantities['stage'].vertex_values 1217 #print "Final Stage" 1218 #print L 1219 1220 C = domain.quantities['stage'].vertex_values 1221 #print C 1222 1223 f = file('test_solution_III.out', 'w') 1224 for i in range(N): 1225 f.write(str(C[i,1])) 1226 f.write("\n") 1227 f.close 1228 """ 1184 1229 #------------------------------------------------------------- 1185 1230 if __name__ == "__main__":
Note: See TracChangeset
for help on using the changeset viewer.