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/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
+