Changeset 3697
 Timestamp:
 Oct 5, 2006, 12:21:52 PM (16 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.