Randy LeVeque's proposed refactor

rp2_pwconst_acoustics_module.f90

! Sample rp2 module for an application where the material parameters
! are piecewise constant, with one set of values for x < xdiscont
! and another set for x > xdiscont.   

! Rather than storing values in aux arrays, the rpaux derived type
! could be modified to contain x values a the cell centers.

! rpt2 not yet written and this hasn't been tried, so may be bugs!

module rp2_pwconst_acoustics_module

    ! Module must contain:
    !   expected values of meqn and mwaves
    !   definition of derived type rpaux_type 
    !   subroutines rpn2 and rpt2

    ! use rpaux_type_default_module  ! Could be used for most problems.
                                     ! Then would not need to define here.

    save
    integer, parameter :: meqn=2, mwaves=2

    ! Sound speed and impedance of medium:  (c1,z1) for x < xdiscont
    !                                 and   (c2,z2) for x > xdiscont
    ! These could be set in setprob.f90 via a 
    !    use rp2_pwconst_acoustics_module" 

    real(kind=8) :: c1,z1,c2,z2,xdiscont

    type rpaux_type

        ! Derived data type for other data to be passed into Riemann solver.

        ! Default version would have only the fields time and dx_normal:
    
        real(kind=8) :: time, dx_normal
        ! dx_normal is cell width in direction normal to interfaces
        ! (i.e. dx_normal=dx if ixy==1 or dx_normal=dy if ixy==2)

        ! For this example...

        ! Extended for a problem where x and/or y needed in each cell:
        real(kind=8), allocatable, dimension(:) :: xnormal_l, xnormal_r, &
                      xtransverse_l, xtransverse_r

        ! (xnormal,xtransverse) = (x,y) if ixy==1 or (y,x) if ixy==2.
        ! At cell centers of cells where ql and qr are defined.
        ! Note: the routines step2.f and qad.f would have to be modified
        ! to set these arrays appropriately before each call to rpn2, rpt2.

    end type

contains

!     =====================================================
      subroutine rpn2(ixy,mbc,mx,maux,ql,qr, &
                      auxl,auxr,rpaux,wave,s,amdq,apdq)
!     =====================================================

!     Riemann solver for the acoustics equations in 2d,
!     with piecewise constant coefficients that are computed 
!     using data in rpaux rather than stored in aux arrays.

!     Note that rpaux%time and rpaux%dx_normal are also available.


    implicit none

    integer, intent(in) :: ixy,mbc,mx,maux
    real(kind=8), intent(in) ::   ql(meqn, 1-mbc:mx+mbc)
    real(kind=8), intent(in) ::   qr(meqn, 1-mbc:mx+mbc)
    real(kind=8), intent(in) :: auxl(maux, 1-mbc:mx+mbc)
    real(kind=8), intent(in) :: auxr(maux, 1-mbc:mx+mbc)
    type(rpaux_type), target, intent(in) :: rpaux
    real(kind=8), intent(out) :: wave(meqn, mwaves, 1-mbc:mx+mbc)
    real(kind=8), intent(out) ::    s(mwaves, 1-mbc:mx+mbc)
    real(kind=8), intent(out) :: apdq(meqn, 1-mbc:mx+mbc)
    real(kind=8), intent(out) :: amdq(meqn, 1-mbc:mx+mbc)

!   local variables
!   ---------------
    integer :: mu,mv,i,m
    real(kind=8) :: a1,a2,zim,zi,cim,ci,delta(3)
    real(kind=8), allocatable :: xl(:),xr(:)

    allocate(xl(1-mbc:mx+mbc), xr(1-mbc:mx+mbc))
    ! or could use pointers rather than local arrays and copying!


!     # set mu to point to  the component of the system that corresponds
!     # to velocity in the direction of this slice, mv to the orthogonal
!     # velocity:

    if (ixy == 1) then
        mu = 2
        mv = 3
    else
        mu = 3
        mv = 2
    endif


    ! Set xl, xr to the x-coordinate of the cells holding ql, qr:
    if (ixy==1) then
        xl = rpaux%xnormal_l
        xr = rpaux%xnormal_r
    else
        xl = rpaux%xtransverse_l
        xr = rpaux%xtransverse_r
    endif


!   # split the jump in q at each interface into waves

    do i = 2-mbc, mx+mbc
        ! Set cc and zz based on location in domain:
        if (xr(i-1) < xdiscont) then
            zim = z1
            cim = c1
        else
            zim = z2
            cim = c2
        endif 

        if (xl(i) < xdiscont) then
            zi = z1
            ci = c1
        else
            zi = z2
            ci = c2
        endif 

!       # find a1 and a2, the coefficients of the 2 eigenvectors:
        delta(1) = ql(1,i) - qr(1,i-1)
        delta(2) = ql(mu,i) - qr(mu,i-1)
        a1 = (-delta(1) + zi*delta(2)) / (zim + zi)
        a2 =  (delta(1) + zim*delta(2)) / (zim + zi)
    
!       # Compute the waves:
    
        wave(1,1,i) = -a1*zim
        wave(mu,1,i) = a1
        wave(mv,1,i) = 0.d0
        s(1,i) = -cim
    
        wave(1,2,i) = a2*zi
        wave(mu,2,i) = a2
        wave(mv,2,i) = 0.d0
        s(2,i) = ci
    
     enddo


!     # compute the leftgoing and rightgoing flux differences:
!     # Note s(i,1) < 0   and   s(i,2) > 0.

    forall (m=1:meqn,  i=2-mbc: mx+mbc)
        amdq(m,i) = s(1,i)*wave(m,1,i)
        apdq(m,i) = s(2,i)*wave(m,2,i)
    end forall

    return

    end subroutine rpn2




!   =====================================================
    subroutine rpt2(ixy,imp,mbc,maux,mx,ql,qr,aux1, &
                aux2,aux3,asdq,bmasdq,bpasdq)
!   =====================================================

!   Riemann solver in the transverse direction for the acoustics equations.

!   # Need to add!

    return
    end subroutine rpt2

end module rp2_pwconst_acoustics_module