xref: /phasta/phSolver/incompressible/solfar.f (revision ec121c45c2bc5b2560ecdf76f4d2a36933c1de2c)

      subroutine SolFlow(y,          ac,         u,
     &                   yold,       acold,      uold,
     &                   x,          iBC,
     &                   BC,         res,
     &                   nPermDims,  nTmpDims,  aperm,
     &                   atemp,      iper,
     &                   ilwork,     shp,        shgl,
     &                   shpb,       shglb,      rowp,
     &                   colm,       lhsK,       lhsP,
     &                   solinc,     rerr,       tcorecp,
     &                   GradV,       sumtime,
     &                   svLS_lhs,  svLS_ls,   svLS_nFaces)
c
c----------------------------------------------------------------------
c
c This is the 2nd interface routine to the  linear equation
c solver library that uses the CGP and GMRES methods.
c
c input:
c  y      (nshg,ndof)           : Y-variables at n+alpha_f
c  ac     (nshg,ndof)           : Primvar. accel. variable n+alpha_m
c  yold   (nshg,ndof)           : Y-variables at beginning of step
c  acold   (nshg,ndof)          : Primvar. accel. at beginning of step
c  x      (numnp,nsd)            : node coordinates
c  iBC    (nshg)                : BC codes
c  BC     (nshg,ndofBC)         : BC constraint parameters
c  iper   (nshg)                : periodic nodal information
c
c output:
c  res    (nshg,nflow)           : preconditioned residual
c  y      (nshg,ndof)           : Y-variables at n+alpha_f
c  ac     (nshg,ndof)           : Primvar. accel. variable n+alpha_m
c
c
c The followings are preliminary steps required to use Farzin's
c solver library.  New way of writing has to be used such as
c
c          |  K     G | | du |    | Rmom  |
c          |          | |    | =  |       |
c          | G^t    C | | dp |    | Rcon  |
c
c          |     E    | | dT | =  | Rtemp |
c
c     where
c
c      xKebe : K_ab = dRmom_a/du_b    xTe : E_ab = dRtemp_a/dT_b
c
c              G_ab = dRmom_a/dp_b
c      xGoC  :
c              C_ab = dRcon_a/dp_b
c
c              resf = Rmon Rcon       rest = Rtemp
c
c
c Zdenek Johan,  Winter 1991.  (Fortran 90)
c Juin Kim, Summer 1998. (Incompressible flow solver interface)
c Alberto Figueroa.  CMM-FSI
c----------------------------------------------------------------------
c
      use pointer_data
#ifdef AMG
      use ramg_data
#endif

      include "common.h"
      include "mpif.h"
      include "auxmpi.h"
        include "svLS.h"
c
C
C     Argument variables
C
      INTEGER            npermdims
      INTEGER             ntmpdims
C
C     Local variables
C
      INTEGER              lesid
C
      REAL*8                rdtmp
C
      TYPE(svLS_lhsType), INTENT(INOUT) :: svLS_lhs
      TYPE(svLS_lsType), INTENT(INOUT) ::  svLS_ls

      real*8    y(nshg,ndof),             ac(nshg,ndof),
     &          yold(nshg,ndof),          acold(nshg,ndof),
     &          u(nshg,nsd),              uold(nshg,nsd),
     &          x(numnp,nsd),             BC(nshg,ndofBC),
     &          res(nshg,nflow),          tmpres(nshg,nflow),
     &          flowDiag(nshg,4),
     &          aperm(nshg,nPermDims),    atemp(nshg,nTmpDims),
     &          sclrDiag(nshg,1),
     &          lhsK(9,nnz_tot),          lhsP(4,nnz_tot),
     &          GradV(nshg,nsdsq)
c
      real*8    shp(MAXTOP,maxsh,MAXQPT),
     &          shgl(MAXTOP,nsd,maxsh,MAXQPT),
     &          shpb(MAXTOP,maxsh,MAXQPT),
     &          shglb(MAXTOP,nsd,maxsh,MAXQPT)
c
      integer   usr(100),                 eqnType,temp,
     &          rowp(nshg*nnz),           colm(nshg+1),
     &          iBC(nshg),                ilwork(nlwork),
     &          iper(nshg)
c
      real*8    yAlpha(nshg,ndof),        acAlpha(nshg,ndof),
     &          uAlpha(nshg,nsd),
     &          lesP(nshg,4),             lesQ(nshg,4),
     &          solinc(nshg,ndof),        CGsol(nshg)

      real*8     tcorecp(2)

      real*8    rerr(nshg,10),            rtmp(nshg,4),rtemp

      real*8    msum(4),mval(4),cpusec(10)
      REAL*8 sumtime
      INTEGER dof, svLS_nFaces, i, j, k, l
      INTEGER, ALLOCATABLE :: incL(:)
      REAL*8, ALLOCATABLE :: faceRes(:), Res4(:,:), Val4(:,:)

c
c.... *******************>> Element Data Formation <<******************
c
c
c.... set the parameters for flux and surface tension calculations
c
c
        temp = npro


        idflx = 0
        if(idiff >= 1 )  idflx= (nflow-1) * nsd
        if (isurf == 1) idflx=nflow*nsd
c.... compute solution at n+alpha
c
      call itrYAlpha( uold,    yold,    acold,
     &                u,       y,       ac,
     &                uAlpha,  yAlpha,  acAlpha)

c
c.... form the LHS matrices, the residual vector (at alpha)
c
c      call summary_start()
      impistat=1
      impistat2=1
      telmcp1 = TMRC()
      call ElmGMR (uAlpha,    yAlpha,     acAlpha,    x,
     &             shp,       shgl,       iBC,
     &             BC,        shpb,       shglb,
     &             res,       iper,       ilwork,
     &             rowp,      colm,       lhsK,
     &             lhsP,      rerr,       GradV   )
      telmcp2 = TMRC()
      impistat=0
      impistat2=0
c      call summary_stop()


            tmpres(:,:) = res(:,:)
            iblk = 1
      IF (svLSFlag .EQ. 1) THEN

c####################################################################
!     Here calling svLS

      ALLOCATE(faceRes(svLS_nFaces), incL(svLS_nFaces))
      faceRes=zero  ! function to compute this left out at this time but would be needed to enable adnvanced p vs. Q BC's
      incL = 1
      dof = 4
      IF (.NOT.ALLOCATED(Res4)) THEN
         ALLOCATE (Res4(dof,nshg), Val4(dof*dof,nnz_tot))
      END IF

      DO i=1, nshg
         Res4(1:dof,i) = res(i,1:dof)
      END DO

      DO i=1, nnz_tot
         Val4(1:3,i)   = lhsK(1:3,i)
         Val4(5:7,i)   = lhsK(4:6,i)
         Val4(9:11,i)  = lhsK(7:9,i)
         Val4(13:15,i) = lhsP(1:3,i)
         Val4(16,i)    = lhsP(4,i)
      END DO

      !Val4(4:12:4,:) = -lhsP(1:3,:)^t
      DO i=1, nshg
         Do j=colm(i), colm(i+1) - 1
            k = rowp(j)
            DO l=colm(k), colm(k+1) - 1
               IF (rowp(l) .EQ. i) THEN
                  Val4(4:12:4,l) = -lhsP(1:3,j)
                  EXIT
               END IF
            END DO
         END DO
      END DO
      CALL svLS_SOLVE(svLS_lhs, svLS_ls, dof, Res4, Val4, incL,
     2   faceRes)

      if(myrank.eq.master) write(*,*) 'svLS outer iterations', svLS_ls%RI%itr
      statsflow(1)=1.0*svLS_ls%GM%itr
      statsflow(4)=1.0*svLS_ls%CG%itr
      DO i=1, nshg
         solinc(i,1:dof) = Res4(1:dof,i)
      END DO

c####################################################################
      ELSE
c
c.... lesSolve : main matrix solver
c
      lesId   = numeqns(1)
      eqnType = 1

c      call summary_start()
      impistat=1
      impistat2=1
      tlescp1 = TMRC()
#ifdef AMG
      ! Initial Time Profiling
      call cpu_time(cpusec(1))
      if (irun_amg_prec.gt.0) then
          call ramg_control(colm,rowp,lhsK,lhsP,
     &         ilwork,BC,iBC,iper)
      end if

      call cpu_time(cpusec(6))
      if (irun_amg_prec.gt.0) then
      ramg_flag = 1
      if (irun_amg_prec.eq.2) then ! Some setup work (mode a)
        ramg_window = 1.0
        ramg_redo = 0
      endif
      do while (ramg_flag.le.irun_amg_prec)
      ! if smart solve, possible run solve twice
      ! restart only if meets plateau
#endif

c
c.... setup the linear algebra solver
c
      rtmp = res(:,1:4)
      call usrNew ( usr,        eqnType,          aperm,
     &              atemp,      rtmp,             solinc,
     &              flowDiag,   sclrDiag,         lesP,
     &              lesQ,       iBC,              BC,
     &              iper,       ilwork,           numpe,
     &              nshg,       nshl,             nPermDims,
     &              nTmpDims,   rowp,             colm,
     &              lhsK,       lhsP,             rdtmp,
     &              nnz_tot,    CGsol )
c
c.... solve linear system
c
      call myfLesSolve ( lesId, usr )
#ifdef AMG
      ramg_flag = ramg_flag + 2 ! Default no second run in mode a
      if (irun_amg_prec.eq.3) then
          if (maxIters.gt.int(statsflow(4))) then
          ramg_flag = ramg_flag + 1 ! Default no second run in mode b
          endif
      endif
      enddo
      else
c
c.... setup the linear algebra solver
c
      rtmp = res(:,1:4)
      call usrNew ( usr,        eqnType,          aperm,
     &              atemp,      rtmp,             solinc,
     &              flowDiag,   sclrDiag,         lesP,
     &              lesQ,       iBC,              BC,
     &              iper,       ilwork,           numpe,
     &              nshg,       nshl,             nPermDims,
     &              nTmpDims,   rowp,             colm,
     &              lhsK,       lhsP,             rdtmp,
     &              nnz_tot,    CGsol )

          call myfLesSolve( lesId, usr )
      endif

      call cpu_time(cpusec(3))

      ramg_time(1) = ramg_time(1) + cpusec(3)-cpusec(1)
      ramg_time(7) = ramg_time(7) + cpusec(6)-cpusec(1)

      ! ramg_time: 1 : local total
      !            4 : local VG-cycle
      !            7 : local setup time
      !           11 : Ap-product level 1
      !           12 : Ap-product level >1
      !           13 : Prolongation/Restriction
      !           20 : local boundary MLS time

      if (myrank.eq.master) then
      write(*,*)
      endif
      call ramg_print_time(" == AMG == Total ACUSIM Solver:",
     &                    ramg_time(1))
      call ramg_print_time(" == AMG == Setup: ",ramg_time(7))
      call ramg_print_time(" == AMG == Prec'd cycle: ",ramg_time(4))
      call ramg_print_time(" == AMG == Ap product(level=1): ",
     &                      ramg_time(11))
      call ramg_print_time(" == AMG == Ap product(level>=2): ",
     &                      ramg_time(12))
      call ramg_print_time(" == AMG == Pro/Restr ",
     &                      ramg_time(13))
      call ramg_print_time(" == AMG == Boundary Ap (GS only)",
     &                      ramg_time(20))
      if (myrank.eq.master) then
      write(*,*)
      endif

#endif

      ! End Time profiling output

      call getSol ( usr, solinc )

      if (numpe > 1) then
         call commu ( solinc, ilwork, nflow, 'out')
      endif
      ENDIF ! end of selection between solvers.
      tlescp2 = TMRC()
      impistat=0
      impistat2=0
c      call summary_stop()

      tcorecp(1) = tcorecp(1) + telmcp2-telmcp1 ! elem. formation
      tcorecp(2) = tcorecp(2) + tlescp2-tlescp1 ! linear alg. solution
      call rstatic (res, y, solinc) ! output flow stats
c
c.... end
c
      return
      end

      subroutine SolSclr(y,          ac,         u,
     &                   yold,       acold,      uold,
     &                   x,          iBC,
     &                   BC,         nPermDimsS,  nTmpDimsS,
     &                   apermS,     atempS,     iper,
     &                   ilwork,     shp,        shgl,
     &                   shpb,       shglb,      rowp,
     &                   colm,       lhsS,       solinc,
     &                   tcorecpscal,
     &                   svLS_lhs,  svLS_ls,   svLS_nFaces)
c
c----------------------------------------------------------------------
c
c This is the 2nd interface routine to the Farzin's linear equation
c solver library.
c
c input:
c  y      (nshg,ndof)           : Y-variables at n+alpha_f
c  ac     (nshg,ndof)           : Primvar. accel. variable n+alpha_m
c  yold   (nshg,ndof)           : Y-variables at beginning of step
c  x      (numnp,nsd)            : node coordinates
c  iBC    (nshg)                : BC codes
c  BC     (nshg,ndofBC)         : BC constraint parameters
c  iper   (nshg)                : periodic nodal information
c
c output:
c  y      (nshg,ndof)           : Y-variables at n+alpha_f
c  ac     (nshg,ndof)           : Primvar. accel. variable n+alpha_m
c
c
c The followings are preliminary steps required to use Farzin's
c solver library.  New way of writing has to be used such as
c
c          |     E    | | dS | =  | RScal |
c
c----------------------------------------------------------------------
c
      use pointer_data

      include "common.h"
      include "mpif.h"
      include "auxmpi.h"
        include "svLS.h"
c
      real*8    y(nshg,ndof),             ac(nshg,ndof),
     &          yold(nshg,ndof),          acold(nshg,ndof),
     &          u(nshg,nsd),              uold(nshg,nsd),
     &          x(numnp,nsd),             BC(nshg,ndofBC),
     &          res(nshg,1),
     &          flowDiag(nshg,4),
     &          sclrDiag(nshg,1),           lhsS(nnz_tot),
     &          apermS(nshg,nPermDimsS),  atempS(nshg,nTmpDimsS)

c
      real*8    shp(MAXTOP,maxsh,MAXQPT),
     &          shgl(MAXTOP,nsd,maxsh,MAXQPT),
     &          shpb(MAXTOP,maxsh,MAXQPT),
     &          shglb(MAXTOP,nsd,maxsh,MAXQPT)
c
      integer   usr(100),                 eqnType,
     &          rowp(nshg*nnz),           colm(nshg+1),
     &          iBC(nshg),                ilwork(nlwork),
     &          iper(nshg)
c
      real*8    yAlpha(nshg,ndof),        acAlpha(nshg,ndof),
     &          uAlpha(nshg,nsd),
     &          lesP(nshg,1),               lesQ(nshg,1),
     &          solinc(nshg,1),           CGsol(nshg),
     &          tcorecpscal(2)
      TYPE(svLS_lhsType), INTENT(INOUT) :: svLS_lhs
      TYPE(svLS_lsType), INTENT(INOUT) ::  svLS_ls
      REAL*8 sumtime
      INTEGER dof, svLS_nFaces, i, j, k, l
      INTEGER, ALLOCATABLE :: incL(:)
      REAL*8, ALLOCATABLE :: faceRes(:), Res1(:,:), Val1(:,:)

c
c.... *******************>> Element Data Formation <<******************
c
c.... compute solution at n+alpha
c
      call itrYAlpha( uold,    yold,    acold,
     &                u,       y,       ac,
     &                uAlpha,  yAlpha,  acAlpha)
c
c.... form the LHS matrices, the residual vector (at alpha)
c
      impistat=2
      impistat2=1
      telmcp1 = TMRC()
      call ElmGMRSclr(yAlpha,acAlpha,    x,
     &             shp,       shgl,       iBC,
     &             BC,        shpb,       shglb,
     &             res,       iper,       ilwork,
     &             rowp,      colm,       lhsS   )
      telmcp2 = TMRC()
      impistat=0
      impistat2=0
      statssclr(1)=0
      IF (svLSFlag .EQ. 1) THEN

c####################################################################
!     Here calling svLS

      ALLOCATE(faceRes(svLS_nFaces), incL(svLS_nFaces))
      faceRes=zero
      incL = 1
      dof = 1
      IF (.NOT.ALLOCATED(Res1)) THEN
         ALLOCATE (Res1(dof,nshg), Val1(dof*dof,nnz_tot))
      END IF

      DO i=1, nshg
         Res1(1,i) = res(i,1)
      END DO

      DO i=1, nnz_tot
         Val1(1,i)    = lhsS(i)
      END DO

      CALL svLS_SOLVE(svLS_lhs, svLS_ls, dof, Res1, Val1, incL,
     2   faceRes)
      statssclr(1)=1.0*svLS_ls%RI%itr
      DO i=1, nshg
         solinc(i,1) = Res1(1,i)
      END DO

c####################################################################
      ELSE
c
c.... lesSolve : main matrix solver
c
      lesId   = numeqns(1+nsolt+isclr)
      eqnType = 2
c
c.... setup the linear algebra solver
c

      impistat=2
      impistat2=1
      tlescp1 = TMRC()
      call usrNew ( usr,        eqnType,          apermS,
     &              atempS,     res,              solinc,
     &              flowDiag,   sclrDiag,         lesP,
     &              lesQ,       iBC,              BC,
     &              iper,       ilwork,           numpe,
     &              nshg,       nshl,             nPermDimsS,
     &              nTmpDimsS,  rowp,             colm,
     &              rlhsK,      rlhsP,            lhsS,
     &              nnz_tot,    CGsol )
c
c.... solve linear system
c
      call myfLesSolve ( lesId, usr )
      call getSol ( usr, solinc )

      if (numpe > 1) then
         call commu ( solinc, ilwork, 1, 'out')
      endif
      ENDIF ! decision between solvers
      tlescp2 = TMRC()
      impistat=0
      impistat2=0

      tcorecpscal(1) = tcorecpscal(1) + telmcp2-telmcp1 ! elem. formation
      tcorecpscal(2) = tcorecpscal(2) + tlescp2-tlescp1 ! linear alg. solution

      nsolsc=5+isclr
      call rstaticSclr (res, y, solinc, nsolsc) ! output scalar stats
c
c.... end
c
      return
      end