Changeset 3697
- Timestamp:
- Oct 5, 2006, 12:21:52 PM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/examples/run_sw_rectangle.py
r3696 r3697 1 1 #!/usr/bin/env python 2 ######################################################### 3 # 4 # Example showing improved wet dry limiting. 5 # 6 # 7 # Authors: Steve Roberts 8 # October 9 # 10 # 11 # 12 ######################################################### 2 """ 3 Example showing improved wet dry limiting. 4 This script runs the same simulation twice and stores the outputs in 5 old_limiter.sww and new_limiter.sww respectively. 6 7 Authors: Steve Roberts 8 October 9 """ 13 10 14 11 #----------------------------------------------------------------- 12 # Common structures 13 #----------------------------------------------------------------- 14 import time 15 15 from Numeric import array 16 16 from anuga.shallow_water import Domain 17 18 19 # mesh partition routines 17 from anuga.shallow_water import Reflective_boundary 20 18 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular 21 22 23 M = 3024 points, vertices, boundary = rectangular(M, M, len1 = 1.0, len2 = 1.0)25 domain = Domain(points, vertices, boundary)26 print 'number of triangles = ', domain.number_of_elements27 28 29 #---------------------------30 #Boundaries31 #---------------------------32 from anuga.shallow_water import Reflective_boundary33 34 R = Reflective_boundary(domain)35 36 37 domain.set_boundary( {'left': R, 'right': R, 'bottom': R, 'top': R} )38 domain.check_integrity()39 19 40 20 class Set_Stage: … … 53 33 return self.h0 + self.h*((x>self.x0)&(x<self.x1)&(y>self.y0)&(y<self.y1)) 54 34 35 M = 30 36 points, vertices, boundary = rectangular(M, M, len1 = 1.0, len2 = 1.0) 37 38 yieldstep = 0.002 39 finaltime = 0.06 40 rect = [0.0, 0.0, 1.0, 1.0] 41 42 43 #----------------------------------------------------------------- 44 # Create domain for "old limiter" scenario 45 #----------------------------------------------------------------- 46 domain = Domain(points, vertices, boundary) 47 domain.set_name('old_limiter_second_order') 48 print 'Number of triangles =', len(domain) 49 50 # Turn on the visualisation 51 try: 52 domain.initialise_visualiser() 53 except: 54 pass 55 56 57 #----------------------------------------------------------------- 58 # Boundaries and Initial conditions 59 #----------------------------------------------------------------- 60 R = Reflective_boundary(domain) 61 domain.set_boundary( {'left': R, 'right': R, 'bottom': R, 'top': R} ) 55 62 domain.set_quantity('stage', Set_Stage(0.2, 0.4, 0.25, 0.75, 1.0, 0.00)) 56 63 57 import time 58 t0 = time.time() 59 60 61 # Turn on the visualisation 62 63 rect = [0.0, 0.0, 1.0, 1.0] 64 domain.initialise_visualiser() 65 66 yieldstep = 0.002 67 finaltime = 0.05 68 69 70 #=============================================================================== 71 #Old Limiter 72 #=============================================================================== 73 64 #----------------------------------------------------------------- 65 # Values for old limiter 66 #----------------------------------------------------------------- 74 67 domain.default_order = 2 75 68 domain.beta_w = 0.9 … … 80 73 domain.beta_vh_dry = 0.9 81 74 82 #Check that the boundary value gets propagated to all elements 75 #----------------------------------------------------------------- 76 # Evolve 77 #----------------------------------------------------------------- 78 t0 = time.time() 83 79 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): 84 80 domain.write_time() 85 81 86 87 82 print 'That took %.2f seconds' %(time.time()-t0) 88 83 print 'Note the small timesteps and the irregular flow' 89 raw_input('press return to continue')84 #raw_input('press return to continue') 90 85 91 #===============================================================================92 #New Limiter93 #===============================================================================94 86 95 t0 = time.time() 87 #----------------------------------------------------------------- 88 # Create domain for "new limiter" scenario (2 order) 89 #----------------------------------------------------------------- 90 domain = Domain(points, vertices, boundary) 91 domain.set_name('new_limiter_second_order') 92 print 'Number of triangles =', len(domain) 93 94 # Turn on the visualisation 95 try: 96 domain.initialise_visualiser() 97 except: 98 pass 99 100 #----------------------------------------------------------------- 101 # Boundaries and Initial conditions 102 #----------------------------------------------------------------- 103 R = Reflective_boundary(domain) 104 domain.set_boundary( {'left': R, 'right': R, 'bottom': R, 'top': R} ) 96 105 domain.set_quantity('stage', Set_Stage(0.2, 0.4, 0.25, 0.75, 1.0, 0.00)) 97 domain.set_quantity('xmomentum', 0.0) 98 domain.set_quantity('ymomentum', 0.0) 99 domain.time = 0.0 100 domain.default_order = 2 106 107 #----------------------------------------------------------------- 108 # Values for new limiter 109 #----------------------------------------------------------------- 110 domain.set_default_order(2) 101 111 domain.minimum_allowed_height = 0.001 102 112 domain.beta_w = 1.0 … … 107 117 domain.beta_vh_dry = 0.2 108 118 109 #Check that the boundary value gets propagated to all elements 119 #----------------------------------------------------------------- 120 # Evolve 121 #----------------------------------------------------------------- 122 t0 = time.time() 110 123 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): 111 124 domain.write_time() 112 125 113 114 126 print 'That took %.2f seconds' %(time.time()-t0) 115 127 print 'Note more uniform and large timesteps' 116 raw_input('press return to continue') 117 118 #=============================================================================== 119 #First Order 120 #=============================================================================== 121 122 t0 = time.time() 123 domain.set_quantity('stage', Set_Stage(0.2, 0.4, 0.25, 0.75, 1.0, 0.00)) 124 domain.set_quantity('xmomentum', 0.0) 125 domain.set_quantity('ymomentum', 0.0) 126 domain.time = 0.0 127 128 domain.default_order = 1 128 #raw_input('press return to continue') 129 129 130 130 131 #Check that the boundary value gets propagated to all elements 131 #----------------------------------------------------------------- 132 # Create domain for "new limiter" scenario (1 order) 133 #----------------------------------------------------------------- 134 domain = Domain(points, vertices, boundary) 135 domain.set_name('new_limiter_first_order') 136 print 'Number of triangles =', len(domain) 137 138 # Turn on the visualisation 139 try: 140 domain.initialise_visualiser() 141 except: 142 pass 143 144 #----------------------------------------------------------------- 145 # Boundaries and Initial conditions 146 #----------------------------------------------------------------- 147 R = Reflective_boundary(domain) 148 domain.set_boundary( {'left': R, 'right': R, 'bottom': R, 'top': R} ) 149 domain.set_quantity('stage', Set_Stage(0.2, 0.4, 0.25, 0.75, 1.0, 0.00)) 150 151 #----------------------------------------------------------------- 152 # Values for new limiter first order 153 #----------------------------------------------------------------- 154 domain.set_default_order(1) 155 domain.minimum_allowed_height = 0.001 156 domain.beta_w = 1.0 157 domain.beta_w_dry = 0.2 158 domain.beta_uh = 1.0 159 domain.beta_uh_dry = 0.2 160 domain.beta_vh = 1.0 161 domain.beta_vh_dry = 0.2 162 163 #----------------------------------------------------------------- 164 # Evolve 165 #----------------------------------------------------------------- 166 t0 = time.time() 132 167 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): 133 168 domain.write_time() 134 169 135 136 170 print 'That took %.2f seconds' %(time.time()-t0) 137 171 172
Note: See TracChangeset
for help on using the changeset viewer.