geoclaw_riemann_utils_vec.f90

    ! This code has been adapted for vectorization. 
! For more details on how to use it, see the Makefile of the chile2010 example.

!-----------------------------------------------------------------------
      subroutine riemann_aug_JCP(maxiter,meqn,mwaves,hL,hR,huL,huR, &
        hvL,hvR,bL,bR,uL,uR,vL,vR,phiL,phiR,pL,pR,sE1,sE2,drytol,g,rho, &
        sw1,sw2,sw3,fw11,fw12,fw13,fw21,fw22,fw23,fw31,fw32,fw33)
      !dir$ attributes forceinline :: riemann_aug_JCP
      !dir$ attributes vector: uniform(maxiter,meqn,mwaves,drytol,g,rho) :: riemann_aug_JCP

      ! solve shallow water equations given single left and right states
      ! This solver is described in J. Comput. Phys. (6): 3089-3113, March 2008
      ! Augmented Riemann Solvers for the Shallow Equations with Steady States and Inundation

      ! To use the original solver call with maxiter=1.

      ! This solver allows iteration when maxiter > 1. The iteration seems to help with
      ! instabilities that arise (with any solver) as flow becomes transcritical over variable topo
      ! due to loss of hyperbolicity.

      implicit none

      !input
      integer meqn,mwaves,maxiter
      
      !  Due to what seems to be a bug with the Intel Compilers (versions 16, 17 and 18 at least),
      !  these arrays declared as PRIVATE in a !$OMP SIMD loop (in rpn2_geoclaw_vec.f90) had 
      !  to be decomposed into scalars.
      !  This may not be necessary for future versions of the Intel Compilers.
      
      !real(kind=8) fw(meqn,mwaves)
      !real(kind=8) sw(mwaves)
      real(kind=8) sw1,sw2,sw3,fw11,fw12,fw13,fw21,fw22,fw23,fw31,fw32,fw33
      real(kind=8) hL,hR,huL,huR,bL,bR,uL,uR,phiL,phiR,sE1,sE2
      real(kind=8) hvL,hvR,vL,vR,pL,pR
      real(kind=8) drytol,g,rho


      !local
      integer m,mw,k,iter
      real(kind=8) A(3,3)
      real(kind=8) r(3,3)
      real(kind=8) lambda(3)
      real(kind=8) del(3)
      real(kind=8) beta(3)

      real(kind=8) delh,delhu,delphi,delb,delnorm
      real(kind=8) rare1st,rare2st,sdelta,raremin,raremax
      real(kind=8) criticaltol,convergencetol,raretol
      real(kind=8) criticaltol_2, hustar_interface
      real(kind=8) s1s2bar,s1s2tilde,hbar,hLstar,hRstar,hustar
      real(kind=8) huRstar,huLstar,uRstar,uLstar,hstarHLL
      real(kind=8) deldelh,deldelphi,delP
      real(kind=8) s1m,s2m,hm
      real(kind=8) det1,det2,det3,determinant

      logical rare1,rare2,rarecorrector,rarecorrectortest,sonic
      
      ! used to pre-compute square roots
      real(kind=8) sqrt_ghL, sqrt_ghR, sqrt_ghm

      !determine del vectors
      delh = hR-hL
      delhu = huR-huL
      delphi = phiR-phiL
      delb = bR-bL
      delP = pR - pL
      delnorm = delh**2 + delphi**2

      !DIR$ FORCEINLINE
      call riemanntype(hL,hR,uL,uR,hm,s1m,s2m,rare1,rare2, &
                                               1,drytol,g)

      ! pre-compute square roots:
      sqrt_ghL = sqrt(g*hL)
      sqrt_ghR = sqrt(g*hR)
      sqrt_ghm = sqrt(g*hm)

      lambda(1)= min(sE1,s2m) !Modified Einfeldt speed
      lambda(3)= max(sE2,s1m) !Modified Eindfeldt speed
      sE1=lambda(1)
      sE2=lambda(3)
      lambda(2) = 0.d0  ! ### Fix to avoid uninitialized value in loop on mw -- Correct?? ###

      
      hstarHLL = max((huL-huR+sE2*hR-sE1*hL)/(sE2-sE1),0.d0) ! middle state in an HLL solve

      !determine the middle entropy corrector wave------------------------
      rarecorrectortest=.false.
      rarecorrector=.false.
      if (rarecorrectortest) then
         sdelta=lambda(3)-lambda(1)
         raremin = 0.5d0
         raremax = 0.9d0
         if (rare1.and.sE1*s1m.lt.0.d0) raremin=0.2d0
         if (rare2.and.sE2*s2m.lt.0.d0) raremin=0.2d0
         if (rare1.or.rare2) then
            !see which rarefaction is larger
            rare1st=3.d0*(sqrt_ghL-sqrt_ghm)
            rare2st=3.d0*(sqrt_ghR-sqrt_ghm)
            if (max(rare1st,rare2st).gt.raremin*sdelta.and. &
              max(rare1st,rare2st).lt.raremax*sdelta) then
                  rarecorrector=.true.
               if (rare1st.gt.rare2st) then
                  lambda(2)=s1m
               elseif (rare2st.gt.rare1st) then
                  lambda(2)=s2m
               else
                  lambda(2)=0.5d0*(s1m+s2m)
               endif
            endif
         endif
         if (hstarHLL.lt.min(hL,hR)/5.d0) rarecorrector=.false.
      endif

!     ## Is this correct 2-wave when rarecorrector == .true. ??
      do mw=1,mwaves
         r(1,mw)=1.d0
         r(2,mw)=lambda(mw)
         r(3,mw)=(lambda(mw))**2
      enddo
      if (.not.rarecorrector) then
         lambda(2) = 0.5d0*(lambda(1)+lambda(3))
         !lambda(2) = max(min(0.5d0*(s1m+s2m),sE2),sE1)
         r(1,2)=0.d0
         r(2,2)=0.d0
         r(3,2)=1.d0
      endif
      !---------------------------------------------------


      !determine the steady state wave -------------------
      !criticaltol = 1.d-6
      ! MODIFIED:
      criticaltol = max(drytol*g, 1d-6)
      criticaltol_2 = sqrt(criticaltol)
      deldelh = -delb
      deldelphi = -0.5d0 * (hR + hL) * (g * delb + delp / rho)

      !determine a few quanitites needed for steady state wave if iterated
      hLstar=hL
      hRstar=hR
      uLstar=uL
      uRstar=uR
      huLstar=uLstar*hLstar
      huRstar=uRstar*hRstar

      !iterate to better determine the steady state wave
      convergencetol=1.d-6
!      do iter=1,maxiter  ! since maxiter=1, no loop needed!
         !determine steady state wave (this will be subtracted from the delta vectors)
         if (min(hLstar,hRstar).lt.drytol.and.rarecorrector) then
            rarecorrector=.false.
            hLstar=hL
            hRstar=hR
            uLstar=uL
            uRstar=uR
            huLstar=uLstar*hLstar
            huRstar=uRstar*hRstar
            lambda(2) = 0.5d0*(lambda(1)+lambda(3))
            !lambda(2) = max(min(0.5d0*(s1m+s2m),sE2),sE1)
            r(1,2)=0.d0
            r(2,2)=0.d0
            r(3,2)=1.d0
         endif

         hbar =  max(0.5d0*(hLstar+hRstar),0.d0)
         s1s2bar = 0.25d0*(uLstar+uRstar)**2 - g*hbar
         s1s2tilde= max(0.d0,uLstar*uRstar) - g*hbar

         !find if sonic problem
         ! MODIFIED from 5.3.1 version
         sonic = .false.
         if (abs(s1s2bar) <= criticaltol) then
            sonic = .true.
         else if (s1s2bar*s1s2tilde <= criticaltol**2) then
            sonic = .true.
         else if (s1s2bar*sE1*sE2 <= criticaltol**2) then
            sonic = .true.
         else if (min(abs(sE1),abs(sE2)) < criticaltol_2) then
            sonic = .true.
         else if (sE1 <  criticaltol_2 .and. s1m > -criticaltol_2) then
            sonic = .true.
         else if (sE2 > -criticaltol_2 .and. s2m <  criticaltol_2) then
            sonic = .true.
         else if ((uL+sqrt_ghL)*(uR+sqrt_ghR) < 0.d0) then
            sonic = .true.
         else if ((uL- sqrt_ghL)*(uR- sqrt_ghR) < 0.d0) then
            sonic = .true.
         end if

         !find jump in h, deldelh
         if (sonic) then
            deldelh =  -delb
         else
            deldelh = delb*g*hbar/s1s2bar
         endif
         !find bounds in case of critical state resonance, or negative states
         if (sE1.lt.-criticaltol.and.sE2.gt.criticaltol) then
            deldelh = min(deldelh,hstarHLL*(sE2-sE1)/sE2)
            deldelh = max(deldelh,hstarHLL*(sE2-sE1)/sE1)
         elseif (sE1.ge.criticaltol) then
            deldelh = min(deldelh,hstarHLL*(sE2-sE1)/sE1)
            deldelh = max(deldelh,-hL)
         elseif (sE2.le.-criticaltol) then
            deldelh = min(deldelh,hR)
            deldelh = max(deldelh,hstarHLL*(sE2-sE1)/sE2)
         endif

         !find jump in phi, deldelphi
         if (sonic) then
            deldelphi = -g*hbar*delb
         else
            deldelphi = -delb*g*hbar*s1s2tilde/s1s2bar
         endif
         !find bounds in case of critical state resonance, or negative states
         deldelphi=min(deldelphi,g*max(-hLstar*delb,-hRstar*delb))
         deldelphi=max(deldelphi,g*min(-hLstar*delb,-hRstar*delb))
         deldelphi = deldelphi - hbar * delp / rho

         del(1)=delh-deldelh
         del(2)=delhu
         del(3)=delphi-deldelphi

         !Determine determinant of eigenvector matrix========
         det1=r(1,1)*(r(2,2)*r(3,3)-r(2,3)*r(3,2))
         det2=r(1,2)*(r(2,1)*r(3,3)-r(2,3)*r(3,1))
         det3=r(1,3)*(r(2,1)*r(3,2)-r(2,2)*r(3,1))
         determinant=det1-det2+det3

         !solve for beta(k) using Cramers Rule=================
         do k=1,3
            do mw=1,3
                  A(1,mw)=r(1,mw)
                  A(2,mw)=r(2,mw)
                  A(3,mw)=r(3,mw)
            enddo
            A(1,k)=del(1)
            A(2,k)=del(2)
            A(3,k)=del(3)
            det1=A(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2))
            det2=A(1,2)*(A(2,1)*A(3,3)-A(2,3)*A(3,1))
            det3=A(1,3)*(A(2,1)*A(3,2)-A(2,2)*A(3,1))
            beta(k)=(det1-det2+det3)/determinant
         enddo

         !exit if things aren't changing -> removed to allow vectorization!
         !if (abs(del(1)**2+del(3)**2-delnorm).lt.convergencetol) exit
         delnorm = del(1)**2+del(3)**2
         !find new states qLstar and qRstar on either side of interface
         hLstar=hL
         hRstar=hR
         uLstar=uL
         uRstar=uR
         huLstar=uLstar*hLstar
         huRstar=uRstar*hRstar
         do mw=1,mwaves
            if (lambda(mw).lt.0.d0) then
               hLstar= hLstar + beta(mw)*r(1,mw)
               huLstar= huLstar + beta(mw)*r(2,mw)
            endif
         enddo
         do mw=mwaves,1,-1
            if (lambda(mw).gt.0.d0) then
               hRstar= hRstar - beta(mw)*r(1,mw)
               huRstar= huRstar - beta(mw)*r(2,mw)
            endif
         enddo

         if (hLstar.gt.drytol) then
            uLstar=huLstar/hLstar
         else
            hLstar=max(hLstar,0.d0)
            uLstar=0.d0
         endif
         if (hRstar.gt.drytol) then
            uRstar=huRstar/hRstar
         else
            hRstar=max(hRstar,0.d0)
            uRstar=0.d0
         endif

  !    enddo ! end iteration on Riemann problem  ! since maxiter=1, no loop needed!

!       do mw=1,mwaves
!          sw(mw)=lambda(mw)
!          fw(1,mw)=beta(mw)*r(2,mw)
!          fw(2,mw)=beta(mw)*r(3,mw)
!          fw(3,mw)=beta(mw)*r(2,mw)
!       enddo
      sw1 = lambda(1)
      sw2 = lambda(2)
      sw3 = lambda(3)
      !mw=1
      fw11 = beta(1)*r(2,1)
      fw21 = beta(1)*r(3,1)
      fw31 = beta(1)*r(2,1)
      !mw=2
      fw12 = beta(2)*r(2,2)
      fw22 = beta(2)*r(3,2)
      fw32 = beta(2)*r(2,2)
      !mw=3
      fw13 = beta(3)*r(2,3)
      fw23 = beta(3)*r(3,3)
      fw33 = beta(3)*r(2,3)

      !find transverse components (ie huv jumps).
      ! MODIFIED from 5.3.1 version
!       fw(3,1)=fw(3,1)*vL
!       fw(3,3)=fw(3,3)*vR
!       fw(3,2)= 0.d0
      fw31=fw31*vL
      fw33=fw33*vR
      fw32=0.d0
 
!       hustar_interface = huL + fw(1,1)   ! = huR - fw(1,3)
!       if (hustar_interface <= 0.0d0) then
!           fw(3,1) = fw(3,1) + (hR*uR*vR - hL*uL*vL - fw(3,1)- fw(3,3))
!         else
!           fw(3,3) = fw(3,3) + (hR*uR*vR - hL*uL*vL - fw(3,1)- fw(3,3))
!       end if
      hustar_interface = huL + fw11   ! = huR - fw(1,3)
      if (hustar_interface <= 0.0d0) then
          fw31 = fw31 + (hR*uR*vR - hL*uL*vL - fw31- fw33)
        else
          fw33 = fw33 + (hR*uR*vR - hL*uL*vL - fw31- fw33)
      end if


      return

      end !subroutine riemann_aug_JCP-------------------------------------------------


!=============================================================================
      subroutine riemanntype(hL,hR,uL,uR,hm,s1m,s2m,rare1,rare2, &
                  maxiter,drytol,g)

      !determine the Riemann structure (wave-type in each family)


      implicit none

      !input
      real(kind=8) hL,hR,uL,uR,drytol,g
      integer maxiter

      !output
      real(kind=8) s1m,s2m
      logical rare1,rare2

      !local
      real(kind=8) hm,u1m,u2m,um,delu
      real(kind=8) h_max,h_min,h0,F_max,F_min,dfdh,F0,slope,gL,gR
      integer iter
      real(kind=8) sqrt_ghL, sqrt_ghR



      !Test for Riemann structure

      h_min=min(hR,hL)
      h_max=max(hR,hL)
      delu=uR-uL
      
      ! pre-compute square roots:
      sqrt_ghL = sqrt(g*hL)
      sqrt_ghR = sqrt(g*hR)

      if (h_min.le.drytol) then
         hm=0.d0
         um=0.d0
         s1m=uR+uL-2.d0*sqrt_ghR+2.d0*sqrt_ghL
         s2m=uR+uL-2.d0*sqrt_ghR+2.d0*sqrt_ghL
         if (hL.le.0.d0) then
            rare2=.true.
            rare1=.false.
         else
            rare1=.true.
            rare2=.false.
         endif

      else
         F_min= delu+2.d0*(sqrt(g*h_min)-sqrt(g*h_max))
         F_max= delu + &
              (h_max-h_min)*(sqrt(.5d0*g*(h_max+h_min)/(h_max*h_min))) 

         if (F_min.gt.0.d0) then !2-rarefactions

            hm=(1.d0/(16.d0*g))* &
                    max(0.d0,-delu+2.d0*(sqrt_ghL+sqrt_ghR))**2
            um=sign(1.d0,hm)*(uL+2.d0*(sqrt_ghL-sqrt(g*hm)))

            s1m=uL+2.d0*sqrt_ghL-3.d0*sqrt(g*hm)
            s2m=uR-2.d0*sqrt_ghR+3.d0*sqrt(g*hm)

            rare1=.true.
            rare2=.true.

         elseif (F_max.le.0.d0) then !2 shocks

            !root finding using a Newton iteration on sqrt(h)===
            h0=h_max
            do iter=1,maxiter
               gL=sqrt(.5d0*g*(1/h0 + 1/hL))
               gR=sqrt(.5d0*g*(1/h0 + 1/hR))
               F0=delu+(h0-hL)*gL + (h0-hR)*gR
               dfdh=gL-g*(h0-hL)/(4.d0*(h0**2)*gL)+ &
                        gR-g*(h0-hR)/(4.d0*(h0**2)*gR)
               slope=2.d0*sqrt(h0)*dfdh
               h0=(sqrt(h0)-F0/slope)**2
            enddo
               hm=h0
               u1m=uL-(hm-hL)*sqrt((.5d0*g)*(1/hm + 1/hL))
               u2m=uR+(hm-hR)*sqrt((.5d0*g)*(1/hm + 1/hR))
               um=.5d0*(u1m+u2m)

               s1m=u1m-sqrt(g*hm)
               s2m=u2m+sqrt(g*hm)
               rare1=.false.
               rare2=.false.

         else !one shock one rarefaction
            h0=h_min

            do iter=1,maxiter
               F0=delu + 2.d0*(sqrt(g*h0)-sqrt(g*h_max)) &
                       + (h0-h_min)*sqrt(.5d0*g*(1/h0+1/h_min))
               slope=(F_max-F0)/(h_max-h_min)
               h0=h0-F0/slope
            enddo

            hm=h0
            if (hL.gt.hR) then
               um=uL+2.d0*sqrt_ghL-2.d0*sqrt(g*hm)
               s1m=uL+2.d0*sqrt_ghL-3.d0*sqrt(g*hm)
               s2m=uL+2.d0*sqrt_ghL-sqrt(g*hm)
               rare1=.true.
               rare2=.false.
            else
               s2m=uR-2.d0*sqrt_ghR+3.d0*sqrt(g*hm)
               s1m=uR-2.d0*sqrt_ghR+sqrt(g*hm)
               um=uR-2.d0*sqrt_ghR+2.d0*sqrt(g*hm)
               rare2=.true.
               rare1=.false.
            endif
         endif
      endif

      return

      end ! subroutine riemanntype----------------------------------------------------------------