Main Page
src/generic/Generic_RKstep.f90

SUBROUTINE OneDRKstep(rk,ux,uxold,dt,cfl,fxi,dcflag,ifilter)

  ! ----  Shared Variables
  USE mesh
  USE array_bounds
  USE method_parms
  USE Interp_coeffs

  ! ----  Shared Procedures
  USE Generic_GetRHS
  USE Generic_CFL
  USE Generic_Transport
  USE Generic_ModifyDCflag

  IMPLICIT NONE

  include 'cles.i'

  DOUBLE PRECISION, INTENT(INOUT) :: uxold(ncomps,ixlo:ixhi)
  DOUBLE PRECISION, INTENT(IN) :: dt
  INTEGER, INTENT(IN) :: rk, ifilter

  DOUBLE PRECISION, INTENT(INOUT) :: ux(ncomps,ixlo:ixhi)
  DOUBLE PRECISION, INTENT(OUT) :: fxi(ncomps,ixlo:ixhi,1)
  INTEGER, INTENT(INOUT) :: dcflag(1:nx+1,1)
  DOUBLE PRECISION, INTENT(OUT) :: cfl

  DOUBLE PRECISION, ALLOCATABLE :: vx(:,:)
  DOUBLE PRECISION, ALLOCATABLE :: rhs(:,:)

  DOUBLE PRECISION :: c_rk1, c_rk2, c_rk3, f_rk

  INTEGER :: iret

  ALLOCATE ( vx(nvars,ixlo:ixhi) , stat = iret )
  if ( iret .ne. 0 ) then
     call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
          'OneDRKstep/vx', iret)
  endif
  ALLOCATE ( rhs(nvars,1:nx), stat=iret )
  if ( iret .ne. 0 ) then
     call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
          'OneDRKstep/rhs', iret)
  endif

  if ( ifilter .eq. CLES_TRUE ) use_dcflag = CLES_FALSE

  IF (useViscous.eq.CLES_TRUE) Call AllocateTransport()

  ! ---- Set the constants for this (number rk) step
  CALL SetRKConstants(rk, c_rk1, c_rk2, c_rk3, f_rk)

  ! ---- course fine boundary flags from AMROC are 
  ! ---- re-interperated in the standard DCflag paradigm
  if ( rk .eq. 1 ) then
     CALL ModifyDCflagFromAMR(dcflag)
  else
     dcflag(:,1) = uxold(nvars+2,1:nx+1)
  endif

  ! ---- Get the right hand side (RHS) 
  CALL GetRHS(rk, ux,vx,rhs,fxi,dcflag,ifilter)

  ! ---- note, this only updates the physical (interior) points
  ! ---- and leaves the ghostcells alone
  
  ux(1:nvars,1:nx) = c_rk1*dt*rhs(1:nvars,:) &
       &           + c_rk2*ux(1:nvars,1:nx)    &
       &           + c_rk3*uxold(1:nvars,1:nx) 

  ! --- weight the flux so that the 'fix-up' in amroc will be 
  ! --- correct.  This flux is exposed to AMROC.
  if ( noTimeRefine .eq. 1 ) then
     fxi = fxi * c_rk1
  else
     fxi = fxi * f_rk
  endif

  IF (rk.EQ.3) THEN
     Call Eval_CFL(dt,cfl,ux,vx)
     
     ux(nvars+2,1:nx) = 0.5D0*(dcflag(1:nx,1) + dcflag(2:nx+1,1))
  else if ( rk .eq. 1 ) then
     uxold(nvars+2,1:nx+1) = dcflag(:,1)+0.5d0
  END IF

  DEALLOCATE ( rhs, stat=iret )
  if ( iret .ne. 0 ) then
     call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
          'OneDRKstep/rhs', iret)
  endif
  DEALLOCATE(vx, stat=iret )
  if ( iret .ne. 0 ) then
     call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
          'OneDRKstep/vx', iret)
  endif
  
  IF (useViscous.eq.CLES_TRUE) Call DeallocateTransport()

  ! to be safe, we deactivate extended output to avoid a residual
  ! effect of calling this subroutine and then the output one but
  ! on a different patch.
  useExOutput = CLES_FALSE

END SUBROUTINE OneDRKstep


SUBROUTINE TwoDRKstep(rk,ux,uxold,dt,cfl,fxi,dcflag,ifilter)

  ! ----  Shared Variables
  USE mesh
  USE array_bounds
  USE method_parms
  USE Interp_coeffs

  ! ----  Shared Procedures
  USE Generic_GetRHS
  USE Generic_CFL
  USE Generic_Transport
  USE Generic_ModifyDCflag

  IMPLICIT NONE

  include 'cles.i'

  DOUBLE PRECISION, INTENT(INOUT) :: uxold(ncomps,ixlo:ixhi,iylo:iyhi)
  DOUBLE PRECISION, INTENT(IN) :: dt
  INTEGER, INTENT(IN) :: rk, ifilter

  DOUBLE PRECISION, INTENT(OUT) :: ux(ncomps,ixlo:ixhi,iylo:iyhi)
  DOUBLE PRECISION, INTENT(OUT) :: fxi(ncomps,ixlo:ixhi,iylo:iyhi,2)
  INTEGER, INTENT(INOUT) :: dcflag(1:nx+1,1:ny+1,2)
  DOUBLE PRECISION, INTENT(OUT) :: cfl

  DOUBLE PRECISION, ALLOCATABLE :: vx(:,:,:)
  DOUBLE PRECISION, ALLOCATABLE :: rhs(:,:,:)

  DOUBLE PRECISION :: c_rk1, c_rk2, c_rk3, f_rk
  
  INTEGER :: iret

  ALLOCATE ( vx(nvars,ixlo:ixhi,iylo:iyhi), stat=iret ) 
  if ( iret .ne. 0 ) then
     call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
          'TwoDRKstep/vx', iret)
  endif
  ALLOCATE ( rhs(nvars,1:nx,1:ny), stat=iret )
  if ( iret .ne. 0 ) then
     call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
          'TwoDRKstep/rhs', iret)
  endif

  if ( ifilter .eq. CLES_TRUE ) use_dcflag = CLES_FALSE

  IF (useViscous.eq.CLES_TRUE) Call AllocateTransport()

  CALL SetRKConstants(rk, c_rk1, c_rk2, c_rk3, f_rk)

  ! ---- course fine boundary flags from AMROC are 
  ! ---- re-interperated in the standard DCflag paradigm
  if ( rk .eq. 1 ) then
     CALL ModifyDCflagFromAMR(dcflag)
  else
     dcflag(:,:,1) = uxold(nvars+2,1:nx+1,1:ny+1)
     dcflag(:,:,2) = MOD(dcflag(:,:,1)/2,2)
     dcflag(:,:,1) = MOD(dcflag(:,:,1),2)
  endif
  
  ! ---- Evaluate the right hand side (RHS)
  CALL GetRHS(rk, ux,vx,rhs,fxi,dcflag,ifilter)

  ! ---- note, this only updates the physical (interior) points
  ! ---- and leaves the ghostcells alone

  ux(1:nvars,1:nx,1:ny) = c_rk1*dt*rhs(1:nvars,:,:)&
       &                 + c_rk2*ux(1:nvars,1:nx,1:ny) & 
       &                 + c_rk3*uxold(1:nvars,1:nx,1:ny) 
  
  ! --- weight the flux so that the 'fix-up' in amroc will be 
  ! --- correct.  This flux is exposed to AMROC.
  if ( noTimeRefine .eq. 1 ) then
     fxi = fxi * c_rk1
  else
     fxi = fxi * f_rk
  endif

  IF (rk.EQ.3) THEN
     Call Eval_CFL(dt,cfl,ux,vx)   

     ux(nvars+2,1:nx,1:ny) = 0.5D0*(dcflag(1:nx,1:ny,1) & 
          + dcflag(2:nx+1,1:ny,1) + dcflag(1:nx,1:ny,2) &
          + dcflag(1:nx,2:ny+1,2) )
  else if ( rk .eq. 1 ) then
     uxold(nvars+2,1:nx+1,1:ny+1) = dcflag(:,:,1)+dcflag(:,:,2)*2+0.5d0
  END IF

  DEALLOCATE ( vx , stat=iret )
  if ( iret .ne. 0 ) then
     call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
          'TwoDRKstep/vx', iret)
  endif
  DEALLOCATE ( rhs, stat=iret )
  if ( iret .ne. 0 ) then
     call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
          'TwoDRKstep/rhs', iret)
  endif

  IF (useViscous.eq.CLES_TRUE) Call DeallocateTransport()

  ! to be safe, we deactivate extended output to avoid a residual
  ! effect of calling this subroutine and then the output one but
  ! on a different patch.
  useExOutput = CLES_FALSE

END SUBROUTINE TwoDRKstep


SUBROUTINE ThreeDRKstep(rk,ux,uxold,dt,cfl,fxi,dcflag, ifilter)

  ! ----  Shared Variables
  USE mesh
  USE array_bounds
  USE method_parms
  USE Interp_coeffs

  ! ----  Shared Procedures
  USE Generic_GetRHS
  USE Generic_CFL
  USE Generic_Transport
  USE Generic_ModifyDCflag

  IMPLICIT NONE

  include 'cles.i'

  DOUBLE PRECISION, INTENT(INOUT) :: uxold(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi)
  DOUBLE PRECISION, INTENT(IN) :: dt
  INTEGER, INTENT(IN) :: rk, ifilter

  DOUBLE PRECISION, INTENT(OUT) :: ux(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi)
  DOUBLE PRECISION, INTENT(OUT) :: fxi(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi,3)
  INTEGER, INTENT(INOUT) :: dcflag(1:nx+1,1:ny+1,1:nz+1,3)
  DOUBLE PRECISION, INTENT(OUT) :: cfl
  
  DOUBLE PRECISION, ALLOCATABLE :: vx(:,:,:,:)
  DOUBLE PRECISION, ALLOCATABLE :: rhs(:,:,:,:)

  DOUBLE PRECISION :: c_rk1, c_rk2, c_rk3,f_rk
  INTEGER :: iret
  
  call cleslog_log_enter('ThreeDRKstep')

  ALLOCATE ( vx(nvars,ixlo:ixhi,iylo:iyhi,izlo:izhi), stat=iret )
  if ( iret .ne. 0 ) then
     call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
          'ThreeDRKstep/vx', iret)
  endif
  ALLOCATE( rhs(nvars,1:nx,1:ny,1:nz), stat=iret )
  if ( iret .ne. 0 ) then
     call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
          'ThreeDRKstep/rhs', iret)
  endif

  if ( ifilter .eq. CLES_TRUE ) use_dcflag = CLES_FALSE
  
  IF (useViscous.eq.CLES_TRUE) Call AllocateTransport()

  ! ---- Create The LES Arrays
  IF (useLES.eq.CLES_TRUE) CALL AllocateLES()

  CALL SetRKConstants(rk, c_rk1, c_rk2, c_rk3, f_rk)

  ! ---- course-fine boundary flags from AMROC are 
  ! ---- re-interperated in the standard DCflag paradigm
  if ( rk .eq. 1 ) then
     CALL ModifyDCflagFromAMR(dcflag)
  else
     dcflag(:,:,:,1) = uxold(nvars+2,1:nx+1,1:ny+1,1:nz+1)
     dcflag(:,:,:,3) = MOD(dcflag(:,:,:,1)/4,2)
     dcflag(:,:,:,2) = MOD(dcflag(:,:,:,1)/2,2)
     dcflag(:,:,:,1) = MOD(dcflag(:,:,:,1),2)
  endif
  
  ! ---- Evaluate the right hand sid (RHS)
  CALL GetRHS(rk, ux,vx,rhs,fxi,dcflag,ifilter)
  
  ! ---- note, this only updates the physical (interior) points
  ! ---- and leaves the ghostcells alone
  
  ux(1:nvars,1:nx,1:ny,1:nz) = c_rk3*uxold(1:nvars,1:nx,1:ny,1:nz) &
       &                     + c_rk2*ux(1:nvars,1:nx,1:ny,1:nz) &
       &                     + c_rk1*dt*rhs(1:nvars,:,:,:)
  
  ! --- weight the flux so that the 'fix-up' in amroc will be 
  ! --- correct.  This flux is exposed to AMROC.
  if ( noTimeRefine .eq. 1 ) then
     fxi = fxi * c_rk1
  else
     fxi = fxi * f_rk
  endif

  if (rk.EQ.3) then
     Call Eval_CFL(dt,cfl,ux,vx)
    
     !---- Extended vector of state:
     ux(nvars+2,1:nx,1:ny,1:nz) = 0.5D0*(dcflag(1:nx,1:ny,1:nz,1) &
          + dcflag(2:nx+1,1:ny,1:nz,1) + dcflag(1:nx,1:ny,1:nz,2) &
          + dcflag(1:nx,2:ny+1,1:nz,2) + dcflag(1:nx,1:ny,1:nz,3) &
          + dcflag(1:nx,1:ny,2:nz+1,3) )
  else if ( rk .eq. 1 ) then
     uxold(nvars+2,1:nx+1,1:ny+1,1:nz+1) = dcflag(:,:,:,1)+dcflag(:,:,:,2)*2 &
          +dcflag(:,:,:,3)*4 + 0.5d0
  endif

  DEALLOCATE ( vx ,stat=iret ) 
  if ( iret .ne. 0 ) then
     call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
          'ThreeDRKstep/vx', iret)
  endif
  DEALLOCATE ( rhs , stat=iret )
  if ( iret .ne. 0 ) then
     call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
          'ThreeDRKstep/rhs', iret)
  endif

  ! --- deallocate the LES arrays
  IF (useLES.eq.CLES_TRUE) CALL DeallocateLES()

  IF (useViscous.eq.CLES_TRUE) Call DeallocateTransport()

  ! to be safe, we deactivate extended output to avoid a residual
  ! effect of calling this subroutine and then the output one but
  ! on a different patch.
  useExOutput = CLES_FALSE

  call cleslog_log_exit('ThreeDRKstep')

END SUBROUTINE ThreeDRKstep