convect.F90

Go to the documentation of this file.

!-----------------------------------------------------------------------
!    Stability limits:
!    AB3:    .7236                     w/safety (1.2):   .603
!    RK3:    1.73   (sqrt 3)           w/safety (1.2):   1.44
!    RK4:    2.828                     w/safety (1.2):   2.36
!    SEM Safety factor:  1.52 for N=3
!                     <  1.20 for N=16
!                     ~  1.16 for N=256
!-----------------------------------------------------------------------
#define INLINE_INTP
subroutine setup_convect(igeom)
  use dealias, only : vxd, vyd, vzd
  use input, only : param, ifchar, ifcons, ifpert
  use soln, only : vx, vy, vz
  implicit none

  integer, intent(in) :: igeom

  if (igeom == 1) return
  if (param(99) < 0) return ! no dealiasing

  if (ifchar) then
    write(*,*) "Oops: ifchar"
#if 0
    nelc = nelv
    if (ifmhd) nelc = max(nelv,nelfld(ifldmhd))
    if (ifmhd) call exitti('no characteristics for mhd yet$',istep)

    ifnew = .true. 
    if (igeom > 2) ifnew = .false. 
    call set_conv_char(ct_vx,c_vx,vx,vy,vz,nelc,time,ifnew)
#endif

  else

    if ( .NOT. ifpert) then
        if (ifcons) then
!max                call set_convect_cons (vxd,vyd,vzd,vx,vy,vz)
        else
            call set_convect_new(vxd,vyd,vzd,vx,vy,vz)
        endif
    endif

  endif

  return
end subroutine setup_convect
!-----------------------------------------------------------------------
subroutine convect_new(bdu,u,ifuf,cx,cy,cz,ifcf)
  use kinds, only : dp
  use size_m, only : nelv
  use size_m, only : lxd, lyd, lzd, ldim
  use size_m, only : nx1, ny1, nz1, nxd, nyd, nzd
  use input, only : if3d
  use ctimer, only : tscn, dnekclock
  use interp, only : jgl, jgt, dgl, dgt
  use interp, only : get_int_ptr, get_dgl_ptr
  implicit none

  real(DP) :: bdu(*),u(*),cx(*),cy(*),cz(*)
  logical :: ifuf,ifcf            ! u and/or c already on fine mesh?

  integer, parameter :: ltd=lxd*lyd*lzd
  real(DP) :: ur(ltd), us(ltd), ut(ltd), tr(ltd,3), uf(ltd)

  integer :: e, iu, ic, ib, i, iptr, iptr2
  integer :: nxyz1, nxyzd, nxyzu, nxyzc
  real(DP) :: etime
  real(DP) :: w((2*lxd)**ldim,2)
  integer, parameter :: ldw = 2*(2*lxd)**ldim

  etime = dnekclock()
!max  call set_dealias_rx()

  nxyz1 = nx1*ny1*nz1
  nxyzd = nxd*nyd*nzd

  nxyzu = nxyz1
  if (ifuf) nxyzu = nxyzd

  nxyzc = nxyz1
  if (ifcf) nxyzc = nxyzd

  iu = 1    ! pointer to scalar field u
  ic = 1    ! pointer to vector field C
  ib = 1    ! pointer to scalar field Bdu

#ifdef INLINE_INTP
  call get_int_ptr(iptr,  nx1,nxd)
  call get_dgl_ptr(iptr2, nxd,nxd)
#endif

  do e=1,nelv

      if (ifcf) then
          call copy(tr(1,1),cx(ic),nxyzd)  ! already in rst form
          call copy(tr(1,2),cy(ic),nxyzd)
          if (if3d) call copy(tr(1,3),cz(ic),nxyzd)

      else  ! map coarse velocity to fine mesh (C-->F)
        write(*,*) "Oops: ifcf"
#if 0
          call intp_rstd(fx,cx(ic),nx1,nxd,if3d,0) ! 0 --> forward
          call intp_rstd(fy,cy(ic),nx1,nxd,if3d,0) ! 0 --> forward
          if (if3d) call intp_rstd(fz,cz(ic),nx1,nxd,if3d,0) ! 0 --> forward

          if (if3d) then  ! Convert convector F to r-s-t coordinates

              do i=1,nxyzd
                  tr(i,1)=rx(i,1,e)*fx(i)+rx(i,2,e)*fy(i)+rx(i,3,e)*fz(i)
                  tr(i,2)=rx(i,4,e)*fx(i)+rx(i,5,e)*fy(i)+rx(i,6,e)*fz(i)
                  tr(i,3)=rx(i,7,e)*fx(i)+rx(i,8,e)*fy(i)+rx(i,9,e)*fz(i)
              enddo

          else

              do i=1,nxyzd
                  tr(i,1)=rx(i,1,e)*fx(i)+rx(i,2,e)*fy(i)
                  tr(i,2)=rx(i,3,e)*fx(i)+rx(i,4,e)*fy(i)
              enddo

          endif
#endif
      endif

!      etime = etime - dnekclock()
#ifndef INLINE_INTP
      if (ifuf) then
          call grad_rst(ur,us,ut,u(iu),nxd,if3d)
      else
          call intp_rstd(uf,u(iu),nx1,nxd,if3d,0) ! 0 --> forward
          call grad_rst(ur,us,ut,uf,nxd,if3d)
      endif
#else
      if (ifuf) then
          call local_grad3(ur,us,ut,u(iu),nxd-1,1,dgl(iptr2),dgt(iptr2))
      else
          call specmpn(uf,nxd,u(iu),nx1,jgl(iptr),jgt(iptr),if3d,w,ldw)
          call local_grad3(ur,us,ut,uf,nxd-1,1,dgl(iptr2),dgt(iptr2))
      endif
#endif
!      etime = etime + dnekclock()

      if (if3d) then
          do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
              uf(i) = tr(i,1)*ur(i)+tr(i,2)*us(i)+tr(i,3)*ut(i)
          enddo
      else
          do i=1,nxyzd ! mass matrix included, per DFM (4.8.5)
              uf(i) = tr(i,1)*ur(i)+tr(i,2)*us(i)
          enddo
      endif
!      etime = etime - dnekclock()
#ifndef INLINE_INTP
      call intp_rstd(bdu(ib),uf,nx1,nxd,if3d,1) ! Project back to coarse
#else
      call specmpn(bdu(ib),nx1,uf,nxd,jgt(iptr),jgl(iptr),if3d,w,ldw)
#endif
!      etime = etime + dnekclock()

      ic = ic + nxyzc
      iu = iu + nxyzu
      ib = ib + nxyz1

  enddo
  tscn = tscn + (dnekclock() - etime)

  return
  end subroutine convect_new

!-----------------------------------------------------------------------
subroutine set_convect_new(cr,cs,ct,ux,uy,uz)
  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1, lxd, lyd, lzd, ldim
  use size_m, only : nx1, ny1, nz1, nxd, nyd, nzd, nelv
  use geom, only : rx
  use input, only : if3d
  use mesh, only : if_ortho
  use ctimer, only : tscn, nscn, dnekclock
  use interp, only : jgl, jgt
  use interp, only : get_int_ptr

  implicit none

  integer, parameter :: lxy=lx1*ly1*lz1, ltd=lxd*lyd*lzd

  real(DP) :: cr(ltd,*),cs(ltd,*),ct(ltd,*)
  real(DP) :: ux(lxy,*),uy(lxy,*),uz(lxy,*)

  real(DP) :: fx(ltd), fy(ltd), fz(ltd)!, ur, us, ut, tr, uf
  real(DP) :: w((2*lxd)**ldim,2)
  integer, parameter :: ldw = 2*(2*lxd)**ldim

  real(DP) :: etime

  integer :: e, nxyz1, nxyzd, ic, i, j, iptr
  nscn = nscn + 1
  etime = dnekclock()
  call set_dealias_rx()

  nxyz1 = nx1*ny1*nz1
  nxyzd = nxd*nyd*nzd

  ic = 1    ! pointer to vector field C
#ifdef INLINE_INTP
  call get_int_ptr(iptr, nx1, nxd) 
#endif

  do e=1,nelv

  !        Map coarse velocity to fine mesh (C-->F)

      !call specmpn(fx,nxd,ux(1,e),nx1,jgl(i),jgt(i),if3d,w,ldw)
      !call specmpn(fy,nxd,uy(1,e),nx1,jgl(i),jgt(i),if3d,w,ldw)
      !call specmpn(fz,nxd,uz(1,e),nx1,jgl(i),jgt(i),if3d,w,ldw)
!      etime = etime - dnekclock()
#ifdef INLINE_INTP
      call specmpn(fx,nxd,ux(1,e),nx1,jgl(iptr),jgt(iptr),if3d,w,ldw)
      call specmpn(fy,nxd,uy(1,e),nx1,jgl(iptr),jgt(iptr),if3d,w,ldw)
      if (if3d) call specmpn(fz,nxd,uz(1,e),nx1,jgl(iptr),jgt(iptr),if3d,w,ldw)
#else
      call intp_rstd(fx,ux(1,e),nx1,nxd,if3d,0) ! 0 --> forward
      call intp_rstd(fy,uy(1,e),nx1,nxd,if3d,0) ! 0 --> forward
      if (if3d) call intp_rstd(fz,uz(1,e),nx1,nxd,if3d,0) ! 0 --> forward
#endif
!      etime = etime + dnekclock()

  !        Convert convector F to r-s-t coordinates

      if (if3d) then

          if (if_ortho) then
            do i=1,nxyzd
                cr(i,e)=rx(i,1,e)*fx(i)
                cs(i,e)=rx(i,2,e)*fy(i)
                ct(i,e)=rx(i,3,e)*fz(i)
            enddo
          else
            do i=1,nxyzd
                cr(i,e)=rx(i,1,e)*fx(i)+rx(i,2,e)*fy(i)+rx(i,3,e)*fz(i)
                cs(i,e)=rx(i,4,e)*fx(i)+rx(i,5,e)*fy(i)+rx(i,6,e)*fz(i)
                ct(i,e)=rx(i,7,e)*fx(i)+rx(i,8,e)*fy(i)+rx(i,9,e)*fz(i)
            enddo
          endif
      else
#if 0
          do i=1,nxyzd
              cr(i,e)=rx(i,1,e)*fx(i)+rx(i,2,e)*fy(i)
              cs(i,e)=rx(i,3,e)*fx(i)+rx(i,4,e)*fy(i)
          enddo
#endif
      endif
  enddo
  tscn = tscn + (dnekclock() - etime)

  return
end subroutine set_convect_new
!-----------------------------------------------------------------------