Changeset 5153
 Timestamp:
 Mar 11, 2008, 4:21:13 AM (16 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/demos/netherlands.py
r4713 r5153 6 6 """ 7 7 8 ###################### 9 # Module imports 10 # 8 # 9 # Import necessary modules 10 # 11 12 from anuga.shallow_water import Domain 13 from anuga.shallow_water import Reflective_boundary, Dirichlet_boundary 14 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 15 import os 16 17 #from anuga.visualiser import RealtimeVisualiser 11 18 #import rpdb 12 19 #rpdb.set_active() 13 20 14 from anuga.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\15 Transmissive_boundary, Constant_height, Constant_stage16 21 17 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 18 from Numeric import array 19 from anuga.visualiser import RealtimeVisualiser 22 # 23 # Setup computational domain 24 # 25 26 N = 150 # size = 45000 27 N = 130 # size = 33800 28 N = 600 # Size = 720000 29 N = 100 30 31 points, elements, boundary = rectangular_cross(N, N) 32 domain = Domain(points, elements, boundary, use_inscribed_circle=True) 33 34 domain.check_integrity() 35 36 domain.set_name(os.path.splitext(__file__)[0]) 37 domain.set_timestepping_method('rk3') 38 domain.set_default_order(2) 39 domain.set_store_vertices_uniquely(True) # Store as internally represented 40 domain.tight_slope_limiters = True 41 print domain.statistics() 42 43 44 # Setup order and all the beta's for the limiters (these should become defaults 45 46 #domain.beta_w = 1.0 47 #domain.beta_w_dry = 0.2 48 #domain.beta_uh = 1.0 49 #domain.beta_uh_dry = 0.2 50 #domain.beta_vh = 1.0 51 #domain.beta_vh_dry = 0.2 52 53 #domain.alpha_balance = 100.0 54 55 56 57 # 58 # Setup initial conditions 59 # 20 60 21 61 class Weir: … … 35 75 z = zeros(N, Float) 36 76 for i in range(N): 37 z[i] = x[i]/20 # General slope77 z[i] = x[i]/20 # General slope 38 78 39 # Flattish bit to the left79 # Flattish bit to the left 40 80 if x[i] <= 0.3: 41 81 #z[i] = x[i]/5 … … 43 83 44 84 45 # Weir85 # Weir 46 86 if x[i] > 0.3 and x[i] < 0.4: 47 87 z[i] = x[i]/20+1.2 48 88 49 # Dip89 # Dip 50 90 #if x[i] > 0.6 and x[i] < 0.9: 51 91 # z[i] = x[i]/200.5 #y[i]/5 52 92 53 # Hole in weir93 # Hole in weir 54 94 #if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4: 55 95 if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.4 and y[i] < 0.6: … … 57 97 z[i] = x[i]/20 58 98 59 # Poles99 # Poles 60 100 #if x[i] > 0.65 and x[i] < 0.8 and y[i] > 0.55 and y[i] < 0.65 or\ 61 101 # x[i] > 0.75 and x[i] < 0.9 and y[i] > 0.35 and y[i] < 0.45: … … 67 107 68 108 69 # Wall109 # Wall 70 110 if x[i] > 0.995: 71 111 z[i] = x[i]/20+0.3 … … 74 114 75 115 116 inflow_stage = 0.5 117 manning = 0.0 76 118 77 ###################### 78 # Domain 79 # 119 domain.set_quantity('elevation', Weir(inflow_stage)) 120 domain.set_quantity('friction', manning) 121 domain.set_quantity('stage', expression='elevation + 0.0') 80 122 81 123 82 N = 150 #size = 45000 83 N = 130 #size = 33800 84 N = 600 #Size = 720000 85 N = 100 124 # 125 # Setup boundary conditions 126 # 127 128 Br = Reflective_boundary(domain) 129 Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0]) # Constant inflow 130 domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br}) 86 131 87 132 88 #N = 15 133 # 134 # Evolve system through time 135 # 89 136 90 print 'Creating domain'91 #Create basic mesh92 points, elements, boundary = rectangular_cross(N, N)93 137 94 #Create shallow water domain 95 domain = Domain(points, elements, boundary, use_inscribed_circle=True) 138 if N <= 150: 139 # Initialise realtime visualiser 96 140 97 domain.check_integrity() 98 99 domain.set_timestepping_method('rk3') 100 domain.set_default_order(2) 101 102 #Setup order and all the beta's for the limiters (these should become defaults 103 104 ## domain.beta_w = 1.0 105 ## domain.beta_w_dry = 0.2 106 ## domain.beta_uh = 1.0 107 ## domain.beta_uh_dry = 0.2 108 ## domain.beta_vh = 1.0 109 ## domain.beta_vh_dry = 0.2 110 111 ## domain.alpha_balance = 100.0 112 113 #Output params 114 domain.smooth = False 115 domain.reduction = min #Looks a lot better on top of steep slopes 116 117 print "Number of triangles = ", len(domain) 118 print "Extent = ", domain.get_extent() 119 120 #Set bedslope and friction 121 inflow_stage = 0.5 122 manning = 0.0 123 Z = Weir(inflow_stage) 124 125 if N > 150: 126 domain.checkpoint = False 127 domain.store = True #Store for visualisation purposes 128 domain.format = 'sww' #Native netcdf visualisation format 129 import sys, os 130 #FIXME: This was os.path.splitext but caused weird filenames based on root 131 base = os.path.basename(sys.argv[0]) 132 basename, _ = os.path.splitext(base) 133 domain.set_name(basename) 134 else: 135 domain.checkpoint = False 136 domain.store = False 137 vis = RealtimeVisualiser(domain) 138 vis.render_quantity_height("elevation", dynamic=False) 139 vis.render_quantity_height("stage", dynamic=True) 140 vis.colour_height_quantity('stage', (0.0, 0.0, 0.8)) 141 vis.start() 141 pass 142 #vis = RealtimeVisualiser(domain) 143 #vis.render_quantity_height("elevation", dynamic=False) 144 #vis.render_quantity_height("stage", dynamic=True) 145 #vis.colour_height_quantity('stage', (0.0, 0.0, 0.8)) 146 #vis.start() 142 147 143 148 144 149 145 print 'Field values'146 domain.set_quantity('elevation', Z)147 domain.set_quantity('friction', manning)148 149 150 ######################151 # Boundary conditions152 #153 print 'Boundaries'154 Br = Reflective_boundary(domain)155 Bt = Transmissive_boundary(domain)156 157 #Constant inflow158 Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0])159 160 161 #Set boundary conditions162 domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br})163 164 165 ######################166 #Initial condition167 #168 print 'Initial condition'169 domain.set_quantity('stage', expression='elevation + 0.0')170 171 #Evolve172 150 import time 173 151 t0 = time.time() 174 152 175 153 for t in domain.evolve(yieldstep = 0.005, finaltime = None): 176 domain.write_time() 177 #domain.write_boundary() 178 vis.update() 154 print domain.timestepping_statistics() 179 155 print domain.quantities['stage'].get_values(location='centroids', 180 156 indices=[0]) 157 158 #vis.update() 181 159 #time.sleep(0.1) 182 160 #raw_input('pause>') … … 184 162 #rpdb.set_active() 185 163 #integral_label.text='Integral=%10.5e'%domain.quantities['stage'].get_integral() 186 vis.evolveFinished()164 #vis.evolveFinished() 187 165 188 166 print 'That took %.2f seconds' %(time.time()t0) 189 167 190 191 vis.join() 168 #vis.join()
Note: See TracChangeset
for help on using the changeset viewer.