xref: /phasta/phSolver/incompressible/itrPC.f (revision 595995161822a203c8467e0e4a253d7bd7d6df32)

c-----------------------------------------------------------------------
c
c    Initialize the predictor multicorrector (set up parameters)
c
c    Modified by Alberto Figueroa.  Winter 2004
c-----------------------------------------------------------------------
      subroutine itrSetup ( y,  acold )

      include "common.h"

      real*8     y(nshg,ndof),     acold(nshg,ndof)

c
c  Define the Hulbert parameters
c  second order if between (and including) 0 and 1 otherwise backward Euler
c
      if( rhoinf(itseq).lt.0.or.rhoinf(itseq).gt.1) then ! backward Euler
         almi   = one
         alfi   = one
         gami   = one
         ipred  = 1
      else           !second order family
         almi   = (three-rhoinf(itseq))/(one+rhoinf(itseq))/two
         alfi   = one/(one+rhoinf(itseq))
         gami   = pt5+almi-alfi
         if(ideformwall.eq.1) then
            betai=1.0
         else
            betai=0.0
         endif
      endif
c
c.... set the jacobian type
c
      Jactyp=0
      if (impl(itseq) .eq. 3) then
         Jactyp = 1
         impl(itseq) = 2
      endif
c
c.... same_Dy predictor special case
c
      if (ipred.eq.4 .and. itseq .eq. 1 ) then
         y=y-(one-alfi)*Delt(1)*acold
         if ( rhoinf(itseq) .eq. 0.0 ) then
            ipred = 3
         endif
      endif
c
c.... set the global time increment and the CFL data
c
      Dtgl   = one / Delt(itseq)  ! caution: inverse of time step
      CFLfld = CFLfl(itseq)
      CFLsld = CFLsl(itseq)

      return
      end


c-----------------------------------------------------------------------
c
c    Predict solution at time n+1
c
c-----------------------------------------------------------------------
      subroutine itrPredict (yold,  y,  acold,   ac,   uold,   u, iBC)

      include "common.h"

      real*8        y(nshg,ndof),               ac(nshg,ndof),
     &              u(nshg,nsd),                yold(nshg,ndof),
     &              acold(nshg,ndof),           uold(nshg,nsd)
      dimension     iBC(nshg)


      if ( ipred.eq.1) then     ! yn+1_pred=yn
         fct = (gami-one)/gami
         y  = yold
         ac = acold * fct
         if(ideformwall.eq.1) then
            do isd = 1, nsd
               where(btest(iBC,13))
                  u(:,isd) = uold(:,isd) + Delt(itseq)*yold(:,isd) +
     &                 pt5*((gami-two*betai)/gami)*
     &                 Delt(itseq)*Delt(itseq)*acold(:,isd)
               endwhere
            enddo
         endif
      endif
c
      if ( ipred.eq.2) then     ! an+1_pred=0
         y  = yold + (one - gami)/Dtgl * acold
         ac = 0.0
      endif
c
      if(ipred.eq.3 ) then      ! an+1_pred=an
         y  = yold+alfi*Delt(itseq)*acold
         ac = acold
      endif
c
      if ( ipred.eq.4 ) then    ! protect from DC=4 rho=0, same dV
         fct1 = alfi/(one-alfi)
         fct2 = one-almi/gami
         fct3 = almi/gami/alfi*Dtgl
         y    = yold+fct1*(yold-y)
         ac   = acold*fct2+(y-yold)*fct3
      endif
c

      return
      end

c-----------------------------------------------------------------------
c
c    Correct solution at time n+1
c
c-----------------------------------------------------------------------
      subroutine itrCorrect ( y,     ac,   u,   solinc ,iBC)

      include "common.h"

      real*8      y(nshg,ndof),               ac(nshg,ndof),
     &            u(nshg,nsd),                solinc(nshg,4)
      dimension   iBC(nshg)

      fct1 = gami*Delt(itseq)
      fct2 = gami*alfi*Delt(itseq)

      y(:,1:3)  = y(:,1:3)  + fct1 * solinc(:,1:3)
      y(:,4  )  = y(:,4  )  + fct2 * solinc(:,4  )
      ac(:,1:3) = ac(:,1:3) + solinc(:,1:3)
      if(ideformwall.eq.1) then
         do isd = 1,nsd
            where (btest(iBC,13))
               u(:,isd)  = u(:,isd)  +
     &              Delt(itseq)*Delt(itseq)*betai*solinc(:,isd)
            endwhere
         enddo
      endif
c
      return
      end

c-----------------------------------------------------------------------
c
c    Correct solution at time n+1
c
c-----------------------------------------------------------------------
      subroutine itrCorrectSclr ( y,     ac,   solinc )

      include "common.h"

      real*8      y(nshg,ndof),               ac(nshg,ndof),
     &            solinc(nshg)

      fct1 = gami*Delt(itseq)
      is=5+isclr
      y(:,is)  = y(:,is)  + fct1 * solinc(:)
      ac(:,is) = ac(:,is) + solinc(:)
c
      return
      end

c-----------------------------------------------------------------------
c
c    Correct solution at time n+1, protecting against negative values
c
c-----------------------------------------------------------------------
      subroutine itrCorrectSclrPos ( y,     ac,   solinc )

      include "common.h"

      real*8      y(nshg,ndof),               ac(nshg,ndof),
     &            solinc(nshg),               posinc(nshg)

      fct1 = gami*Delt(itseq)
      if(fct1.eq.0) then
         call itrCorrectSclr( y, ac, solinc)
      else
         turbUpdFct = 1.0
         updFct = turbUpdFct
         updFct2 = 1.0
c         if(topHdTscFlag.and.PRM_TCS_PARTIAL_UPDATE)then
c            updFct2 = max(MTH_SQRT_EPS,min(topHdUpdFct,1.0))
c         endif
         fctCr = fct1*updFct*updFct2
         c0 = 0.01
         c1 = -(1.0-c0)/fctCr
         is=5+isclr
c         if(any(updFct*updFct2*solinc(:).lt.c1*y(:,is))) then
c            write(*,*) 'INSTEAD OF GETTING NEGATIVE FIELD'
c            write(*,*) 'BROUGHT FIELD DOWN TO 1 PERCENT'
c            write(*,*) 'FOR SCALAR NUMBER ',isclr
c            write(*,*) '(SEE itrCorrectSclr in itrPC.f)'
c         endif
         posinc = max(updFct*updFct2*solinc,c1*y(:,is))
         y(:,is)  = y(:,is)  + fct1 * posinc(:)
         ac(:,is) = ac(:,is) + posinc(:)
      endif
c
      return
      end


c-----------------------------------------------------------------------
c
c    Compute solution and acceleration at n+alpha
c
c-----------------------------------------------------------------------
      subroutine itrYAlpha ( uold,        yold,        acold,
     &                       u,           y,           ac,
     &                       uAlpha,      yAlpha,      acAlpha )

c      use readarrays       !reads in uold and acold
      include "common.h"

      real*8        yold(nshg,ndof),            acold(nshg,ndof),
     &              y(nshg,ndof),               ac(nshg,ndof),
     &              yAlpha(nshg,ndof),          acAlpha(nshg,ndof),
     &              u(nshg,nsd),                uold(nshg,nsd),
     &              uAlpha(nshg,nsd)

      acAlpha(:,4) = zero  !pressure acceleration is never used but....

      yAlpha (:,1:3) = yold(:,1:3)
     &                  + alfi * (y(:,1:3) - yold(:,1:3))

      acAlpha(:,1:3) = acold(:,1:3)
     &                  + almi * (ac(:,1:3) - acold(:,1:3))

      yAlpha (:,4  ) = y(:,4)

      if(ideformwall.eq.1) uAlpha (:,1:3) = uold(:,1:3)
     &                  + alfi * (u(:,1:3) - uold(:,1:3))

      if(ndof.ge.5) then
c
c  Now take care of temperature, turbulence, what have you
c


         yAlpha (:,5:ndof  ) = yold(:,5:ndof)
     &                       + alfi * (y(:,5:ndof) - yold(:,5:ndof))
         acAlpha(:,5:ndof  ) = acold(:,5:ndof)
     &                       + almi * (ac(:,5:ndof) - acold(:,5:ndof))

      endif
      return
      end

c-----------------------------------------------------------------------
c
c    Update solution at end of time step
c
c-----------------------------------------------------------------------
      subroutine itrUpdate( yold,          acold,        uold,
     &                      y,             ac,           u )

c      use readarrays            !reads in uold and acold

      include "common.h"

      real*8        yold(nshg,ndof),            acold(nshg,ndof),
     &              y(nshg,ndof),               ac(nshg,ndof),
     &              u(nshg,nsd),                uold(nshg,nsd)


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!	Allow for redistance of the level set function, AD 5/8/00
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!      if we are re-distancing the levelset function the we need
!      to update the levelset variable with the corrected value
!      stored in the "next" scalar position.

        if(iLSet .eq.2 .and. isclr .eq. 2) then

           y(:,6) = y(:,7) ! will need to also fix the yold ac etc.?????
     	   ac(:,7)=zero
        endif



      yold  = y
      acold = ac
      if(ideformwall.eq.1)  uold  = u

      return
      end