rp1ac.f.html CLAWPACK  
 Source file:   rp1ac.f
 Directory:   /Users/rjl/git/rjleveque/clawpack-4.x/doc/sphinx/example-acoustics-1d
 Converted:   Fri Dec 28 2012 at 15:50:55   using clawcode2html
 This documentation file will not reflect any later changes in the source file. 

 

c
c
c     =====================================================
      subroutine rp1(maxm,meqn,mwaves,mbc,mx,ql,qr,auxl,auxr,
     &                  wave,s,amdq,apdq)
c     =====================================================
c
CLAWPACK Riemann solver for constant coefficient acoustics in 1 space dimension. \[ q_t + A q_x = 0. \] where \[ q(x,t) = \vector{ p(x,t)\\ u(x,t)} \] and the coefficient matrix is \[ A = \begin{matrix} 0 & K\\ 1/\rho & 0 \end{matrix}. \]

The parameters $\rho = $ density and $K =$ bulk modulus are set in setprob.data [.html] and the sound speed $c$ and impedance $Z$ are determined from these in setprob.f [.html].

On input:

  • ql contains the state vector at the left edge of each cell,
  • qr contains the state vector at the right edge of each cell,
  • auxl, auxr are not used in this Riemann solver.

Note that the i'th Riemann problem has left state qr(i-1,:) and right state ql(i,:). From the basic clawpack routine step1, rp1 is called with ql = qr = q.

On output:

  • wave contains the waves,
  • s the speeds,
  • amdq the left-going flux difference $\A^-\Delta Q$,
  • apdq the right-going flux difference $\A^+\Delta Q$

For additional documentation on Riemann solvers rp1, see [claw/doc/rp1]

For details on solution of the Riemann problem for acoustics, see Chapter 3 of FVMHP . 

 
c
c
c
      implicit double precision (a-h,o-z)
c
      dimension wave(1-mbc:maxm+mbc, meqn, mwaves)
      dimension    s(1-mbc:maxm+mbc, mwaves)
      dimension   ql(1-mbc:maxm+mbc, meqn)
      dimension   qr(1-mbc:maxm+mbc, meqn)
      dimension apdq(1-mbc:maxm+mbc, meqn)
      dimension amdq(1-mbc:maxm+mbc, meqn)
c
c     local arrays
c     ------------
      dimension delta(2)
c
c     # density, bulk modulus, and sound speed, and impedence of medium:
c     # (should be set in setprob.f)
      common /cparam/ rho,bulk,cc,zz   
c
c
Split the jump in $Q$ at each interface into waves. First find $\alpha^1$ and $\alpha^2$, the coefficients of the 2 eigenvectors: \[ \delta = \alpha^1 \vector{ -Z \\ 1} + \alpha^2 \vector{ Z \\ 1} \]

Note that the eigendecomposition of $A$ is $A = R \Lambda R^{-1}$, with \[ R = \begin{matrix} -Z & Z \\ 1 & 1 \end{matrix}, \quad \Lambda = \begin{matrix} -c & 0 \\ 0 & c \end{matrix}, \quad R^{-1} = \frac{1}{2Z} \begin{matrix} -1 & Z \\ 1 & Z \end{matrix}, \quad \] 

 
c
      do 20 i = 2-mbc, mx+mbc
         delta(1) = ql(i,1) - qr(i-1,1)
         delta(2) = ql(i,2) - qr(i-1,2)
         a1 = (-delta(1) + zz*delta(2)) / (2.d0*zz)
         a2 =  (delta(1) + zz*delta(2)) / (2.d0*zz)
c
c        # Compute the waves.
c
         wave(i,1,1) = -a1*zz
         wave(i,2,1) = a1
         s(i,1) = -cc
c
         wave(i,1,2) = a2*zz
         wave(i,2,2) = a2
         s(i,2) = cc
c
   20    continue
c
c
Compute the leftgoing and rightgoing fluctuations:

For this problem we always have $s^1 =-c \lt 0$ and $s^2 = c\gt 0$, so \[ \A^-\Delta Q = s^1\W^1, \quad \A^+\Delta Q = s^2\W^2. \] 

 
c
      do 220 m=1,meqn
         do 220 i = 2-mbc, mx+mbc
            amdq(i,m) = s(i,1)*wave(i,m,1)
            apdq(i,m) = s(i,2)*wave(i,m,2)
  220       continue
c
      return
      end