Merge branch 'master' of git://github.com/clawpack/pyclaw into boxlib…

…-merge Conflicts: examples/acoustics_1d_homogeneous/acoustics_1d.py examples/acoustics_1d_homogeneous/test_acoustics.py examples/advection_1d/test_advection.py
clawpack · Nov 10, 2014 · 414250e · 414250e
2 parents 2fbf81e + 9426db9
commit 414250e
Show file tree

Hide file tree

Showing 62 changed files with 1,962 additions and 674 deletions.
diff --git a/.travis_build_petsc4py.sh b/.travis_build_petsc4py.sh
@@ -2,6 +2,10 @@
 
 # Get <#>
 git clone https://github.com/hashdist/hashdist
+( 
+    cd hashdist
+    git checkout v0.3
+)
 export PATH=`pwd`/hashdist/bin:$PATH
 
 # a reasonable profile for building petsc4py on Travis

diff --git a/CHANGES.md b/CHANGES.md
@@ -1,4 +1,12 @@
-# Changes since 5.1.0
+#5.2.1 release
+
+- **New 3D Euler examples**: Sedov blast, shock tube, and shock-bubble interaction, with IPython notebooks showing how to visualize.  The Sedov example has a test using HDF5 files.
+- **New interface to boundary condition routines**: the routines for setting q and aux both get both q and aux now.  The old interface is supported for backward compatibility until 5.3.  All examples have been updated to use the new interface.
+- New 2D Euler example: flow over a forward-facing step.  It runs, but the solution blows up with the Roe solvers; this is a known issue.
+- In SharpClaw, the total fluctuation solver can now be specified in a manner similar to how the Riemann solver is specified.  The 1D Euler examples demonstrate this.
+- The Riemann solvers in clawpack.riemann now have defined constants for each conserved quantity.  Many of the PyClaw examples have been updated to use those for indexing into the q array.
+
+# 5.2.0 release
 
 - **Linear multistep methods for timestepping in SharpClaw**: 
   Previously, only Runge-Kutta time stepping was included in SharpClaw.
@@ -9,12 +17,16 @@
   The user must provide a routine that computes the eigenvectors.  See 
   `clawpack/pyclaw/examples/euler_1d/`.
 - **Logging control**: it is now easy to modify the logging levels interactively,
-  without modifying the logger files.
+  without modifying the logger files.  All logging configuration is
+  set by default with `pyclaw/log.config`; the file `petclaw/log.config` is
+  no longer used.
 - It is no longer necessary to compile dummy transverse Riemann solvers when using
   an algorithm with no transverse wave propagation.
 - Add functions to compute cell centers of ghost cells.
 - Tests run in parallel on Travis.
-- Writing of HDF5 and netcdf files now works in serial and parallel.
+- All tests now run in under two minutes on most systems.
+- PyWENO-generated code updated to match latest PyWENO release.
+- Writing and reading of HDF5 and netcdf files now works in serial and parallel.
 - Improvements to examples.
 - Miscellaneous bug fixes.
 

diff --git a/README.md b/README.md
@@ -1,3 +1,7 @@
+[![Build Status](https://travis-ci.org/clawpack/pyclaw.svg?branch=master)](https://travis-ci.org/clawpack/pyclaw)
+[![Coverage Status](https://img.shields.io/coveralls/clawpack/pyclaw.svg)](https://coveralls.io/r/clawpack/pyclaw?branch=master)
+
+
 Quick start:
 
     git clone [email protected]:clawpack/clawpack.git

diff --git a/examples/__init__.py b/examples/__init__.py
@@ -12,6 +12,7 @@
 from kpp import kpp
 from psystem_2d import psystem_2d
 from shallow_1d import dam_break
+from shallow_1d import sill
 from shallow_2d import radial_dam_break
 from shallow_sphere import Rossby_wave
 from stegoton_1d import stegoton

diff --git a/examples/acoustics_1d_homogeneous/acoustics_1d.py b/examples/acoustics_1d_homogeneous/acoustics_1d.py
@@ -7,8 +7,8 @@
 
 Solve the (linear) acoustics equations:
 
-.. math:: 
-    p_t + K u_x & = 0 \\ 
+.. math::
+    p_t + K u_x & = 0 \\
     u_t + p_x / \rho & = 0.
 
 Here p is the pressure, u is the velocity, K is the bulk modulus,
@@ -22,8 +22,10 @@
 from numpy import sqrt, exp, cos
 from clawpack import riemann
 
-def setup(use_petsc=False,use_boxlib=False,kernel_language='Fortran',solver_type='classic',
-          outdir='./_output',weno_order=5,time_integrator='SSP104', disable_output=False):
+def setup(use_petsc=False, use_boxlib=False, kernel_language='Fortran', solver_type='classic',
+              outdir='./_output', ptwise=False, weno_order=5,
+              time_integrator='SSP104',
+              disable_output=False):
 
     if use_petsc:
         import clawpack.petclaw as pyclaw
@@ -33,7 +35,11 @@ def setup(use_petsc=False,use_boxlib=False,kernel_language='Fortran',solver_type
         from clawpack import pyclaw
 
     if kernel_language == 'Fortran':
-        riemann_solver = riemann.acoustics_1D
+        if ptwise:
+            riemann_solver = riemann.acoustics_1D_ptwise
+        else:
+            riemann_solver = riemann.acoustics_1D
+
     elif kernel_language=='Python':
         riemann_solver = riemann.acoustics_1D_py.acoustics_1D
 
@@ -63,7 +69,7 @@ def setup(use_petsc=False,use_boxlib=False,kernel_language='Fortran',solver_type
     state.problem_data['bulk']=bulk
     state.problem_data['zz']=sqrt(rho*bulk) # Impedance
     state.problem_data['cc']=sqrt(bulk/rho) # Sound speed
- 
+
     xc=domain.grid.x.centers
     beta=100; gamma=0; x0=0.75
     state.q[0,:] = exp(-beta * (xc-x0)**2) * cos(gamma * (xc - x0))
@@ -86,11 +92,11 @@ def setup(use_petsc=False,use_boxlib=False,kernel_language='Fortran',solver_type
 
 
 def setplot(plotdata):
-    """ 
+    """
     Specify what is to be plotted at each frame.
     Input:  plotdata, an instance of visclaw.data.ClawPlotData.
     Output: a modified version of plotdata.
-    """ 
+    """
     plotdata.clearfigures()  # clear any old figures,axes,items data
 
     # Figure for pressure
@@ -108,7 +114,7 @@ def setplot(plotdata):
     plotitem.plotstyle = '-o'
     plotitem.color = 'b'
     plotitem.kwargs = {'linewidth':2,'markersize':5}
-    
+
     # Set up for axes in this figure:
     plotaxes = plotfigure.new_plotaxes()
     plotaxes.axescmd = 'subplot(212)'
@@ -122,7 +128,7 @@ def setplot(plotdata):
     plotitem.plotstyle = '-'
     plotitem.color = 'b'
     plotitem.kwargs = {'linewidth':3,'markersize':5}
-    
+
     return plotdata
 
 

diff --git a/examples/acoustics_1d_homogeneous/test_acoustics.py b/examples/acoustics_1d_homogeneous/test_acoustics.py
@@ -16,7 +16,7 @@ def acoustics_verify(claw):
             qfinal=claw.frames[claw.num_output_times].state.get_q_global()
 
             # and q_global is only returned on process 0
-            if q0 != None and qfinal != None:
+            if q0 is not None and qfinal is not None:
                 q0 = q0.reshape([-1])
                 qfinal = qfinal.reshape([-1])
                 dx=claw.solution.domain.grid.delta[0]
@@ -31,6 +31,11 @@ def acoustics_verify(claw):
     classic_tests = gen_variants(acoustics_1d.setup, verify_expected(0.00104856594174),
                                  kernel_languages=('Python','Fortran'), solver_type='classic', disable_output=True)
 
+    ptwise_tests = gen_variants(acoustics_1d.setup, verify_expected(0.00104856594174),
+                                 kernel_languages=('Fortran',), ptwise=True,
+                                 solver_type='classic', disable_output=True)
+
+
     sharp_tests_rk   = gen_variants(acoustics_1d.setup, verify_expected(0.000298879563857),
                                  kernel_languages=('Python','Fortran'), solver_type='sharpclaw',
                                  time_integrator='SSP104', disable_output=True)
@@ -44,7 +49,7 @@ def acoustics_verify(claw):
                                  time_integrator='SSP104', weno_order=17, disable_output=True)
 
     from itertools import chain
-    for test in chain(classic_tests, sharp_tests_rk, sharp_tests_lmm, weno_tests):
+    for test in chain(classic_tests, ptwise_tests, sharp_tests_rk, sharp_tests_lmm, weno_tests):
         yield test
 
 

diff --git a/examples/acoustics_2d_homogeneous/acoustics_2d.py b/examples/acoustics_2d_homogeneous/acoustics_2d.py
@@ -17,8 +17,9 @@
 from clawpack import riemann
 import numpy as np
 
-def setup(kernel_language='Fortran',use_petsc=False,outdir='./_output',solver_type='classic',
-        time_integrator='SSP104', disable_output=False):
+def setup(kernel_language='Fortran', use_petsc=False, outdir='./_output', 
+              solver_type='classic', time_integrator='SSP104', ptwise=False,
+              disable_output=False):
     """
     Example python script for solving the 2d acoustics equations.
     """
@@ -27,8 +28,11 @@ def setup(kernel_language='Fortran',use_petsc=False,outdir='./_output',solver_ty
     else:
         from clawpack import pyclaw
 
-    if solver_type=='classic':
-        solver=pyclaw.ClawSolver2D(riemann.acoustics_2D)
+    if solver_type == 'classic':
+        if ptwise:
+            solver = pyclaw.ClawSolver2D(riemann.acoustics_2D_ptwise)
+        else:
+            solver = pyclaw.ClawSolver2D(riemann.acoustics_2D)
         solver.dimensional_split=True
         solver.cfl_max = 0.5
         solver.cfl_desired = 0.45

diff --git a/examples/acoustics_2d_homogeneous/test_2d_acoustics.py b/examples/acoustics_2d_homogeneous/test_2d_acoustics.py
@@ -27,6 +27,9 @@ def verify(claw):
     classic_tests = gen_variants(acoustics_2d.setup, verify_data('verify_classic.txt'),
                                  kernel_languages=('Fortran',), solver_type='classic', disable_output=True)
 
+    ptwise_tests = gen_variants(acoustics_2d.setup, verify_data('verify_classic.txt'),
+                                 kernel_languages=('Fortran',), ptwise=True, solver_type='classic', disable_output=True)
+
     sharp_tests_rk   = gen_variants(acoustics_2d.setup, verify_data('verify_sharpclaw.txt'),
                                  kernel_languages=('Fortran',), solver_type='sharpclaw', 
                                  time_integrator='SSP104', disable_output=True)
@@ -36,5 +39,10 @@ def verify(claw):
                                  time_integrator='SSPMS32', disable_output=True)
 
     from itertools import chain
-    for test in chain(classic_tests, sharp_tests_rk, sharp_tests_lmm):
+    for test in chain(classic_tests, ptwise_tests, sharp_tests_rk, sharp_tests_lmm):
         yield test
+
+
+if __name__=="__main__":
+    import nose
+    nose.main()
diff --git a/examples/acoustics_2d_variable/test_acoustics_2d_variable.py b/examples/acoustics_2d_variable/test_acoustics_2d_variable.py
@@ -13,7 +13,7 @@ def verify_acoustics(controller, solver_type='classic'):
         dx, dy = controller.solution.domain.grid.delta
         test_q = state.get_q_global()
 
-        if test_q != None:
+        if test_q is not None:
             thisdir = os.path.dirname(__file__)
             expected_pressure = np.loadtxt(os.path.join(thisdir,
                                                'pressure_%s.txt' % solver_type))
@@ -48,3 +48,8 @@ def verify_acoustics(controller, solver_type='classic'):
     from itertools import chain
     for test in chain(classic_tests, sharp_tests_rk, sharp_tests_lmm):
         yield test
+
+
+if __name__=="__main__":
+    import nose
+    nose.main()
diff --git a/examples/acoustics_2d_variable/test_acoustics_2d_variable_io.py b/examples/acoustics_2d_variable/test_acoustics_2d_variable_io.py
@@ -71,3 +71,8 @@ def verify_acoustics_io(controller):
             shutil.rmtree(tempdir )
         except OSError as (errno, strerror):
             print ERROR_STR % {'path' : tempdir, 'error': strerror }
+
+
+if __name__=="__main__":
+    import nose
+    nose.main()
diff --git a/examples/acoustics_3d_variable/test_3d_acoustics.py b/examples/acoustics_3d_variable/test_3d_acoustics.py
@@ -69,3 +69,8 @@ def acoustics_verify_heterogeneous(claw):
     for test in chain(classic_homogeneous_tests, classic_heterogeneous_tests, sharp_homogeneous_tests,
                       sharp_heterogeneous_tests):
         yield test
+
+
+if __name__=="__main__":
+    import nose
+    nose.main()
diff --git a/examples/advection_1d/test_advection.py b/examples/advection_1d/test_advection.py
@@ -14,7 +14,7 @@ def advection_verify(claw):
             q0=claw.frames[0].state.get_q_global()
             qfinal=claw.frames[claw.num_output_times].state.get_q_global()
 
-            if q0 != None and qfinal != None:
+            if q0 is not None and qfinal is not None:
                 dx=claw.solution.domain.grid.delta[0]
                 test = dx*np.linalg.norm(qfinal-q0,1)
                 return check_diff(expected, test, reltol=1e-4)
@@ -37,7 +37,7 @@ def advection_verify(claw):
                                  solver_type='sharpclaw',time_integrator='SSPMS32', outdir=None)
 
     weno_tests = gen_variants(advection_1d.setup, verify_expected(7.489618e-06),
-                                 kernel_languages=('Fortran',), solver_type='sharpclaw', 
+                                 kernel_languages=('Fortran',), solver_type='sharpclaw',
                                  time_integrator='SSP104', weno_order=17,
                                  outdir=None)
 
@@ -57,4 +57,3 @@ def advection_verify(claw):
 
     for test in test_1d_advection():
         test()
-
diff --git a/examples/advection_2d_annulus/advection_annulus.py b/examples/advection_2d_annulus/advection_annulus.py
@@ -58,7 +58,7 @@ def qinit(state):
                    + A2*np.exp(-beta2*(np.square(R-r2) + np.square(Theta-theta2)))
 
 
-def ghost_velocities_upper(state,dim,t,auxbc,num_ghost):
+def ghost_velocities_upper(state,dim,t,qbc,auxbc,num_ghost):
     """
     Set the velocities for the ghost cells outside the outer radius of the annulus.
     In the computational domain, these are the cells at the top of the grid.
@@ -77,7 +77,7 @@ def ghost_velocities_upper(state,dim,t,auxbc,num_ghost):
         raise Exception('Custom BC for this boundary is not appropriate!')
 
 
-def ghost_velocities_lower(state,dim,t,auxbc,num_ghost):
+def ghost_velocities_lower(state,dim,t,qbc,auxbc,num_ghost):
     """
     Set the velocities for the ghost cells outside the inner radius of the annulus.
     In the computational domain, these are the cells at the bottom of the grid.

diff --git a/examples/euler_1d/Makefile b/examples/euler_1d/Makefile
@@ -2,17 +2,21 @@ F2PY = f2py
 
 SHARPCLAW = ../../src/pyclaw/sharpclaw
 
-ONE_D_SHARPCLAW_SOURCES = evec.f90 tfluct.f90 $(SHARPCLAW)/flux1.f90 $(SHARPCLAW)/ClawParams.f90 $(SHARPCLAW)/workspace.f90 $(SHARPCLAW)/weno.f90 $(SHARPCLAW)/reconstruct.f90
+ONE_D_SHARPCLAW_SOURCES = evec.f90 $(SHARPCLAW)/flux1.f90 $(SHARPCLAW)/ClawParams.f90 $(SHARPCLAW)/workspace.f90 $(SHARPCLAW)/weno.f90 $(SHARPCLAW)/reconstruct.f90
 
 all: 
-	make sharpclaw1.so
+    make sharpclaw1.so
+    make euler_tfluct.so
 
 sharpclaw1.so: $(ONE_D_SHARPCLAW_SOURCES) 
-	${F2PY} -m sharpclaw1 -c $(ONE_D_SHARPCLAW_SOURCES) 
+    ${F2PY} -m sharpclaw1 -c $(ONE_D_SHARPCLAW_SOURCES) 
+
+euler_tfluct.so: euler_tfluct.f90
+    f2py -m $(basename $(notdir $@)) -c $^
 
 clean:
-	rm -f *.o *.so *.pyc *.log
+    rm -f *.o *.so *.pyc *.log
 
 clobber: clean
-	rm -rf _output/
-	rm -rf _plots/
+    rm -rf _output/
+    rm -rf _plots/
diff --git a/examples/euler_1d/tfluct.f90 → examples/euler_1d/euler_tfluct.f90 b/examples/euler_1d/tfluct.f90 → examples/euler_1d/euler_tfluct.f90
@@ -1,7 +1,6 @@
 ! ============================================================================
-subroutine tfluct(ixyz,maxnx,num_eqn,num_waves,num_ghost,mx,ql,qr,auxl,auxr,adq)
+subroutine tfluct1(maxnx,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,adq)
 ! ============================================================================
-!   
 !   "Internal" Riemann solver for the euler equations in 1D.
 !   The riemann problem is solved by assuming a discontinuity at the
 !   center of the i'th cell.
@@ -11,8 +10,8 @@ subroutine tfluct(ixyz,maxnx,num_eqn,num_waves,num_ghost,mx,ql,qr,auxl,auxr,adq)
 !             auxl contains the auxiliary vector at the left edge of each cell
 !             auxr contains the state vector at the right edge of each cell
 !             maxnx is the number of physical points (without ghost cells)
-!             num_ghost is the number of ghost cells
-!             num_eqn is the number of equations
+!             mbc is the number of ghost cells
+!             meqn is the number of equations
 !             ixyz is the dimension index
 !             mx is the size of the patch for the dimension corresponding
 !               to the value of ixyz
@@ -28,18 +27,23 @@ subroutine tfluct(ixyz,maxnx,num_eqn,num_waves,num_ghost,mx,ql,qr,auxl,auxr,adq)
 !             |_  ( gamma*q3 - 0.5(gamma-1)q2^2/q1 ) * q2/g1  _|
 
 
-    implicit double precision (a-h,o-z)
-    integer,          intent(in) :: maxnx, num_eqn, num_waves, num_ghost, mx, ixyz
-    double precision, intent(in) :: ql(num_eqn,1-num_ghost:maxnx+num_ghost)
-    double precision, intent(in) :: qr(num_eqn,1-num_ghost:maxnx+num_ghost)
-    double precision, intent(out) :: adq(num_eqn,1-num_ghost:maxnx+num_ghost)
+    implicit none
+    integer,          intent(in)  :: maxnx, mx, meqn, mwaves, maux, mbc
+    double precision, intent(in)  :: auxl(maux,1-mbc:maxnx+mbc)
+    double precision, intent(in)  :: auxr(maux,1-mbc:maxnx+mbc)
+    double precision, intent(in)  ::   ql(meqn,1-mbc:maxnx+mbc)
+    double precision, intent(in)  ::   qr(meqn,1-mbc:maxnx+mbc)
+
+    double precision, intent(out) :: adq(meqn,1-mbc:maxnx+mbc)
+
+    integer :: i
+    double precision :: gamma, gamma1
 
 !   local storage
 !   ---------------
-    common /cparam/  gamma1
-    double precision :: gamma
+    common /cparam/  gamma
 
-    gamma = gamma1 + 1.d0
+    gamma1 = gamma - 1.d0
 
     do i = 1,mx
         adq(1,i) = qr(2,i) - ql(2,i)
@@ -50,5 +54,4 @@ subroutine tfluct(ixyz,maxnx,num_eqn,num_waves,num_ghost,mx,ql,qr,auxl,auxr,adq)
     enddo
 
     return
-end subroutine tfluct
-
+end subroutine tfluct1
diff --git a/examples/euler_1d/evec.f90 b/examples/euler_1d/evec.f90
@@ -22,11 +22,13 @@ subroutine evec(maxnx,num_eqn,num_ghost,mx,q,auxl,auxr,evl,evr)
 
 !   local storage
 !   ---------------
-    common /cparam/ gamma1
+    common /cparam/ gamma
     integer :: mx2
 
     mx2 = size(q,2)
 
+    gamma1 = gamma - 1.d0
+
     do i=2,mx2
         ! Compute velocity, speed and enthalpy
         rhsqrtl = dsqrt(q(1,i-1))