Merge pull request #111 from rjleveque/am574

fvmbook updates

fvmbook/chap21/corner/Makefile

@@ -0,0 +1,63 @@
+
+# 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
+
+# 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 ?= 
+
+# ---------------------------------
+# package sources for this program:
+# ---------------------------------
+
+include $(CLAW)/classic/src/2d/Makefile.classic_2d
+
+# ---------------------------------------
+# package sources specifically to exclude
+# (i.e. if a custom replacement source 
+# under a different name is provided)
+# ---------------------------------------
+
+EXCLUDE_MODULES = \
+
+EXCLUDE_SOURCES = \
+
+
+# ---------------------------------
+# List of custom sources for this program:
+# ---------------------------------
+
+MODULES = \
+
+SOURCES = \
+ qinit.f \
+ setaux.f90 \
+ setprob.f \
+ fdisc.f \
+ $(CLAW)/classic/src/2d/cellave.f \
+ $(CLAW)/riemann/src/rpn2_vc_acoustics.f90 \
+ $(CLAW)/riemann/src/rpt2_vc_acoustics.f90 
+
+#-------------------------------------------------------------------
+# Include Makefile containing standard definitions and make options:
+include $(CLAWMAKE)

fvmbook/chap21/corner/README.rst

@@ -0,0 +1,23 @@
+.. _fvmbook_chap21_corner:
+
+2D acoustics with variable coefficients
+--------------------------------------------------
+
+
+Example [book/chap21/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 2023.
+
+Two materials with a single interface, specified in `fdisc.f`
+(which is called from the library routine `$CLAW/classic/src/2d/cellave.f`
+to compute the fraction of each grid cell lying in the left and right states).
+
+`Figure showing interface <interface.png>`__ created by `makegridfig.py` 
+
+The sound speed and impedance are stored in the aux array, specified in
+`setaux.f`.
+For cells cut by the interface, the arithmetic average of rho and harmonic
+average of the bulk modulus are used to determine these parameters for this
+cell.

fvmbook/chap21/corner/fdisc.f

@@ -0,0 +1,29 @@
+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 or coefficients that
+c # have a discontinuity along some curve. 
+
+c # fdisc should be negative to the "left" of the curve and 
+c # positive to the "right".
+
+c # The cellave routine can then be used to compute the fraction wl of
+c # a grid cell that lies to the "left" of the curve.
+
+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/chap21/corner/makegridfig.py

@@ -0,0 +1,16 @@
+
+# make the figure for the book showing the interface.
+
+from pylab import *
+
+figure(5)
+clf()
+plot([0,0,1],[-1,0,0.55],'g',linewidth=2)
+text(.4,-.6,'right',fontsize=20)
+text(-.6,0,'left',fontsize=20)
+axis('scaled')
+xlim([-1,1])
+ylim([-1,1])
+savefig('interface.png', bbox_inches='tight')
+print("Created interface.png")
+

fvmbook/chap21/corner/qinit.f90

@@ -0,0 +1,33 @@
+subroutine qinit(meqn,mbc,mx,my,xlower,ylower,dx,dy,q,maux,aux)
+
+ ! Set initial conditions for the q array.
+
+ ! Right-going plane wave in acoustics, assuming Z=1
+
+ implicit none
+
+ integer, intent(in) :: meqn,mbc,mx,my,maux
+ real(kind=8), intent(in) :: xlower,ylower,dx,dy
+ real(kind=8), intent(in) :: aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc)
+ real(kind=8), intent(inout) :: q(meqn,1-mbc:mx+mbc,1-mbc:my+mbc)
+
+ ! local variables
+ integer :: i,j
+ real(kind=8) :: xi,yj
+
+
+ do i=1,mx
+ xi = xlower + (i-0.5d0)*dx
+ do j=1,my
+ yj = ylower + (j-0.5d0)*dy
+ if (xi .gt. -0.35d0 .and. xi.lt.-0.2d0) then
+ q(1,i,j) = 1.d0
+ else
+ q(1,i,j) = 0.d0
+ endif
+ q(2,i,j) = q(1,i,j)
+ q(3,i,j) = 0.d0
+ enddo
+ enddo
+
+end subroutine qinit

fvmbook/chap21/corner/setaux.f90

@@ -0,0 +1,54 @@
+subroutine setaux(mbc,mx,my,xlower,ylower,dx,dy,maux,aux)
+
+ ! Called at start of computation before calling qinit, and
+ ! when AMR is used, also called every time a new grid patch is created.
+ ! Use to set auxiliary arrays aux(1,1:maux, 1-mbc:mx+mbc, 1-mbc:my+mbc).
+ ! Note that ghost cell values may need to be set if the aux arrays
+ ! are used by the Riemann solver(s).
+ !
+ ! For variable-coefficient acoustics in a heterogeneous medium
+
+ implicit none
+ integer, intent(in) :: mbc,mx,my,maux
+ real(kind=8), intent(in) :: xlower,ylower,dx,dy
+ real(kind=8), intent(inout) :: aux(1:maux,1-mbc:mx+mbc,1-mbc:my+mbc)
+
+ ! local variables:
+ integer :: i,j
+ real(kind=8) :: rhol,rhor,bulkl,bulkr,wl,wr,xl,yl,rho,bulk
+
+ ! global variables:
+ real(kind=8) :: zl,cl,zr,cr
+ common /comaux/ zl,cl,zr,cr
+
+ ! density and bulk moduli:
+ rhol = zl/cl
+ rhor = zr/cr
+ bulkl = cl**2 * rhol
+ bulkr = cr**2 * rhor
+
+ do j=1-mbc,my+mbc
+ do i=1-mbc,mx+mbc
+ xl = xlower + (i-1.0d0)*dx
+ yl = ylower + (j-1.0d0)*dy
+
+ ! determine what fraction wl of this cell lies to the 
+ ! "left" of the interface:
+ call cellave(xl,yl,dx,dy,wl)
+ wr = 1.d0 - wl
+
+ ! average density:
+ rho = wl*rhol + wr*rhor
+
+ ! harmonic average bulk modulus:
+ bulk = 1.d0 / (wl/bulkl + wr/bulkr)
+
+ ! average sound speed:
+ aux(2,i,j) = dsqrt(bulk/rho)
+
+ ! average impedance:
+ aux(1,i,j) = rho*aux(2,i,j)
+ enddo
+ enddo
+
+end subroutine setaux

fvmbook/chap21/corner/setplot.py

@@ -0,0 +1,131 @@
+
+""" 
+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.
+ 
+""" 
+
+import os
+import numpy as np
+
+
+#--------------------------
+def setplot(plotdata):
+#--------------------------
+
+ """ 
+ Specify what is to be plotted at each frame.
+ Input: plotdata, an instance of pyclaw.plotters.data.ClawPlotData.
+ Output: a modified version of plotdata.
+ 
+ """ 
+
+
+ from clawpack.visclaw import colormaps
+
+ plotdata.clearfigures() # clear any old figures,axes,items data
+
+ def add_interface(current_data):
+ # plot the interface between materials
+ from pylab import plot
+ plot([0,0,1],[-1,0,0.55],'g',linewidth=2)
+
+ def fix_layout(current_data):
+ from pylab import tight_layout
+ tight_layout()
+
+ plotdata.afterframe = fix_layout
+
+ # Figure with 4 subplots
+ plotfigure = plotdata.new_plotfigure(name='solution', figno=0)
+ plotfigure.kwargs = {'figsize':(11,8)}
+
+ # pressure:
+ plotaxes = plotfigure.new_plotaxes()
+ plotaxes.axescmd = 'subplot(2,2,1)'
+ plotaxes.xlimits = [-1,1]
+ plotaxes.ylimits = [-1,1]
+ plotaxes.title = 'pressure'
+ plotaxes.scaled = True
+ plotaxes.afteraxes = add_interface
+
+ plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor')
+ plotitem.plot_var = 0
+ plotitem.pcolor_cmap = colormaps.red_yellow_blue
+ plotitem.pcolor_cmin = -0.2
+ plotitem.pcolor_cmax = 0.2
+ #plotitem.add_colorbar = False
+ plotitem.add_colorbar = True
+ plotitem.colorbar_shrink = 0.7
+
+
+ # Contours of pressure:
+ plotaxes = plotfigure.new_plotaxes()
+ plotaxes.axescmd = 'subplot(2,2,2)'
+ plotaxes.xlimits = [-1,1]
+ plotaxes.ylimits = [-1,1]
+ plotaxes.title = 'pressure'
+ plotaxes.scaled = True
+ plotaxes.afteraxes = add_interface
+
+ plotitem = plotaxes.new_plotitem(plot_type='2d_contour')
+ plotitem.plot_var = 0
+ plotitem.contour_levels = np.linspace(-0.9,0.9,20)
+ plotitem.kwargs = {'linestyles': '-'} # all solid rather than dashed
+
+
+ # u-velocity
+ plotaxes = plotfigure.new_plotaxes()
+ plotaxes.axescmd = 'subplot(2,2,3)'
+ plotaxes.xlimits = [-1,1]
+ plotaxes.ylimits = [-1,1]
+ plotaxes.title = 'u'
+ plotaxes.scaled = True
+ plotaxes.afteraxes = add_interface
+
+ plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor')
+ plotitem.plot_var = 1
+ plotitem.pcolor_cmap = colormaps.red_yellow_blue
+ plotitem.pcolor_cmin = -0.2
+ plotitem.pcolor_cmax = 0.2
+ plotitem.add_colorbar = True
+ plotitem.colorbar_shrink = 0.7
+
+
+ # v-velocity
+ plotaxes = plotfigure.new_plotaxes()
+ plotaxes.axescmd = 'subplot(2,2,4)'
+ plotaxes.xlimits = [-1,1]
+ plotaxes.ylimits = [-1,1]
+ plotaxes.title = 'v'
+ plotaxes.scaled = True
+ plotaxes.afteraxes = add_interface
+
+ plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor')
+ plotitem.plot_var = 2
+ plotitem.pcolor_cmap = colormaps.red_yellow_blue
+ plotitem.pcolor_cmin = -0.2
+ plotitem.pcolor_cmax = 0.2
+ plotitem.add_colorbar = True
+ plotitem.colorbar_shrink = 0.7
+
+
+ # Parameters used only when creating html and/or latex hardcopy
+ # e.g., via pyclaw.plotters.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
+
+