Merge pull request clawpack#29 from mandli/multi-layer

Update multilayer SWE application for 5.0

multilayer/1d/dry_state.py

@@ -5,7 +5,7 @@
 
 import sys
 
-import clawpack.riemann as riemann
+from clawpack.riemann import layered_shallow_water_1D
 import clawpack.clawutil.runclaw as runclaw
 import clawpack.pyclaw.plot as plot
 
@@ -15,7 +15,7 @@ def dry_state(num_cells,eigen_method,entropy_fix,**kargs):
     r"""Run and plot a multi-layer dry state problem"""
 
     # Construct output and plot directory paths
-    name = 'dry_state'
+    name = 'multilayer/dry_state'
     prefix = 'ml_e%s_m%s_fix' % (eigen_method,num_cells)
 
     if entropy_fix:
@@ -61,7 +61,7 @@ def dry_state(num_cells,eigen_method,entropy_fix,**kargs):
     solver.aux_bc_upper[0] = 1
 
     # Set the Riemann solver
-    solver.rp = riemann.layered_shallow_water_1D
+    solver.rp = layered_shallow_water_1D
 
     # Set the before step function
     solver.before_step = lambda solver,solution:ml.step.before_step(solver,solution)
@@ -137,7 +137,7 @@ def dry_state(num_cells,eigen_method,entropy_fix,**kargs):
     # ============
     plot_kargs = {'rho':solution.state.problem_data['rho'],
                   'dry_tolerance':solution.state.problem_data['dry_tolerance']}
-    plot.plot(setplot_path="./setplot_drystate.py", outdir=outdir, 
+    plot.plot(setplot="./setplot_drystate.py", outdir=outdir, 
               plotdir=plotdir, htmlplot=kargs.get('htmlplot',False), 
               iplot=kargs.get('iplot',False),
               file_format=controller.output_format,

multilayer/1d/internal_lapping.py

@@ -5,7 +5,7 @@
 
 import sys
 
-import clawpack.riemann as riemann
+from clawpack.riemann import layered_shallow_water_1D
 import clawpack.clawutil.runclaw as runclaw
 from clawpack.pyclaw.plot import plot
 
@@ -16,7 +16,7 @@ def internal_lapping(num_cells,eigen_method,**kargs):
 
     # Construct output and plot directory paths
     prefix = 'ml_e%s_n%s' % (eigen_method,num_cells)
-    name = 'lapping'
+    name = 'multilayer/lapping'
     outdir,plotdir,log_path = runclaw.create_output_paths(name,prefix,**kargs)
 
     # Redirect loggers
@@ -55,7 +55,7 @@ def internal_lapping(num_cells,eigen_method,**kargs):
     solver.aux_bc_upper[0] = 1
 
     # Set the Riemann solver
-    solver.rp = riemann.layered_shallow_water_1D
+    solver.rp = layered_shallow_water_1D
 
     # Set the before step functioning including the wind forcing
     solver.before_step = lambda solver,solution:ml.step.before_step(solver,solution)
@@ -126,7 +126,7 @@ def internal_lapping(num_cells,eigen_method,**kargs):
     # ============
     plot_kargs = {'rho':solution.state.problem_data['rho'],
                   'dry_tolerance':solution.state.problem_data['dry_tolerance']}
-    plot(setplot_path="./setplot_lapping.py",outdir=outdir,plotdir=plotdir,
+    plot(setplot="./setplot_lapping.py",outdir=outdir,plotdir=plotdir,
          htmlplot=kargs.get('htmlplot',False),iplot=kargs.get('iplot',False),
          file_format=controller.output_format,**plot_kargs)
 

multilayer/1d/multilayer/step.py

@@ -11,7 +11,6 @@
 
 import numpy as np
 
-import pyclaw.solution as solution
 from aux import set_no_wind,kappa_index
 
 class NegativeDepthError(Exception):
@@ -107,23 +106,30 @@ def before_step(solver,state,wind_func=set_no_wind,dry_tolerance=1e-3,
 
 def friction_source(solver,state,dt,TOLERANCE=1e-30):
     r""""""
-    manning = state.problem_data['manning']
-    g = state.problem_data['g']
-    rho = state.problem_data['rho']
-    dry_tolerance = state.problem_data['dry_tolerance']
-
-    if manning > TOLERANCE:
-        for i in xrange(state.q.shape[1]):
-            h = state.q[2,i] / rho[1]
-            if h < dry_tolerance:
-                h = state.q[0,i] / rho[0]
-                u = state.q[1,i] / rho[0]
-                layer_index = 0
-            else:
-                u = state.q[2,i] / rho[1]
-                layer_index = 1
-
-            gamma = u * g * manning**2 / h**(4/3)
-            dgamma = 1.0 + dt * gamma
-            hu_index = 2 * (layer_index) + 1
-            state.q[hu_index,i] = state.q[hu_index,i] / dgamma * rho[layer_index]
+
+    if state.problem_data['manning'] != 0.0:
+        if isinstance(solver, ClawSolver):
+            manning = state.problem_data['manning']
+            g = state.problem_data['g']
+            rho = state.problem_data['rho']
+            dry_tolerance = state.problem_data['dry_tolerance']
+
+            if manning > TOLERANCE:
+                for i in xrange(state.q.shape[1]):
+                    h = state.q[2,i] / rho[1]
+                    if h < dry_tolerance:
+                        h = state.q[0,i] / rho[0]
+                        u = state.q[1,i] / rho[0]
+                        layer_index = 0
+                    else:
+                        u = state.q[2,i] / rho[1]
+                        layer_index = 1
+
+                    gamma = u * g * manning**2 / h**(4/3)
+                    dgamma = 1.0 + dt * gamma
+                    hu_index = 2 * (layer_index) + 1
+                    state.q[hu_index,i] = state.q[hu_index,i] / dgamma * rho[layer_index]
+        elif isinstance(solver, SharpClawSolver):
+            pass
+        else:
+            raise ValueError("Solver type %s not supported." % type(solver))

multilayer/1d/oscillatory.py

@@ -3,7 +3,7 @@
 
 r""" Run the suite of tests for the 1d two-layer equations"""
 
-import clawpack.riemann as riemann
+from clawpack.riemann import layered_shallow_water_1D
 import clawpack.clawutil.runclaw as runclaw
 from clawpack.pyclaw.plot import plot
 
@@ -14,7 +14,7 @@ def oscillatory_wind(num_cells,eigen_method,**kargs):
 
     # Construct output and plot directory paths
     prefix = 'ml_e%s_n%s' % (eigen_method,num_cells)
-    name = 'oscillatory_wind'
+    name = 'multilayer/oscillatory_wind'
     outdir,plotdir,log_path = runclaw.create_output_paths(name,prefix,**kargs)
 
     # Redirect loggers
@@ -57,7 +57,7 @@ def oscillatory_wind(num_cells,eigen_method,**kargs):
     solver.aux_bc_upper[0] = 1
 
     # Set the Riemann solver
-    solver.rp = riemann.layered_shallow_water_1D
+    solver.rp = layered_shallow_water_1D
 
     # Set the before step functioning including the wind forcing
     wind_func = lambda state:ml.aux.set_oscillatory_wind(state,
@@ -137,7 +137,7 @@ def oscillatory_wind(num_cells,eigen_method,**kargs):
                   'xupper':solution.state.grid.x.upper,
                   'rho':solution.state.problem_data['rho'],
                   'dry_tolerance':solution.state.problem_data['dry_tolerance']}
-    plot(setplot_path="./setplot_oscillatory.py",outdir=outdir,
+    plot(setplot="./setplot_oscillatory.py",outdir=outdir,
          plotdir=plotdir,htmlplot=kargs.get('htmlplot',False),
          iplot=kargs.get('iplot',False),file_format=controller.output_format,
          **plot_kargs)

multilayer/1d/rarefaction.py

@@ -3,9 +3,7 @@
 
 r""" Run the suite of tests for the 1d two-layer equations"""
 
-import sys
-
-import clawpack.riemann as riemann
+from clawpack.riemann import layered_shallow_water_1D
 import clawpack.clawutil.runclaw as runclaw
 from clawpack.pyclaw.plot import plot
 
@@ -20,7 +18,7 @@ def rarefaction(num_cells,eigen_method,entropy_fix,**kargs):
         prefix = "".join((prefix,"T"))
     else:
         prefix = "".join((prefix,"F"))
-    name = 'all_rare'
+    name = 'multilayer/all_rare'
     outdir,plotdir,log_path = runclaw.create_output_paths(name,prefix,**kargs)
 
     # Redirect loggers
@@ -60,7 +58,7 @@ def rarefaction(num_cells,eigen_method,entropy_fix,**kargs):
     solver.aux_bc_upper[0] = 1
 
     # Set the Riemann solver
-    solver.rp = riemann.layered_shallow_water_1D
+    solver.rp = layered_shallow_water_1D
 
     # Set the before step function
     solver.before_step = lambda solver,solution:ml.step.before_step(solver,solution)
@@ -146,7 +144,7 @@ def rarefaction(num_cells,eigen_method,entropy_fix,**kargs):
     # ============
     plot_kargs = {'rho':solution.state.problem_data['rho'],
                   'dry_tolerance':solution.state.problem_data['dry_tolerance']}
-    plot(setplot_path="./setplot_drystate.py",outdir=outdir,plotdir=plotdir,
+    plot(setplot="./setplot_drystate.py",outdir=outdir,plotdir=plotdir,
          htmlplot=kargs.get('htmlplot',False),iplot=kargs.get('iplot',False),
          file_format=controller.output_format,**plot_kargs)
 

multilayer/1d/setplot_drystate.py

@@ -32,9 +32,9 @@
 # Need to do this after the above
 import matplotlib.pyplot as mpl
 
-import clawpack.pyclaw.solution
+from clawpack.pyclaw.solution import Solution
 
-from multilayer.aux import bathy_index, kappa_index
+from multilayer.aux import bathy_index, kappa_index, wind_index
 import multilayer.plot as plot
 
 # matplotlib.rcParams['figure.figsize'] = [6.0,10.0]
@@ -59,8 +59,7 @@ def jump_afteraxes(current_data):
         mpl.title('Layer Velocities')
 
     # Load bathymetery
-    b = clawpack.pyclaw.solution.Solution(0, path=plotdata.outdir, 
-                                         read_aux=True).state.aux[bathy_index,:]
+    b = Solution(0, path=plotdata.outdir,read_aux=True).state.aux[bathy_index,:]
 
     def bathy(cd):
         return b

multilayer/1d/setplot_lapping.py

@@ -7,7 +7,6 @@
     
 """ 
 
-import os
 import numpy as np
 
 # Plot customization
@@ -52,8 +51,6 @@ def setplot(plotdata,rho,dry_tolerance):
     
     """
 
-    bathy_ref_lines = []
-
     # Load bathymetry
     b = Solution(0,path=plotdata.outdir,read_aux=True).state.aux[bathy_index,:]
 

multilayer/1d/setplot_oscillatory.py

@@ -7,8 +7,6 @@
     
 """ 
 
-import os
-
 import numpy as np
 
 # Plot customization

multilayer/1d/setplot_shelf.py

@@ -7,8 +7,6 @@
     
 """ 
 
-import os
-
 import numpy as np
 
 # Plot customization

multilayer/1d/setplot_wave_family.py

@@ -7,7 +7,6 @@
     
 """ 
 
-import os
 import numpy as np
 
 # Plot customization
@@ -52,8 +51,6 @@ def setplot(plotdata,wave_family,rho,dry_tolerance):
     
     """
 
-    bathy_ref_lines = []
-
     # Load bathymetry
     b = Solution(0,path=plotdata.outdir,read_aux=True).state.aux[bathy_index,:]
 

multilayer/1d/shelf.py

@@ -5,7 +5,7 @@
 
 import sys
 
-import clawpack.riemann as riemann
+from clawpack.riemann import layered_shallow_water_1D
 import clawpack.clawutil.runclaw as runclaw
 from clawpack.pyclaw.plot import plot
 
@@ -16,7 +16,7 @@ def jump_shelf(num_cells,eigen_method,**kargs):
 
     # Construct output and plot directory paths
     prefix = 'ml_e%s_n%s' % (eigen_method,num_cells)
-    name = 'jump_shelf'
+    name = 'multilayer/jump_shelf'
     outdir,plotdir,log_path = runclaw.create_output_paths(name,prefix,**kargs)
 
     # Redirect loggers
@@ -57,7 +57,7 @@ def jump_shelf(num_cells,eigen_method,**kargs):
     solver.aux_bc_upper[0] = 1
 
     # Set the Riemann solver
-    solver.rp = riemann.layered_shallow_water_1D
+    solver.rp = layered_shallow_water_1D
 
     # Set the before step function
     solver.before_step = lambda solver,solution:ml.step.before_step(solver,
@@ -112,11 +112,11 @@ def jump_shelf(num_cells,eigen_method,**kargs):
     controller.solver = solver
 
     # Output parameters
-    # controller.output_style = 1
-    # controller.tfinal = 7200.0
-    # controller.num_output_times = 300
-    controller.output_style = 2
-    controller.out_times = [0.0,720.0,2400.0,4800.0,7200.0]
+    controller.output_style = 1
+    controller.tfinal = 7200.0
+    controller.num_output_times = 300
+    # controller.output_style = 2
+    # controller.out_times = [0.0,720.0,2400.0,4800.0,7200.0]
     controller.write_aux_init = True
     controller.outdir = outdir
     controller.write_aux = True
@@ -134,7 +134,7 @@ def jump_shelf(num_cells,eigen_method,**kargs):
                   "g":solution.state.problem_data['g'],
                   "dry_tolerance":solution.state.problem_data['dry_tolerance'],
                   "bathy_ref_lines":[-30e3]}
-    plot(setplot_path="./setplot_shelf.py",outdir=outdir,plotdir=plotdir,
+    plot(setplot="./setplot_shelf.py",outdir=outdir,plotdir=plotdir,
          htmlplot=kargs.get('htmlplot',False),iplot=kargs.get('iplot',False),
          file_format=controller.output_format,**plot_kargs)
 
@@ -144,7 +144,7 @@ def sloped_shelf(num_cells,eigen_method,**kargs):
 
     # Construct output and plot directory paths
     prefix = 'ml_e%s_n%s' % (eigen_method,num_cells)
-    name = 'sloped_shelf'
+    name = 'multilayer/sloped_shelf'
     outdir,plotdir,log_path = runclaw.create_output_paths(name,prefix,**kargs)
 
     # Redirect loggers
@@ -186,7 +186,7 @@ def sloped_shelf(num_cells,eigen_method,**kargs):
     solver.aux_bc_upper[0] = 1
 
     # Set the Riemann solver
-    solver.rp = riemann.layered_shallow_water_1D
+    solver.rp = layered_shallow_water_1D
 
     # Set the before step function
     solver.before_step = lambda solver,solution:ml.step.before_step(solver,solution)
@@ -264,7 +264,7 @@ def sloped_shelf(num_cells,eigen_method,**kargs):
                   "g":solution.state.problem_data['g'],
                   "dry_tolerance":solution.state.problem_data['dry_tolerance'],
                   "bathy_ref_lines":[x0,x1]}
-    plot(setplot_path="./setplot_shelf.py",outdir=outdir,plotdir=plotdir,
+    plot(setplot="./setplot_shelf.py",outdir=outdir,plotdir=plotdir,
          htmlplot=kargs.get('htmlplot',False),iplot=kargs.get('iplot',False),
          file_format=controller.output_format,**plot_kargs)
 
@@ -280,5 +280,5 @@ def sloped_shelf(num_cells,eigen_method,**kargs):
 
     for method in eig_methods:
         jump_shelf(2000,method,iplot=False,htmlplot=True)
-    # for method in eig_methods:
-    #     sloped_shelf(2000,method,iplot=False,htmlplot=True)
+    for method in eig_methods:
+        sloped_shelf(2000,method,iplot=False,htmlplot=True)