Changeset 8990

Timestamp:

Sep 26, 2013, 11:27:01 AM (11 years ago)

Author:

davies

Message:

First-draft Audusse-type solver in anuga_work/gareth/development/experimental/bal_and

Location:

trunk/anuga_work/development/gareth

Files:

• experimental/bal_and (added)
• experimental/bal_and/plot_utils.py (added)
• experimental/bal_and/swb2_boundary_conditions.py (added)
• experimental/bal_and/swb2_domain.py (added)
• experimental/bal_and/swb2_domain_ext.c (added)
• tests/channel_floodplain/channel_floodplain1.py (modified) (1 diff)
• tests/channel_floodplain/plotme.py (modified) (5 diffs)
• tests/parabolic/parabolic.py (modified) (2 diffs)
• tests/runup/runup.py (modified) (6 diffs)
• tests/runup_sinusoid/runup_sinusoid.py (modified) (1 diff)
• tests/runup_sinusoid/runup_sinusoidplot.py (modified) (1 diff)
• tests/shallow_steep_slope/channel_SU_sparse.py (modified) (1 diff)
• tests/wave/run_wave.py (modified) (3 diffs)

trunk/anuga_work/development/gareth/tests/channel_floodplain/channel_floodplain1.py

@@ -14 +14 @@
 #from anuga import *
 #from balanced_basic import *
-from balanced_dev import *
+#from balanced_dev import *
+from bal_and import *
 #from anuga_tsunami import *
 #from balanced_basic.swb2_domain import *

trunk/anuga_work/development/gareth/tests/channel_floodplain/plotme.py

@@ -1 @@
-import util
+#import util
+from bal_and import plot_utils as util
 from matplotlib import pyplot as pyplot
 #import pylab
 
 # Time-index to plot outputs from
-index=900
+index=150
+#index=900
 
 #p2 = util.get_output('channel_floodplain1_bal_dev.sww', minimum_allowed_height=0.01)
@@ -28 +30 @@
 
 # Numerical results along a central channel 'slice'
-V1 = p.stage[index,v] - p.elev[v]
+V1 = p.stage[index,v] - p.elev[index,v]
 V2 = p.yvel[index,v]
 V3 = p.xvel[index,v]
@@ -90 +92 @@
 
 pyplot.clf()
-analytical_stage = min(p.elev[v1]) + dc_analytical
-analytic_vel = ( (1./300.)*(analytical_stage-p.elev[v1])**(4./3.)*(1./0.03)**2.)**0.5
-analytic_vel = analytic_vel*(analytical_stage>p.elev[v1])
+analytical_stage = min(p.elev[index,v1]) + dc_analytical
+analytic_vel = ( (1./300.)*(analytical_stage-p.elev[index,v1])**(4./3.)*(1./0.03)**2.)**0.5
+analytic_vel = analytic_vel*(analytical_stage>p.elev[index,v1])
 pyplot.figure(figsize=(12.,8.))
 pyplot.plot(p.x[v1], p.yvel[index,v1],'o', label='computed velocity (m/s)')
 pyplot.plot(p.x[v1], analytic_vel,'o', label='analytical velocity (m/s)')
-pyplot.plot(p.x[v1],p.elev[v1],'o', label='bed elevation (m)')
+pyplot.plot(p.x[v1],p.elev[index,v1],'o', label='bed elevation (m)')
 pyplot.plot(p.x[v1],p.stage[index,v1],'o', label='computed stage (m)')
 pyplot.plot(p.x[v1],p.stage[index,v1]*0. + analytical_stage,'o', label='analytical stage (m)')
@@ -109 +111 @@
 
 pyplot.clf()
-analytical_stage = min(p.elev[v1]) + dc_analytical
-analytic_vel = ( (1./300.)*(analytical_stage-p.elev[v1])**(4./3.)*(1./0.03)**2.)**0.5
-analytic_vel = analytic_vel*(analytical_stage>p.elev[v1])
+analytical_stage = min(p.elev[index,v1]) + dc_analytical
+analytic_vel = ( (1./300.)*(analytical_stage-p.elev[index,v1])**(4./3.)*(1./0.03)**2.)**0.5
+analytic_vel = analytic_vel*(analytical_stage>p.elev[index,v1])
 pyplot.figure(figsize=(12.,8.))
 pyplot.plot(p.x[v1], p.yvel[index,v1],'o', label='computed velocity (m/s)')
 pyplot.plot(p.x[v1], analytic_vel,'o', label='analytical velocity (m/s)')
-pyplot.plot(p.x[v1],p.elev[v1],'o', label='bed elevation (m)')
+pyplot.plot(p.x[v1],p.elev[index,v1],'o', label='bed elevation (m)')
 pyplot.plot(p.x[v1],p.stage[index,v1],'o', label='computed stage (m)')
 pyplot.plot(p.x[v1],p.stage[index,v1]*0. + analytical_stage,'o', label='analytical stage (m)')
@@ -127 +129 @@
 
 pyplot.clf()
-analytical_stage = min(p.elev[v1]) + dc_analytical
-analytic_vel = ( (1./300.)*(analytical_stage-p.elev[v1])**(4./3.)*(1./0.03)**2.)**0.5
-analytic_vel = analytic_vel*(analytical_stage>p.elev[v1])
+analytical_stage = min(p.elev[index,v1]) + dc_analytical
+analytic_vel = ( (1./300.)*(analytical_stage-p.elev[index,v1])**(4./3.)*(1./0.03)**2.)**0.5
+analytic_vel = analytic_vel*(analytical_stage>p.elev[index,v1])
 pyplot.figure(figsize=(12.,8.))
 pyplot.plot(p.x[v1], p.yvel[index,v1],'o', label='computed velocity (m/s)')
 pyplot.plot(p.x[v1], analytic_vel,'o', label='analytical velocity (m/s)')
-pyplot.plot(p.x[v1],p.elev[v1],'o', label='bed elevation (m)')
+pyplot.plot(p.x[v1],p.elev[index,v1],'o', label='bed elevation (m)')
 pyplot.plot(p.x[v1],p.stage[index,v1],'o', label='computed stage (m)')
 pyplot.plot(p.x[v1],p.stage[index,v1]*0. + analytical_stage,'o', label='analytical stage (m)')

trunk/anuga_work/development/gareth/tests/parabolic/parabolic.py

@@ -8 +8 @@
 import numpy
 
+from bal_and import *
+
 #from balanced_dev import *
 #from anuga_tsunami import *
 #from balanced_basic import *
 #from anuga.shallow_water_balanced2.swb2_domain import Domain as Domain
-from anuga.shallow_water.shallow_water_domain import Domain as Domain
+#from anuga.shallow_water.shallow_water_domain import Domain as Domain
 #---------
 #Setup computational domain
@@ -21 +23 @@
 domain.set_name('parabola_v2')                         # Output to file runup.sww
 domain.set_datadir('.')                          # Use current folder
-domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
 domain.set_minimum_allowed_height(0.001)
-domain.set_flow_algorithm('tsunami')
+
+#domain.set_flow_algorithm('tsunami')
 # Time stepping
 #domain.set_timestepping_method('euler') # Default

trunk/anuga_work/development/gareth/tests/runup/runup.py

@@ -18 +18 @@
 #path.append('/home/gareth/storage/anuga_clean/anuga_jan12/trunk/anuga_work/development/gareth/balanced_basic')
 #from balanced_basic import *
-from balanced_dev import *
+from bal_and import *
 #from anuga_tsunami import *
 #---------
@@ -28 +28 @@
 domain.set_name('runup_v2')                         # Output to file runup.sww
 domain.set_datadir('.')                          # Use current folder
-domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
+#domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
 #domain.set_store_vertices_uniquely(True)
 #------------------
@@ -34 +34 @@
 #------------------
 
-scale_me=100.0
+scale_me=1.0
 
 def topography(x,y):
@@ -57 +57 @@
 Br=anuga.Reflective_boundary(domain)            # Solid reflective wall
 Bt=anuga.Transmissive_boundary(domain)          # Continue all values of boundary -- not used in this example
-Bd=anuga.Dirichlet_boundary([-0.2*scale_me,0.,0.])       # Constant boundary values -- not used in this example
+Bd=anuga.Dirichlet_boundary([-0.18*scale_me,0.,0.])       # Constant boundary values -- not used in this example
 #Bw=anuga.Time_boundary(domain=domain,
 #       f=lambda t: [(0.0*sin(t*2*pi)-0.1)*exp(-t)-0.1,0.0,0.0]) # Time varying boundary -- get rid of the 0.0 to do a runup.
@@ -64 +64 @@
 # Associate boundary tags with boundary objects
 #----------------------------------------------
-domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom':Br})
+domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom':Br})
 
 #------------------------------
@@ -70 +70 @@
 #------------------------------
 
-for t in domain.evolve(yieldstep=0.2,finaltime=60.0):
+for t in domain.evolve(yieldstep=0.2,finaltime=10.0):
     print domain.timestepping_statistics()
     uh=domain.quantities['xmomentum'].centroid_values
     vh=domain.quantities['ymomentum'].centroid_values
-    depth=domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values
+    depth=domain.quantities['height'].centroid_values#domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values
     depth=depth*(depth>1.0e-06) + 1.0e-06
     vel=((uh/depth)**2 + (vh/depth)**2)**0.5

trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py

@@ -14 +14 @@
 #from swb2_domain import *
 #from balanced_basic import *
-from balanced_dev import *
+#from balanced_dev import *
+from bal_and import *
 #from anuga_tsunami import *
 

trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoidplot.py

@@ -4 +4 @@
 # Import Modules
 #---------------
-import util
+#import util
+from bal_and import plot_utils as util
 from matplotlib import pyplot as pyplot
 import numpy

trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_sparse.py

@@ -21 +21 @@
 from anuga.structures.inlet_operator import Inlet_operator
 from anuga.shallow_water.shallow_water_domain import Domain as Domain
-from balanced_dev import *
+from bal_and import *
 #from balanced_dev import Domain as Domain
 #from anuga_tsunami import *

trunk/anuga_work/development/gareth/tests/wave/run_wave.py

@@ -15 +15 @@
 
 #from balanced_dev import *
-from anuga_tsunami import *
+#from anuga_tsunami import *
+from bal_and import *
 #from balanced_dev import Domain as Domain
 
@@ -89 +90 @@
 amplitude = 1
 wave_length = 300.0
-Bw = anuga.Time_boundary(domain=domain,     # Time dependent boundary
-## Sine wave
-                  f=lambda t: [(-amplitude*sin((1./wave_length)*t*2*pi)), 0.0, 0.0])
+#Bw = anuga.Time_boundary(domain=domain,     # Time dependent boundary
+### Sine wave
+#                  f=lambda t: [(-amplitude*sin((1./wave_length)*t*2*pi)), 0.0, 0.0])
 ## Sawtooth?
 #                   f=lambda t: [(-8.0*(sin((1./180.)*t*2*pi))+(1./2.)*sin((2./180.)*t*2*pi)+(1./3.)*sin((3./180.)*t*2*pi)), 0.0, 0.0])
@@ -129 +130 @@
 #------------------------------------------------------------------------------
 
-for t in domain.evolve(yieldstep = 10.0, finaltime = 60*60.*3.):
+for t in domain.evolve(yieldstep = 10.0, finaltime = 5*10*60.*3.):
     domain.write_time()
     if interactive_visualisation: