conduct.F90

Go to the documentation of this file.

!-----------------------------------------------------------------------
subroutine cdscal (igeom)
  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1, lelt, nx1, ny1, nz1, nfield, nid
  use helmholtz, only : hsolve, approx_space
  use input, only : ifmodel, ifkeps, ifaxis, ifaziv, iftran, iftmsh, ifprint
  use geom, only : bintm1, binvm1
  use soln, only : t, bq, tmask, tmult
  use tstep, only : nelfld, ifield, nmxnl, imesh, tolht, nmxh
  use ctimer, only : nheat2, theat2, dnekclock
  implicit none

  integer, intent(in) :: igeom

  LOGICAL ::          IFCONV

  real(DP), allocatable :: TA(:,:,:,:), TB(:,:,:,:)
  real(DP), allocatable :: H1(:,:,:,:), H2(:,:,:,:)

  real(DP) :: etime

  type(approx_space), save :: apx
  character(4) ::     name4

  integer :: nel, ntot, nfldt, if1, isd, iter, intype

  nel    = nelfld(ifield)
  ntot   = nx1*ny1*nz1*nel

  if (igeom == 1) then   ! geometry at t^{n-1}
      call makeq
      call lagscal
  else                   ! geometry at t^n

      IF (ifprint) THEN
          IF (ifmodel .AND. ifkeps) THEN
              nfldt = nfield - 1
              IF (ifield == nfldt .AND. nid == 0) THEN
                  WRITE (6,*) ' Turbulence Model - k/epsilon solution'
              ENDIF
          ELSE
              IF (ifield == 2 .AND. nid == 0) THEN
                  WRITE (6,*) ' Temperature/Passive scalar solution'
              ENDIF
          ENDIF
      ENDIF
      if1=ifield-1
      write(name4,1) if1-1
      1 format('PS',i2)
      if(ifield == 2) write(name4,'(A4)') 'TEMP'

  
  !        New geometry
  
      isd = 1
      if (ifaxis .AND. ifaziv .AND. ifield == 2) isd = 2
  !        if (ifaxis.and.ifmhd) isd = 2 !This is a problem if T is to be T!
      etime = dnekclock()
      nheat2 = nheat2 + 1
      allocate(ta(lx1,ly1,lz1,lelt), tb(lx1,ly1,lz1,lelt))
      allocate(h1(lx1,ly1,lz1,lelt), h2(lx1,ly1,lz1,lelt))
      do iter=1,nmxnl ! iterate for nonlin. prob. (e.g. radiation b.c.)

          intype = 0
          IF (iftran) intype = -1
          CALL sethlm(h1,h2,intype)
          CALL bcneusc(ta,-1)
          h2 = h2 + ta
          CALL bcdirsc(t(1,1,1,1,ifield-1))
          etime = etime - dnekclock()
          CALL axhelm(ta,t(1,1,1,1,ifield-1),h1,h2,imesh,isd)
          etime = etime + dnekclock()
          tb = bq(:,:,:,:,ifield-1) - ta
          CALL bcneusc(ta,1)
          tb = tb + ta

      !        CALL HMHOLTZ (name4,TA,TB,H1,H2
      !    $                 ,TMASK(1,1,1,1,IFIELD-1)
      !    $                 ,TMULT(1,1,1,1,IFIELD-1)
      !    $                 ,IMESH,TOLHT(IFIELD),NMXH,isd)

          etime = etime - dnekclock()
          if(iftmsh(ifield)) then
              call hsolve(name4,ta,tb,h1,h2 &
              ,tmask(1,1,1,1,ifield-1) &
              ,tmult(1,1,1,1,ifield-1) &
              ,imesh,tolht(ifield),nmxh,1 &
              ,apx,bintm1)
          else
              call hsolve(name4,ta,tb,h1,h2 &
              ,tmask(1,1,1,1,ifield-1) &
              ,tmult(1,1,1,1,ifield-1) &
              ,imesh,tolht(ifield),nmxh,1 &
              ,apx,binvm1)
          endif
          etime = etime + dnekclock()

          t(:,:,:,:,ifield-1) = t(:,:,:,:,ifield-1) + ta

          call cvgnlps(ifconv) ! Check convergence for nonlinear problem
          if (ifconv) goto 2000

      !        Radiation case, smooth convergence, avoid flip-flop (ER).
          t(:,:,:,:,ifield-1) = t(:,:,:,:,ifield-1) - 0.5_dp*ta

      END DO
      2000 CONTINUE
      deallocate(h1,h2,tb)
      CALL bcneusc(ta,1)
      bq(:,:,:,:,ifield-1) = bq(:,:,:,:,ifield-1) + ta ! no idea why... pf
      deallocate(ta)
      theat2 = theat2 + (dnekclock() - etime)
  endif

  return
end subroutine cdscal

!-----------------------------------------------------------------------
subroutine makeuq
  use kinds, only : dp
  use size_m, only : nx1, ny1, nz1
  use geom, only : bm1
  use soln, only : bq
  use tstep, only : time, dt, ifield, nelfld
  implicit none

  integer :: ntot

  ntot = nx1*ny1*nz1*nelfld(ifield)

  time = time-dt        ! Set time to t^n-1 for user function

  bq(:,:,:,:,ifield-1) = 0._dp
  call setqvol( bq(1,1,1,1,ifield-1)          )
  bq(:,:,:,:,ifield-1) = bq(:,:,:,:,ifield-1) * bm1

  time = time+dt        ! Restore time

  return
end subroutine makeuq

!-----------------------------------------------------------------------
subroutine setqvol(bql)
  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1, lelt, nx1, ny1, nz1
  use input, only : igroup, matype, cpgrp
  use tstep, only : ifield, nelfld
  implicit none

  real(DP) :: bql(lx1*ly1*lz1,lelt)
  real(DP) :: cqvol
  integer :: nel, nxyz1, ntot1, iel, igrp

#ifndef MOAB
  nel   = nelfld(ifield)
  nxyz1 = nx1*ny1*nz1
  ntot1 = nxyz1*nel

  do iel=1,nel
      igrp = igroup(iel)
      if (matype(igrp,ifield) == 1) then ! constant source within a group
          cqvol = cpgrp(igrp,ifield,3)
          bql(:,iel) = cqvol
      else  !  pff 2/6/96 ............ default is to look at userq
          call nekuq(bql,iel)
      endif
  enddo

! 101 FORMAT(' Wrong material type (',I3,') for group',I3,', field',I2
!    $    ,/,' Aborting in SETQVOL.')
#else
! pulling in temperature right now, since we dont have anything else
  call userq2(bql)
#endif

  return
end subroutine setqvol

!------------------------------------------------------------------
!------------------------------------------------------------------
subroutine nekuq (bql,iel)
  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1, lelt, nx1, ny1, nz1
  use nekuse, only : qvol
  use parallel, only : lglel
  implicit none

  real(DP) :: bql(lx1,ly1,lz1,lelt)
  integer :: iel

  integer :: ielg, k, j, i

  ielg = lglel(iel)
  do k=1,nz1
      do j=1,ny1
          do i=1,nx1
              call nekasgn(i,j,k,iel)
              qvol = 0.0
              call userq(i,j,k,ielg)
              bql(i,j,k,iel) = qvol
          enddo
      enddo
  END DO

  return
end subroutine nekuq

!-----------------------------------------------------------------------
!---------------------------------------------------------------
subroutine convab()
  use kinds, only : dp
  use size_m, only : nx1, ny1, nz1, lx1, ly1, lz1, lelt
  use geom, only : bm1
  use soln, only : t, vtrans, bq
  use tstep, only : ifield, nelfld
  implicit none

  real(DP), allocatable :: TA (:,:,:,:)
  integer :: nel, ntot1

  allocate(ta(lx1,ly1,lz1,lelt))

  nel = nelfld(ifield)
  ntot1 = nx1*ny1*nz1*nel
  CALL convop(ta,t(1,1,1,1,ifield-1))
  bq(:,:,:,:,ifield-1) = bq(:,:,:,:,ifield-1) - bm1*ta *vtrans(:,:,:,:,ifield)

  return
end subroutine convab

!-----------------------------------------------------------------------
subroutine makeabq
  use kinds, only : dp
  use size_m, only : nx1, ny1, nz1, lx1, ly1, lz1, lelt
  use soln, only : vgradt1, vgradt2, bq
  use tstep, only : ab, ifield, nelfld
  implicit none

  real(DP), allocatable :: TA (:,:,:,:)
  real(DP) :: ab0, ab1, ab2
  integer :: nel, ntot1

  allocate(ta(lx1,ly1,lz1,lelt))

  ab0   = ab(1)
  ab1   = ab(2)
  ab2   = ab(3)
  nel   = nelfld(ifield)
  ntot1 = nx1*ny1*nz1*nel

  ta = ab1*vgradt1(:,:,:,:,ifield-1) + ab2*vgradt2(:,:,:,:,ifield-1)
  vgradt2(:,:,:,:,ifield-1) = vgradt1(:,:,:,:,ifield-1)
  vgradt1(:,:,:,:,ifield-1) = bq(:,:,:,:,ifield-1)
  bq(:,:,:,:,ifield-1) = ab0*bq(:,:,:,:,ifield-1) + ta

  return
end subroutine makeabq

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
subroutine makebdq()
  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1, lelt
  use size_m, only : nx1, ny1, nz1
  use geom, only : ifgeom
  use geom, only : bm1, bm1lag
  use soln, only : vtrans, t, tlag, bq
  use tstep, only : ifield, nelfld, dt, bd, nbd
  implicit none

  real(DP), allocatable :: TB(:,:,:,:)

  integer :: nel, ntot1, ilag

  allocate(tb(lx1,ly1,lz1,lelt))

  nel   = nelfld(ifield)
  ntot1 = nx1*ny1*nz1*nel

  tb = bd(2) * bm1 * t(:,:,:,:,ifield-1)

  DO ilag=2,nbd
      IF (ifgeom) THEN
          tb = tb + bd(ilag+1) * bm1lag(:,:,:,:,ilag-1) * tlag(:,:,:,:,ilag-1, ifield-1)
      ELSE
          tb = tb + bd(ilag+1) * bm1 * tlag(:,:,:,:,ilag-1,ifield-1)
      ENDIF
  END DO

  tb = (1./dt) * tb * vtrans(:,:,:,:,ifield)
  bq(:,:,:,:,ifield-1) = bq(:,:,:,:,ifield-1) + tb

  return
end subroutine makebdq

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
subroutine lagscal
  use size_m, only : nx1, ny1, nz1
  use soln, only : t, tlag
  use tstep, only : ifield, nelfld, nbdinp
  implicit none

  integer :: ntot1, ilag
  ntot1 = nx1*ny1*nz1*nelfld(ifield)

  DO 100 ilag=nbdinp-1,2,-1
      CALL copy(tlag(1,1,1,1,ilag  ,ifield-1), &
      tlag(1,1,1,1,ilag-1,ifield-1),ntot1)
  100 END DO

  CALL copy(tlag(1,1,1,1,1,ifield-1),t(1,1,1,1,ifield-1),ntot1)

  return
  end subroutine lagscal

!-----------------------------------------------------------------------