Merge pull request clawpack#43 from rjleveque/chap22

converted fvmbook/chap22/corner

fvmbook/chap22/corner/Makefile

@@ -0,0 +1,65 @@
+
+# Makefile for Clawpack code in this directory.
+# This version only sets the local files and frequently changed
+# options, and then includes the standard makefile pointed to by CLAWMAKE.
+CLAWMAKE = $(CLAW)/clawutil/src/Makefile.common
+
+# See the above file for details and a list of make options, or type
+#   make .help
+# at the unix prompt.
+
+
+# Adjust these variables if desired:
+# ----------------------------------
+
+CLAW_PKG = classic                  # Clawpack package to use
+EXE = xclaw                         # Executable to create
+SETRUN_FILE = setrun.py             # File containing function to make data
+OUTDIR = _output                    # Directory for output
+SETPLOT_FILE = setplot.py           # File containing function to set plots
+PLOTDIR = _plots                    # Directory for plots
+
+OVERWRITE ?= True                   # False ==> make a copy of OUTDIR first
+RESTART ?= False                    # Should = clawdata.restart in setrun
+
+# Environment variable FC should be set to fortran compiler, e.g. gfortran
+
+# Compiler flags can be specified here or set as an environment variable
+FFLAGS ?=  
+
+# ---------------------------------
+# List of sources for this program:
+# ---------------------------------
+
+LIB = $(CLAW)/classic/src/2d
+
+MODULES = \
+
+SOURCES = \
+  qinit.f \
+  setprob.f \
+  fdisc.f \
+  setaux.f \
+  $(CLAW)/classic/src/2d/cellave.f \
+  $(CLAW)/riemann/src/rpn2_vc_elasticity.f90 \
+  $(CLAW)/riemann/src/rpt2_vc_elasticity.f90 \
+  $(CLAW)/classic/src/2d/driver.f90 \
+  $(CLAW)/classic/src/2d/claw2ez.f \
+  $(CLAW)/classic/src/2d/claw2.f \
+  $(CLAW)/classic/src/2d/bc2.f \
+  $(CLAW)/classic/src/2d/b4step2.f90 \
+  $(CLAW)/classic/src/2d/step2.f90 \
+  $(CLAW)/classic/src/2d/step2ds.f90 \
+  $(CLAW)/classic/src/2d/dimsp2.f \
+  $(CLAW)/classic/src/2d/flux2.f90 \
+  $(CLAW)/classic/src/2d/copyq2.f \
+  $(CLAW)/classic/src/2d/inlinelimiter.f90 \
+  $(CLAW)/classic/src/2d/src2.f90 \
+  $(CLAW)/classic/src/2d/out2.f \
+  $(CLAW)/classic/src/2d/restart2.f \
+  $(CLAW)/classic/src/2d/opendatafile.f
+
+#-------------------------------------------------------------------
+# Include Makefile containing standard definitions and make options:
+include $(CLAWMAKE)
+

fvmbook/chap22/corner/README.rst

@@ -0,0 +1,24 @@
+
+.. _fvmbook_chap22/corner:
+
+Elasticity equations in a region consisting of two materials with a corner
+--------------------------------------------------------------------------
+
+
+Example [book/chap22/corner] to accompany the book 
+`Finite Volume Methods for Hyperbolic Problems
+<http://www.clawpack.org/book.html>`_
+by R. J. LeVeque.
+
+Converted to Clawpack 5.0 form in 2014.
+
+
+**Note:**
+
+ - `mx=400, my=400` was used for the plots in the book.
+ - The elastic parameters are defined in `setprob.data`
+ - and the inclusion geometry is specified by `fdisc.f`
+ - Elastic parameters in each cell are set in `setaux.f`
+ - The wave is generated by specifying a right-going P-wave in `qinit.f`
+
+

fvmbook/chap22/corner/fdisc.f

@@ -0,0 +1,24 @@
+c
+c
+c
+c     =================================================
+      function fdisc(x,y)
+c     =================================================
+
+      implicit double precision (a-h,o-z)
+
+c
+c     # for computing cell averages for initial data that has a
+c     # discontinuity along some curve.  fdisc should be negative to the 
+c     # left of the curve and positive to the right
+
+c     # half wedge
+      if (x .gt. 0.d0 .and. (y.lt.(0.55d0*x))) then
+          fdisc = 1.d0
+        else
+          fdisc = -1.d0
+        endif
+c
+      return
+      end
+

fvmbook/chap22/corner/qinit.f

@@ -0,0 +1,38 @@
+
+c
+c
+c
+c     =====================================================
+       subroutine qinit(meqn,mbc,mx,my,xlower,ylower,
+     &                   dx,dy,q,maux,aux)
+c     =====================================================
+c
+c     # Set initial conditions for q.
+c
+       implicit double precision (a-h,o-z)
+       dimension q(meqn,1-mbc:mx+mbc,1-mbc:my+mbc)
+       dimension aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc)
+       common /comaux/ rho1,amu1,alam1,rho2,amu2,alam2
+c
+       cp2 = dsqrt((alam2+2.d0*amu2)/rho2)
+       do 20 i=1,mx
+          xi = xlower + (i-0.5d0)*dx
+          do 20 j=1,my
+             yj = ylower + (j-0.5d0)*dy
+             q(1,i,j) = 0.d0
+             q(2,i,j) = 0.d0
+             q(3,i,j) = 0.d0
+             q(4,i,j) = 0.d0
+             q(5,i,j) = 0.d0
+     	     if (xi .gt. -0.35d0 .and. xi.lt.-0.2d0) then
+c                # set to be an eigenvector corresponding to right-going
+c                # P-wave:
+    	         q(1,i,j) = aux(2,i,j) + 2.d0*aux(3,i,j)  
+    	         q(2,i,j) = aux(2,i,j)
+    	         q(4,i,j) = -aux(4,i,j)
+        		 endif
+  20         continue
+
+       return
+       end
+

fvmbook/chap22/corner/setaux.f

@@ -0,0 +1,42 @@
+c     =========================================================
+      subroutine setaux(mbc,mx,my,xlower,ylower,dx,dy,maux,aux)
+c     =========================================================
+c
+c     # set auxiliary arrays 
+c     # variable coefficient acoustics
+c     #  aux(1,i,j) = density rho in (i,j) cell
+c     #  aux(2,i,j) = lambda in (i,j) cell
+c     #  aux(3,i,j) = mu in (i,j) cell
+c     #  aux(4,i,j) = cp in (i,j) cell
+c     #  aux(5,i,j) = cs in (i,j) cell
+c
+c     # Piecewise constant medium
+c     # Material parameters are set in setprob.f
+
+c
+c     
+      implicit double precision (a-h,o-z)
+      dimension aux(5,1-mbc:mx+mbc,1-mbc:my+mbc)
+      common /comaux/ rho1,amu1,alam1,rho2,amu2,alam2
+c
+
+      do 30 j=1-mbc,my+mbc
+       do 20 i=1-mbc,mx+mbc
+          xl = xlower + (i-1.0d0)*dx
+          yl = ylower + (j-1.0d0)*dy
+          call cellave(xl,yl,dx,dy,w1)
+          w2 = 1.d0 - w1
+
+          aux(1,i,j) = w1*rho1 + w2*rho2
+          aux(2,i,j) = w1*alam1 + w2*alam2
+          aux(3,i,j) = w1*amu1 + w2*amu2
+          bulk       = aux(2,i,j) + 2.d0*aux(3,i,j)
+          aux(4,i,j) = dsqrt(bulk/aux(1,i,j))
+          aux(5,i,j) = dsqrt(aux(3,i,j)/aux(1,i,j))
+
+   20     continue
+   30    continue
+
+       return
+       end
+

fvmbook/chap22/corner/setplot.py

@@ -0,0 +1,125 @@
+
+""" 
+Set up the plot figures, axes, and items to be done for each frame.
+
+This module is imported by the plotting routines and then the
+function setplot is called to set the plot parameters.
+    
+""" 
+
+#--------------------------
+def setplot(plotdata):
+#--------------------------
+
+    """ 
+    Specify what is to be plotted at each frame.
+    Input:  plotdata, an instance of clawpack.visclaw.data.ClawPlotData.
+    Output: a modified version of plotdata.
+    
+    """ 
+
+
+    from clawpack.visclaw import colormaps
+    from numpy import linspace
+
+    plotdata.clearfigures()  # clear any old figures,axes,items data
+
+    def plot_corner(current_data):
+        from pylab import plot
+        plot([.0,.0],[-1,0],'r',linewidth=2)
+        plot([.0,1],[0,.55],'r',linewidth=2)
+
+    def sigmatr(current_data):
+        # trace of sigma
+        q = current_data.q
+        return q[0,:,:] + q[1,:,:]
+
+
+
+    # Figure for trace of sigma
+    plotfigure = plotdata.new_plotfigure(name='trace(sigma)', figno=0)
+
+    # Set up for axes in this figure:
+    plotaxes = plotfigure.new_plotaxes()
+    plotaxes.xlimits = 'auto'
+    plotaxes.ylimits = 'auto'
+    plotaxes.title = 'trace(sigma)'
+    plotaxes.scaled = True
+    plotaxes.afteraxes = plot_corner
+
+    # Set up for item on these axes:
+    plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor')
+    plotitem.plot_var = sigmatr
+    plotitem.pcolor_cmap = colormaps.red_yellow_blue
+    plotitem.pcolor_cmin = -1.
+    plotitem.pcolor_cmax = 1.
+    plotitem.add_colorbar = True
+    plotitem.show = True       # show on plot?
+
+
+
+    # Figure for shear stress
+    plotfigure = plotdata.new_plotfigure(name='shear', figno=1)
+
+    # Set up for axes in this figure:
+    plotaxes = plotfigure.new_plotaxes()
+    plotaxes.xlimits = 'auto'
+    plotaxes.ylimits = 'auto'
+    plotaxes.title = 'shear stress'
+    plotaxes.scaled = True
+    plotaxes.afteraxes = plot_corner
+
+    # Set up for item on these axes:
+    plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor')
+    plotitem.plot_var = 2  # sigma_12
+    plotitem.pcolor_cmap = colormaps.red_yellow_blue
+    plotitem.pcolor_cmin = -0.2
+    plotitem.pcolor_cmax = 0.2
+    plotitem.add_colorbar = True
+    plotitem.show = True       # show on plot?
+
+
+    # Figure for contours 
+    plotfigure = plotdata.new_plotfigure(name='contours', figno=2)
+
+    # Set up for axes in this figure:
+    plotaxes = plotfigure.new_plotaxes()
+    plotaxes.xlimits = 'auto'
+    plotaxes.ylimits = 'auto'
+    plotaxes.title = 'pressure(black) and shear(green)'
+    plotaxes.scaled = True
+    plotaxes.afteraxes = plot_corner
+
+    # Set up for item on these axes:
+    plotitem = plotaxes.new_plotitem(plot_type='2d_contour')
+    plotitem.plot_var = sigmatr
+    plotitem.contour_levels = linspace(-2,8,50)
+    plotitem.contour_colors = 'k'
+    plotitem.show = True       # show on plot?
+
+    # Set up for item on these axes:
+    plotitem = plotaxes.new_plotitem(plot_type='2d_contour')
+    plotitem.plot_var = 2  # sigma_12
+    plotitem.contour_levels = linspace(-0.4,0.4,30)
+    plotitem.contour_colors = 'g'
+    plotitem.show = True       # show on plot?
+
+
+
+    # Parameters used only when creating html and/or latex hardcopy
+    # e.g., via clawpack.visclaw.frametools.printframes:
+
+    plotdata.printfigs = True                # print figures
+    plotdata.print_format = 'png'            # file format
+    plotdata.print_framenos = 'all'          # list of frames to print
+    plotdata.print_fignos = 'all'            # list of figures to print
+    plotdata.html = True                     # create html files of plots?
+    plotdata.html_homelink = '../README.html'   # pointer for top of index
+    plotdata.latex = True                    # create latex file of plots?
+    plotdata.latex_figsperline = 2           # layout of plots
+    plotdata.latex_framesperline = 1         # layout of plots
+    plotdata.latex_makepdf = False           # also run pdflatex?
+
+    return plotdata
+
+

fvmbook/chap22/corner/setprob.f

@@ -0,0 +1,31 @@
+      subroutine setprob
+      implicit double precision (a-h,o-z)
+      character*12 fname
+      common /comaux/ rho1,amu1,alam1,rho2,amu2,alam2
+
+c
+c     # Set the material parameters for the elasticity equations
+c
+c
+      iunit = 7
+      fname = 'setprob.data'
+c     # open the unit with new routine from Clawpack 4.4 to skip over
+c     # comment lines starting with #:
+      call opendatafile(iunit, fname)
+
+
+c
+c     # Piecewise constant medium 
+c     # Material parameters
+
+      read(7,*) rho1
+      read(7,*) alam1
+      read(7,*) amu1
+
+      read(7,*) rho2
+      read(7,*) alam2
+      read(7,*) amu2
+
+      return
+      end
+