!========================================================================
! BEARCLAW Boundary Embedded Adaptive Refinement Conservation LAW package
!========================================================================
! (c) Copyright Sorin Mitran, 2000
! Department of Applied Mathematics
! University of Washington
! mitran@amath.washington.edu
! -----------------------------------------------------------------------
! This code may be freely used for educational and research purposes.
! For any other use please contact the author.
! -----------------------------------------------------------------------
! File:             problem.f90
! Purpose:          Definition of problem specific data and procedures
!                   Contains:
!                      setprob  - Initialize problem parameters
!                      qinit    - Initialize field variables
!                      b4step   - Called prior to each time step
!                      setaux   - Initialize auxilliary variable fields
!                      src      - Compute source terms
!                      physflux - Compute the physical flux, related quantities
!                                 Compute wave decomposition
!
! Comments:         This module is problem specific and must be written by the
!                   user. For exampels see the BEARCLAW documentation and also
!                   the examples and applications directories
!                 
!                   The encapsulation of all problem specific data and routines
!                   and data into a module is useful in exposing data that needs
!                   to be shared.
! Application:      Euler equations, shock tube problem
! Contains:
! Revision History: Ver. 1.0 Oct. 2001 Sorin Mitran
! -----------------------------------------------------------------------

!*************
MODULE problem
!*************
  ! User definitions of data structures associated with each node
  USE NodeInfoDef
  IMPLICIT NONE    ! It's safer to require explicit declarations
  SAVE             ! Save module information
  PRIVATE          ! Everything is implicitly private, i.e. accessible only
                   ! to this module. 
  ! These are the public entry points (acessible to other modules)                   
  PUBLIC setprob,qinit,b4step,setaux,src,physflux,problemBC,afterstep,problemIO
  
  ! Problem specific global variables:  
  !  - Problem parameters
  REAL (KIND=qPrec) :: gamma,gamma1
  REAL (KIND=qPrec) :: rhol,rhoul,el,rhor,rhour,er
  
  !
  INTEGER, PARAMETER :: OK=0
  REAL (KIND=qPrec), PARAMETER :: zero=0., half=0.5, one=1., two=2.

  !=======
  CONTAINS
  !=======
  
    SUBROUTINE setprob  
      ! User definitions of data structures associated with each node
      USE NodeInfoDef      
      REAL (KIND=qPrec) :: ul,pl,ur,pr
      ! Set the initial shock tube states
      OPEN(UNIT=7,FILE='setprob.data',STATUS='old',FORM='formatted')  
      READ(7,*) gamma
      gamma1 = gamma - 1.d0
      READ(7,*) rhol,ul,pl
      READ(7,*) rhor,ur,pr
      ! data in left state
      rhoul = rhol*ul
      el = pl/gamma1 + 0.5d0*rhol*ul**2
      ! data in right state
      rhour = rhor*ur
      er = pr/gamma1 + 0.5d0*rhor*ur**2
      CLOSE(7)
    END SUBROUTINE setprob

    SUBROUTINE problemBC(Info)
      ! User definitions of data structures associated with each node
      USE NodeInfoDef  
      ! Interface declarations      
      TYPE (NodeInfo) :: Info  ! Data associated with this grid      
      ! Internal declarations
    END SUBROUTINE problemBC

    SUBROUTINE qinit(Info)
      ! User definitions of data structures associated with each node
      USE NodeInfoDef  
      ! Interface declarations      
      TYPE (NodeInfo) :: Info  ! Data associated with this grid  
      ! Internal declarations
      INTEGER i
      REAL (KIND=xPrec), PARAMETER :: xDiscontinuity = 0.5
      REAL (KIND=xPrec) :: x      
      !
      x=Info%Xlower(1)+Info%dX(1)/2.0
      DO i=1,Info%mX(1)
        IF (x<xDiscontinuity) THEN
          !      x y z t m
          Info%q(i,1,1,1,1) = rhol
          Info%q(i,1,1,1,2) = rhoul
          Info%q(i,1,1,1,3) = el
        ELSE
          !      x y z t m
          Info%q(i,1,1,1,1) = rhor
          Info%q(i,1,1,1,2) = rhour
          Info%q(i,1,1,1,3) = er
        END IF
        x=x+Info%dX(1)
      END DO    
    END SUBROUTINE qinit
    
    SUBROUTINE b4step(Info)  
      ! Interface declarations    
      TYPE (NodeInfo) :: Info        
    END SUBROUTINE b4step
    
    SUBROUTINE afterstep(Info)
      ! Interface declarations    
      TYPE (NodeInfo) :: Info
      ! Include operations that need to be carried out after each time step here:
    END SUBROUTINE afterstep
    
    SUBROUTINE problemIO(nframe,tnow,IOrequest,Info,qmax,qmin)
      ! Interface declarations
      INTEGER :: nframe,IOrequest; REAL (KIND=qPrec) :: tnow
      TYPE (NodeInfo), OPTIONAL :: Info      
      REAL (KIND=qPrec), OPTIONAL, DIMENSION(:) :: qmax,qmin
      INTEGER, PARAMETER :: UserBeforeGridIO=-1, UserAfterGridIO=-2
      INTEGER, PARAMETER :: UserIOMinMax=0
      ! This routine provides for output of computation snapshots in user-defined
      ! formats.
      ! IOrequest is a function parameter with values:
      ! -1 UserBeforeGridIO - indicates routine is being called prior to parsing
      !                       of tree of grids. This is where output files
      !                       should be opened
      ! -2 UserAfterGridIO  - indicates routine is being called after parsing
      !                       of tree of grids. This is where output files
      !                       should be closed
      ! -3 UserOutputq      - indicates routine should output the field
      !                       variables in user specified format
      ! 0  UserIOMinMax     - routine is being called after qmin and qmax have
      !                       been determined. These are the min/max field
      !                       variable values. If output of quantities different
      !                       from the field variables is desired, qmin &
      !                       qmax should be defined with respect to those quantities
      ! 1,2 ...             - routine is being called to output field variable
      !                       1,2,... on the current grid in conjunction with
      !                       HDF support routines in BEARIO. Only valid for
      !                       UserHDF output style      
    END SUBROUTINE problemIO
    
    SUBROUTINE setaux(Info)  
      ! Interface declarations    
      TYPE (NodeInfo) :: Info        
    END SUBROUTINE setaux
    
    SUBROUTINE src(Info)
      ! Interface declarations    
      TYPE (NodeInfo) :: Info        
    END SUBROUTINE src
    
! -----------------------------------------------------------------------
!
! Compute the Riemann problem solution for the 1D Euler equations
!
!    q_t + f(q)_x = 0
!
! Routine should return wave decomposition resulting from solving Riemann
! problem at interfaces j=2-mbc:mxnow+mbc, between left states in
! q1D(1-mbc:mxnow+mbc-1) and right states in q1D(2-mbc:mxnow+mbc)
! -----------------------------------------------------------------------
    SUBROUTINE physflux(ixy,indx,irequest,Info,q,f,s,A)
      ! User definitions of data structures associated with each node
      USE NodeInfoDef  
      ! Interface declarations    
      INTEGER ixy           ! The dimension along which to solve the normal Riemann problem    
      INTEGER indx(MaxDims) ! Indices of this slice, indx(ixy) has no significance, other
                            ! indices give position of this 1D slice in index space
      INTEGER irequest      ! Request code, see NodeInfoGlobal.f90 for definitions
                            ! Additional request codes may be defined in solver modules
      TYPE (NodeInfo) :: Info  ! Data associated with this grid
      REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: q,f,s
      REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: A
      ! Internal declarations      
      LOGICAL, PARAMETER :: efix = .TRUE.
      INTEGER mbc,mx,nq,n,mw,NrVars,mwaves,iError,j
      ! Shorter local names for the wave propagation quantities
      REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: Apdq,Amdq,speed,q1D
      REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: wave
      REAL (KIND=qPrec) :: rho1,rhou1,en1,p1,c1,rho2,rhou2,en2,p2,c2
      REAL (KIND=qPrec) :: s0,s1,s2,s3,sfract,df(3)
      ! Some additional local fields
      ! Cell center quantities
      REAL (KIND=qPrec), ALLOCATABLE, SAVE, DIMENSION (:) :: sqrtrho, &
           pres,rhoRsq2,delta,coef,rho,u,c,h
      ! Roe averaged quantities
      REAL (KIND=qPrec), ALLOCATABLE, SAVE, DIMENSION (:) :: uR,hR,cR
      !
      ! Current number of interior cells, boundary cells, field variables, waves
      mx=Info%mXnow; mbc=Info%mbc; NrVars=Info%NrVars; mwaves=Info%mwaves            
      ! Associate local names with Info components
      q1D=>Info%q1D
      Apdq=>Info%Apdq; Amdq=>Info%Amdq; speed=>Info%speed; wave=>Info%wave
      Apdq=0.; Amdq=0.; speed=0.; wave=0.      
      SELECT CASE (irequest)
        CASE (Initialize)
          ! Allocate local work space
          ALLOCATE(delta(NrVars),coef(mwaves),                                   &
                   sqrtrho(1-mbc:mx+mbc),pres(1-mbc:mx+mbc),uR(1-mbc:mx+mbc), &
                   hR(1-mbc:mx+mbc),cR(1-mbc:mx+mbc),rhoRsq2(1-mbc:mx+mbc),     &
                   h(1-mbc:mx+mbc),c(1-mbc:mx+mbc),rho(1-mbc:mx+mbc),     &
                   u(1-mbc:mx+mbc),     &
                   STAT=iError)
          IF (iError/=OK) THEN
            PRINT *,'Cannot allocate work space in problem module, physflux'
            STOP
          END IF
        CASE (Finalize)
          ! Release local work space
          DEALLOCATE(delta,coef,sqrtrho,pres,uR,hR,cR,rhoRsq2,h,c,rho,u,STAT=iError)
          IF (iError/=OK) THEN
            PRINT *,'Cannot allocate work space in problem module, physflux'
            STOP
          END IF
        CASE (RequestFluxes)
          IF (MINVAL(q1D(1-mbc:mx+mbc,1))<=zero) THEN
            PRINT *,'Error: Negative or zero density in physflux.'
            PRINT '(1x,a6,1x,32(f7.4,1x))','rho=  ',q1D(1-mbc:mx+mbc,1)
            STOP
          END IF
          ! Compute cell centered quantities required in Riemann solve & entropy fix
          rho(1-mbc:mx+mbc) = q1D(1-mbc:mx+mbc,1)
          u(1-mbc:mx+mbc) = q1D(1-mbc:mx+mbc,2)/rho(1-mbc:mx+mbc)
          sqrtrho(1-mbc:mx+mbc) = SQRT(rho(1-mbc:mx+mbc))
          pres(1-mbc:mx+mbc) = gamma1 *                    &
            ( q1D(1-mbc:mx+mbc,3) - half*q1D(1-mbc:mx+mbc,2)*u(1-mbc:mx+mbc) )
          c(1-mbc:mx+mbc) = gamma*pres(1-mbc:mx+mbc)/rho(1-mbc:mx+mbc)
          IF (MINVAL(c(1-mbc:mx+mbc))<=zero) THEN
            PRINT *,'Error: Non-physical centered sound velocity in physflux.'
            PRINT '(1x,a6,1x,32(f7.4,1x))','c**2 =  ',c(1-mbc:mx+mbc)
            STOP
          END IF
          c(1-mbc:mx+mbc) = SQRT(c(1-mbc:mx+mbc))
          h(1-mbc:mx+mbc) = (q1D(1-mbc:mx+mbc,3) + pres(1-mbc:mx+mbc)) / &
                            rho(1-mbc:mx+mbc)
          ! Compute the physical fluxes at cell centers
          f(1-mbc:mx+mbc,1)=q1D(1-mbc:mx+mbc,2)
          f(1-mbc:mx+mbc,2)=f(1-mbc:mx+mbc,1)*u(1-mbc:mx+mbc)+pres(1-mbc:mx+mbc)
          f(1-mbc:mx+mbc,3)=f(1-mbc:mx+mbc,1)*h(1-mbc:mx+mbc)          
        CASE (RequestNormalWaves)          
          ! Compute Roe averages. Cell centered quantities are available from 
          ! previous RequestFluxes code, guaranteed to be executed before
          ! RequestNormalWaves (notation: uR=Roe average, u=cell value, etc.)               
          rhoRsq2(2-mbc:mx+mbc) = one / (sqrtrho(1-mbc:mx+mbc-1)+sqrtrho(2-mbc:mx+mbc))
          uR(2-mbc:mx+mbc) = ( q1D(1-mbc:mx+mbc-1,2)/sqrtrho(1-mbc:mx+mbc-1) + &
                              q1D(2-mbc:mx+mbc,2)/sqrtrho(2-mbc:mx+mbc) )  *  &
                              rhoRsq2(2-mbc:mx+mbc)
          hR(2-mbc:mx+mbc) =                                                          &
            ( (q1D(1-mbc:mx+mbc-1,3)+pres(1-mbc:mx+mbc-1))/sqrtrho(1-mbc:mx+mbc-1) +   &
              (q1D(2-mbc:mx+mbc,3)+pres(2-mbc:mx+mbc))/sqrtrho(2-mbc:mx+mbc)       ) * &
            rhoRsq2(2-mbc:mx+mbc)          
          cR = gamma1*(hR-half*uR**2)          
          IF (MINVAL(cR(2-mbc:mx+mbc))<=zero) THEN
            PRINT *,'Error: Negative or zero Roe average sound velocity in physflux.'            
            STOP
          END IF          
          cR(2-mbc:mx+mbc) = SQRT(cR(2-mbc:mx+mbc))
          ! Construct the eigenbasis 
          wave(2-mbc:mx+mbc,1:NrVars,1:mwaves)=one
          
          speed(2-mbc:mx+mbc,1) = uR(2-mbc:mx+mbc) - cR(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,2,1) = speed(2-mbc:mx+mbc,1)
          wave(2-mbc:mx+mbc,3,1) = hR(2-mbc:mx+mbc) - uR(2-mbc:mx+mbc)*cR(2-mbc:mx+mbc)
          
          speed(2-mbc:mx+mbc,2) = uR(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,2,2) = uR(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,3,2) = half*uR(2-mbc:mx+mbc)**2
          
          speed(2-mbc:mx+mbc,3) = uR(2-mbc:mx+mbc) + cR(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,2,3) = speed(2-mbc:mx+mbc,3)
          wave(2-mbc:mx+mbc,3,3) = hR(2-mbc:mx+mbc) + uR(2-mbc:mx+mbc)*cR(2-mbc:mx+mbc)
          
          ! Decompose jump in q onto eigenbasis
          
          DO j=2-mbc,mx+mbc
            ! Find coef(1) thru coef(3), the coefficients of the 3 eigenvectors
            delta(:) = q1D(j,:) - q1D(j-1,:)
            coef(2) = gamma1/cR(j)**2 * ((hR(j)-uR(j)**2)*delta(1) + uR(j)*delta(2) - delta(3))
            coef(3) = (delta(2) + (cR(j)-uR(j))*delta(1) - cR(j)*coef(2)) / (two*cR(j))
            coef(1) = delta(1) - coef(2) - coef(3)            
            DO mw=1,mwaves
              wave(j,1:NrVars,mw) = coef(mw)*wave(j,1:NrVars,mw)
            END DO
          END DO          
          
          IF (efix) THEN
            ! With entropy fix
            ! Compute flux differences amdq and apdq.
            ! First compute amdq as sum of s*wave for left going waves.
            ! Incorporate entropy fix by adding a modified fraction of wave
            ! if s should change sign.
            DO j=2-mbc,mx+mbc
              ! 1-wave. Check for fully supersonic case
              s0=u(j-1)-c(j-1)
              IF ( s0>=zero .AND. speed(j,1) > zero ) THEN
                ! Everything is right-going
                Amdq(j,1:NrVars) = zero  
                CYCLE
              END IF
              ! u-c to right of 1-wave
              rho1 = q1D(j-1,1) + wave(j,1,1)
              rhou1 = q1D(j-1,2) + wave(j,2,1)
              en1 = q1D(j-1,3) + wave(j,3,1)
              p1 = gamma1*(en1  - 0.5*rhou1**2/rho1)
              c1 = SQRT(gamma*p1/rho1)
              s1 = rhou1/rho1 - c1
              IF ( s0<zero .AND. s1>0 ) THEN
                ! Transonic rarefaction in the 1-wave                
                sfract = s0 * (s1-speed(j,1)) / (s1-s0)
              ELSE IF (speed(j,1) < zero) THEN
                ! Left-going 1-wave
                sfract=speed(j,1)
              ELSE
                ! Right-going 1-wave
                sfract=zero  ! Never should reach this instruction
              END IF
              Amdq(j,1:NrVars) = sfract*wave(j,1:NrVars,1)
              
              ! Check 2-wave
              IF (speed(j,2) >= zero) CYCLE ! 2-wave is right-going
              Amdq(j,1:NrVars) = Amdq(j,1:NrVars) + speed(j,2)*wave(j,1:NrVars,2)
              
              ! Check 3-wave. Compute u+c to left of 3-wave
              rho2 = q1D(j,1) - wave(j,1,3)
              rhou2 = q1D(j,2) - wave(j,2,3)
              en2 = q1D(j,3) - wave(j,3,3)
              p2 = gamma1*(en2 - 0.5d0*rhou2**2/rho2)
              c2 = dsqrt(gamma*p2/rho2)
              s2 = rhou2/rho2 + c2; s3=u(j)+c(j)
              IF ( s2 < zero .AND. s3 > zero ) THEN
                ! Transonic rarefaction in the 3-wave
                sfract = s2 * (s3-speed(j,3)) / (s3-s2)
              ELSE IF (speed(j,3) < zero) THEN
                ! Left-going 3-wave
                sfract = speed(j,3)
              ELSE
                ! Right-going 3-wave
                CYCLE
              END IF
              Amdq(j,1:NrVars) = Amdq(j,1:NrVars) + sfract*wave(j,1:NrVars,3)
              
            END DO

            ! Compute the right-going flux-differences
            DO j=2-mbc,mx+mbc
              df=zero
              DO mw=1,mwaves
                df(1:NrVars) = df(1:NrVars) + speed(j,mw) * wave(j,1:NrVars,mw)
              END DO
              Apdq(j,1:NrVars) = df(1:NrVars) - Amdq(j,1:NrVars)
            END DO            
          !===
          ELSE
          !===
            ! No entropy fix
            ! amdq = SUM s*wave   over left-going waves
            ! apdq = SUM s*wave   over right-going waves            
            DO mw=1,mwaves              
              DO j=2-mbc,mx+mbc
                IF (speed(j,mw) < zero) THEN
                  Amdq(j,1:NrVars) = Amdq(j,1:NrVars) +  speed(j,mw) * wave(j,1:NrVars,mw)
                ELSE
                  Apdq(j,1:NrVars) = Apdq(j,1:NrVars) +  speed(j,mw) * wave(j,1:NrVars,mw)
                END IF              
              END DO
            END DO            
          END IF
        CASE DEFAULT          
          PRINT *,'Error: the solver module Scheme.f90 emitted a request that'
          PRINT *,'is not implemented in the problem module problem.f90'
          STOP
      END SELECT      
    END SUBROUTINE physflux
              
!*****************
END MODULE problem
!*****************