| rp1acv.f.html | mathcode2html |
| Source file: rp1acv.f | |
| Directory: /nfs/aesop01/hw00/d35/rjl/mathcode2html | |
| Converted: Sat Oct 22 2011 at 13:35:54 | |
| This documentation file will not reflect any later changes in the source file. |
$\phantom{******** If you see this on the webpage then the
browser could not locate *********}$
$\phantom{******** the jsMath file load.js *********}$
$\newcommand{\vector}[1]{\left[\begin{array}{c} #1 \end{array}\right]}$ $\newenvironment{matrix}{\left[\begin{array}{cccccccccc}} {\end{array}\right]}$
c
c
c =====================================================
subroutine rp1(maxm,meqn,mwaves,mbc,mx,ql,qr,auxl,auxr,
& wave,s,amdq,apdq)
c =====================================================
c
|
|
CLAWPACK Riemann solver
in 1 space dimension.
Solve Riemann problems for the 1D acoustics equations with variable coefficients (heterogeneous media) \[ q_t + A(x) q_x = 0 \] where \[ q(x,t) = \vector{ p(x,t)\\ u(x,t)} \] and the coefficient matrix is \[ A(x) = \begin{matrix} 0 & K(x)\\ 1/\rho(x) & 0 \end{matrix}. \] Here $p$ is the pressure perturbation and $u$ the velocity. The parameters $K$ (bulk modulus) and $\rho$ (density) are used in the setaux routine to set the impedance and sound speed in each cell. On input:
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 and auxl = auxr = aux. On output:
For additional documentation on Riemann solvers rp1, see www.clawpack.org/rp1.html For details on solution of the Riemann problem for variable coefficient acoustics, see Chapter 9 of FVMHP .
|
c
c
implicit double precision (a-h,o-z)
c
dimension auxl(1-mbc:maxm+mbc, 2)
dimension auxr(1-mbc:maxm+mbc, 2)
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)
common /comlim/ mylim,mrplim(2)
c
c
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.
For the left-going wave we use the eigenvector for a 1-wave propagating with speed $-c_{i-1}$ into the cell on the left, which is \[ r^1_{i-1/2} = \vector{ -Z_{i-1} \\ 1} \] For the right-going wave we use the eigenvector for a 2-wave propagating with speed $+c_i$ into the cell on the right, which is \[ r^2_{i-1/2} = \vector{ Z_i \\ 1} \] Note that this corresponds to solving the Riemann problem \[ q_t + A_{i-1/2} q_x = 0 \] where the eigendecomposition of $A_{i-1/2}$ is $A_{i-1/2} = R_{i-1/2} \Lambda_{i-1/2} R^{-1}_{i-1/2}$, with \[ R_{i-1/2} = \begin{matrix} -Z_{i-1} & Z_i \\ 1 & 1 \end{matrix}, \quad \Lambda_{i-1/2} = \begin{matrix} -c_{i-1} & 0 \\ 0 & c_i \end{matrix}, \] \[ R^{-1}_{i-1/2} = \frac{1}{Z_{i-1}+Z_i} \begin{matrix} -1 & Z_i \\ 1 & Z_{i-1} \end{matrix}, \quad \] |
c
do 20 i = 2-mbc, mx+mbc
c # loop over all cells, including ghost cells
c # impedances:
zi = auxl(i,1)
zim = auxl(i-1,1)
|
|
Set $\delta = Q_i - Q_{i-1}$ |
delta1 = ql(i,1) - qr(i-1,1)
delta2 = ql(i,2) - qr(i-1,2)
|
|
Set $\alpha = R_{i-1/2}^{-1} \delta$ |
a1 = (-delta1 + zi*delta2) / (zim + zi)
a2 = (delta1 + zim*delta2) / (zim + zi)
c
|
|
Compute the waves and speeds: \[ {\cal W}^1 = \alpha^1 r_{i-1/2}^1,~~~~ {\cal W}^2 = \alpha^2 r_{i-1/2}^2, \] \[ s^1 = -c_{i-1},~~~~ s^2 = +c_i. \] |
c
wave(i,1,1) = -a1*zim
wave(i,2,1) = a1
s(i,1) = -auxl(i-1,2)
c
wave(i,1,2) = a2*zi
wave(i,2,2) = a2
s(i,2) = auxl(i,2)
c
20 continue
c
c
|
|
Compute the leftgoing and rightgoing fluctuations:
For this problem we always have $s^1 <0$ and $s^2 > 0$, so \[ {\cal A}^-\Delta Q = s^1{\cal W}^1, \quad {\cal A}^+\Delta Q = s^2{\cal 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