Add old history from app repository

1d/README

@@ -0,0 +1,29 @@
+Multi-Layer Shallow Water Examples - 1D
+
+This directory contains a number of examples for solving the two layer shallow
+water equations in one dimension using PyClaw. If you have a valid install of 
+PyClaw you should be able to run the examples from this directory with some
+tweaks such as where you want the output and plots to go (set the DATA_PATH
+environment variable).
+
+Note that many of these scripts were set up to make nice movies and output 
+much more often than is probably necessary. In fact the output takes up a
+significant ammount of the wall clock time (as does the plotting).
+
+Examples
+========
+ - Simple Waves (wave_family.py): Test which uses a pertubation on a steady
+ state in only one wave family. Also includes a jump in bathymetry at 
+ x = 0.5 which causes the bottom layer to become dry.
+ - Oscillatory Wind in a Bath Tub (oscillatory.py): Test demonstrating loss
+ of hyperbolicity. Starts with the water at rest and a wind field applied
+ over the course of the simulation. The boundaries are implemented as wall
+ boundaries in multilayer/bc.py.
+ - Oceanic Shelf (shelf.py): Tests involving oceanic scale flows with a shelf
+ represented with either a jump discontinuity or sloping transition. The
+ initial condition is a small hump of water away from the shelf. The file
+ shelf_contour.py plots contours of the entire simulation.
+ - All Rarefaction Test (rarefaction.py): Example intialized with an all
+ rarefaction solution as the expected Riemann solution.
+ - Dry State Tests (dry_state.py): Example of the entropy violating
+ rarefaction case.

1d/dry_state.py

@@ -0,0 +1,149 @@
+#!/usr/bin/env python
+# encoding: utf-8
+
+r""" Run the suite of tests for the 1d two-layer equations"""
+
+import sys
+
+from clawpack.riemann import layered_shallow_water_1D
+import clawpack.clawutil.runclaw as runclaw
+import clawpack.pyclaw.plot as plot
+
+import multilayer as ml
+
+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 = 'multilayer/dry_state'
+ prefix = 'ml_e%s_m%s_fix' % (eigen_method,num_cells)
+
+ if entropy_fix:
+ prefix = "".join((prefix,"T"))
+ else:
+ prefix = "".join((prefix,"F"))
+ outdir,plotdir,log_path = runclaw.create_output_paths(name,prefix,**kargs)
+
+ # Redirect loggers
+ # This is not working for all cases, see comments in runclaw.py
+ for logger_name in ['pyclaw.io','pyclaw.solution','plot','pyclaw.solver','f2py','data']:
+ runclaw.replace_stream_handlers(logger_name,log_path,log_file_append=False)
+
+ # Load in appropriate PyClaw version
+ if kargs.get('use_petsc',False):
+ import clawpack.petclaw as pyclaw
+ else:
+ import clawpack.pyclaw as pyclaw
+
+
+ # =================
+ # = Create Solver =
+ # =================
+ if kargs.get('solver_type','classic') == 'classic':
+ solver = pyclaw.ClawSolver1D()
+ else:
+ raise NotImplementedError('Classic is currently the only supported solver.')
+
+ # Solver method parameters
+ solver.cfl_desired = 0.9
+ solver.cfl_max = 1.0
+ solver.max_steps = 5000
+ solver.fwave = True
+ solver.kernel_language = 'Fortran'
+ solver.num_waves = 4
+ solver.limiters = 3
+ solver.source_split = 1
+
+ # Boundary conditions
+ solver.bc_lower[0] = 1
+ solver.bc_upper[0] = 1
+ solver.aux_bc_lower[0] = 1
+ solver.aux_bc_upper[0] = 1
+
+ # Set the Riemann solver
+ solver.rp = layered_shallow_water_1D
+
+ # Set the before step function
+ solver.before_step = lambda solver,solution:ml.step.before_step(solver,solution)
+
+ # Use simple friction source term
+ solver.step_source = ml.step.friction_source
+
+ # ============================
+ # = Create Initial Condition =
+ # ============================
+ num_layers = 2
+
+ x = pyclaw.Dimension('x',0.0,1.0,num_cells)
+ domain = pyclaw.Domain([x])
+ state = pyclaw.State(domain,2*num_layers,3+num_layers)
+ state.aux[ml.aux.kappa_index,:] = 0.0
+
+ # Set physics data
+ state.problem_data['g'] = 9.8
+ state.problem_data['manning'] = 0.0
+ state.problem_data['rho_air'] = 1.15e-3
+ state.problem_data['rho'] = [0.95,1.0]
+ state.problem_data['r'] = state.problem_data['rho'][0] / state.problem_data['rho'][1]
+ state.problem_data['one_minus_r'] = 1.0 - state.problem_data['r']
+ state.problem_data['num_layers'] = num_layers
+
+ # Set method parameters, this ensures it gets to the Fortran routines
+ state.problem_data['eigen_method'] = eigen_method
+ state.problem_data['dry_tolerance'] = 1e-3
+ state.problem_data['inundation_method'] = 2
+ state.problem_data['entropy_fix'] = entropy_fix
+
+ solution = pyclaw.Solution(state,domain)
+ solution.t = 0.0
+
+ # Set aux arrays including bathymetry, wind field and linearized depths
+ ml.aux.set_jump_bathymetry(solution.state,0.5,[-1.0,-1.0])
+ ml.aux.set_no_wind(solution.state)
+ ml.aux.set_h_hat(solution.state,0.5,[0.0,-0.5],[0.0,-1.0])
+
+ # Set sea at rest initial condition
+ q_left = [0.5 * state.problem_data['rho'][0],0.0,
+ 0.5 * state.problem_data['rho'][1],0.0]
+ q_right = [1.0 * state.problem_data['rho'][0],0.0,
+ 0.0,0.0]
+ ml.qinit.set_riemann_init_condition(solution.state,0.5,q_left,q_right)
+
+
+ # ================================
+ # = Create simulation controller =
+ # ================================
+ controller = pyclaw.Controller()
+ controller.solution = solution
+ controller.solver = solver
+
+ # Output parameters
+ controller.output_style = 3
+ controller.nstepout = 1
+ controller.num_output_times = 100
+ controller.write_aux_init = True
+ controller.outdir = outdir
+ controller.write_aux = True
+
+
+ # ==================
+ # = Run Simulation =
+ # ==================
+ state = controller.run()
+
+
+ # ============
+ # = Plotting =
+ # ============
+ plot_kargs = {'rho':solution.state.problem_data['rho'],
+ 'dry_tolerance':solution.state.problem_data['dry_tolerance']}
+ 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,
+ **plot_kargs)
+
+if __name__ == "__main__":
+ # Run test case for eigen method = 2 turning on and off entropy fix
+ # dry_state(500,2,True)
+ dry_state(500,2,False,htmlplot=True)

1d/internal_lapping.py

@@ -0,0 +1,158 @@
+#!/usr/bin/env python
+# encoding: utf-8
+
+r""" Run the suite of tests for the 1d two-layer equations"""
+
+import sys
+
+from clawpack.riemann import layered_shallow_water_1D
+import clawpack.clawutil.runclaw as runclaw
+from clawpack.pyclaw.plot import plot
+
+import multilayer as ml
+
+def internal_lapping(num_cells,eigen_method,**kargs):
+ r"""docstring for oscillatory_wind"""
+
+ # Construct output and plot directory paths
+ prefix = 'ml_e%s_n%s' % (eigen_method,num_cells)
+ name = 'multilayer/lapping'
+ outdir,plotdir,log_path = runclaw.create_output_paths(name,prefix,**kargs)
+
+ # Redirect loggers
+ # This is not working for all cases, see comments in runclaw.py
+ for logger_name in ['pyclaw.io','pyclaw.solution','plot','pyclaw.solver','f2py','data']:
+ runclaw.replace_stream_handlers(logger_name,log_path,log_file_append=False)
+
+ # Load in appropriate PyClaw version
+ if kargs.get('use_petsc',False):
+ import clawpack.petclaw as pyclaw
+ else:
+ import clawpack.pyclaw as pyclaw
+
+ # =================
+ # = Create Solver =
+ # =================
+ if kargs.get('solver_type','classic') == 'classic':
+ solver = pyclaw.ClawSolver1D()
+ else:
+ raise NotImplementedError('Classic is currently the only supported solver.')
+
+ # Solver method parameters
+ solver.cfl_desired = 0.9
+ solver.cfl_max = 1.0
+ solver.max_steps = 5000
+ solver.fwave = True
+ solver.kernel_language = 'Fortran'
+ solver.num_waves = 4
+ solver.limiters = 3
+ solver.source_split = 1
+
+ # Boundary conditions
+ solver.bc_lower[0] = 1
+ solver.bc_upper[0] = 1
+ solver.aux_bc_lower[0] = 1
+ solver.aux_bc_upper[0] = 1
+
+ # Set the Riemann solver
+ 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)
+
+ # Use simple friction source term
+ solver.step_source = ml.step.friction_source
+
+
+ # ============================
+ # = Create Initial Condition =
+ # ============================
+ num_layers = 2
+
+ x = pyclaw.Dimension('x',0.0,1.0,num_cells)
+ domain = pyclaw.Domain([x])
+ state = pyclaw.State(domain,2*num_layers,3+num_layers)
+ state.aux[ml.aux.kappa_index,:] = 0.0
+
+ # Set physics data
+ state.problem_data['g'] = 9.8
+ state.problem_data['manning'] = 0.022
+ state.problem_data['rho_air'] = 1.15e-3
+ state.problem_data['rho'] = [0.95,1.0]
+ state.problem_data['r'] = state.problem_data['rho'][0] / state.problem_data['rho'][1]
+ state.problem_data['one_minus_r'] = 1.0 - state.problem_data['r']
+ state.problem_data['num_layers'] = num_layers
+
+ # Set method parameters, this ensures it gets to the Fortran routines
+ state.problem_data['eigen_method'] = eigen_method
+ state.problem_data['dry_tolerance'] = 1e-3
+ state.problem_data['inundation_method'] = 2
+ state.problem_data['entropy_fix'] = False
+
+ solution = pyclaw.Solution(state,domain)
+ solution.t = 0.0
+
+ # Set aux arrays including bathymetry, wind field and linearized depths
+ ml.aux.set_sloped_shelf_bathymetry(solution.state,0.4,0.6,-1.0,-0.2)
+ ml.aux.set_no_wind(solution.state)
+ ml.aux.set_h_hat(solution.state,0.5,[0.0,-0.6],[0.0,-0.6])
+
+ # Set initial condition
+ ml.qinit.set_gaussian_init_condition(solution.state,0.2,0.2,0.01,
+ internal_layer=True)
+
+ # ================================
+ # = Create simulation controller =
+ # ================================
+ controller = pyclaw.Controller()
+ controller.solution = solution
+ controller.solver = solver
+
+ # Output parameters
+ controller.output_style = 1
+ controller.tfinal = 2.0
+ controller.num_output_times = 50
+ controller.write_aux_init = True
+ controller.outdir = outdir
+ controller.write_aux = True
+
+ # ==================
+ # = Run Simulation =
+ # ==================
+ state = controller.run()
+
+ # ============
+ # = Plotting =
+ # ============
+ plot_kargs = {'rho':solution.state.problem_data['rho'],
+ 'dry_tolerance':solution.state.problem_data['dry_tolerance']}
+ 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)
+
+
+if __name__ == "__main__":
+ # Run the test for the 3rd and 4th wave families for each eigen method
+ if len(sys.argv) > 1:
+ eig_methods = []
+ for value in sys.argv[1:]:
+ eig_methods.append(int(value))
+ else:
+ eig_methods = [1,2,3,4]
+
+ # Resolutions for tests
+
+ # Display runs
+ resolution = 128
+ for method in eig_methods:
+ print "Running internal lapping, eigen=%s resolution=%s" % (method,resolution)
+ internal_lapping(resolution,method,iplot=False,htmlplot=True)
+
+ # # Resolutions for tests
+ # resolutions = [64,128,256,512,1024,5000]
+
+ # # Run for comparison runs
+ # for method in eig_methods:
+ # for resolution in resolutions:
+ # print "Running internal lapping, eigen=%s resolution=%s" % (method,resolution)
+ # internal_lapping(resolution,method,iplot=False,htmlplot=True)