[Nek5000-users] Weakly non linear Navier-Stokes

[nek5000-users at lists.mcs.anl.gov]

Hi Neks !

We are interested in studying the weakly non linear dynamics of a 
perturbation. We are therefore using Nek in perturbation mode, and we 
have implemented the usually omitted 2nd order non linear perturbation 
convection terms in advabp() in perturb.f as follows :

      subroutine advabp
C
C     Eulerian scheme, add convection term to forcing function
C     at current time step.
C
       include 'SIZE'
       include 'INPUT'
       include 'SOLN'
       include 'MASS'
       include 'TSTEP'
C
       COMMON /SCRNS/ TA1 (LX1*LY1*LZ1*LELV)
      $ ,             TA2 (LX1*LY1*LZ1*LELV)
      $ ,             TA3 (LX1*LY1*LZ1*LELV)
      $ ,             TB1 (LX1*LY1*LZ1*LELV)
      $ ,             TB2 (LX1*LY1*LZ1*LELV)
      $ ,             TB3 (LX1*LY1*LZ1*LELV)
C
       ntot1 = nx1*ny1*nz1*nelv
       ntot2 = nx2*ny2*nz2*nelv
c
       if (if3d) then
          call opcopy  (tb1,tb2,tb3,vx,vy,vz)                   ! Save velocity
          call opcopy  (vx,vy,vz,vxp(1,jp),vyp(1,jp),vzp(1,jp)) ! U <-- dU
          call convop  (ta1,tb1)                                ! du.grad U
          call convop  (ta2,tb2)
          call convop  (ta3,tb3)
          call opcopy  (vx,vy,vz,tb1,tb2,tb3)  ! Restore velocity
c
          do i=1,ntot1
             tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
             bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
             bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
             bfzp(i,jp) = bfzp(i,jp)-tmp*ta3(i)
          enddo
c
          call convop  (ta1,vxp(1,jp))       !  U.grad dU
          call convop  (ta2,vyp(1,jp))
          call convop  (ta3,vzp(1,jp))
c
          do i=1,ntot1
             tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
             bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
             bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
             bfzp(i,jp) = bfzp(i,jp)-tmp*ta3(i)
          enddo

      ! ADD NON LINEAR TERM
          call opcopy  (tb1,tb2,tb3,vx,vy,vz)
          call opcopy  (vx,vy,vz,vxp(1,jp),vyp(1,jp),vzp(1,jp))
          call convop  (ta1,vxp(1,jp))
          call convop  (ta2,vyp(1,jp))
          call convop  (ta3,vzp(1,jp))
          call opcopy  (vx,vy,vz,tb1,tb2,tb3)
          
          do i=1,ntot1
             tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
             bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
             bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
             bfzp(i,jp) = bfzp(i,jp)-tmp*ta3(i)
          enddo

c
       else  ! 2D
c
          call opcopy  (tb1,tb2,tb3,vx,vy,vz)                   ! Save velocity
          call opcopy  (vx,vy,vz,vxp(1,jp),vyp(1,jp),vzp(1,jp)) ! U <-- dU
          call convop  (ta1,tb1)                                ! du.grad U
          call convop  (ta2,tb2)
          call opcopy  (vx,vy,vz,tb1,tb2,tb3)  ! Restore velocity
c
          do i=1,ntot1
             tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
             bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
             bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
          enddo
c
          call convop  (ta1,vxp(1,jp))       !  U.grad dU
          call convop  (ta2,vyp(1,jp))
c
          do i=1,ntot1
             tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
             bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
             bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
          enddo

      ! ADD NON LINEAR TERM
          call opcopy  (tb1,tb2,tb3,vx,vy,vz)
          call opcopy  (vx,vy,vz,vxp(1,jp),vyp(1,jp),vzp(1,jp))
          call convop  (ta1,vxp(1,jp))
          call convop  (ta2,vyp(1,jp))
          call opcopy  (vx,vy,vz,tb1,tb2,tb3)
          
          do i=1,ntot1
             tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,ifield)
             bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)
             bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)
          enddo
c
       endif
c
       return
       end


Is this the correct implementation ? Are there any other modifications 
to be made ? An alternative would obviously be to add the non linear 
perturbation convection term in USERF, would that be a safer solution ?

Any help would be greatly appreciated !

Best,

Holly
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20140225/b3a36ea5/attachment.html>