[Nek5000-users] Adjoint perturbation

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

Hi Nek's,


I am interested into the linear adjoint perturbation. The equations read:

-du/dt = (U.Grad) u - transpose(Grad(U)) u + 1/Re Laplacian(u) - grad(p)

div(u) = 0

Assuming t' = T-t, one can rewrite: du/dt' = (U.Grad) u - transpose(Grad(U))
u + 1/Re Laplacian(u) - grad(p)
Only small modifications to the perturbation mode are necessary, mainly to
evaluate the tranpose gradient terms. Here is what I've done to evaluate
these:

      subroutine conv_adj(conv_x,conv_y,conv_z,u_x,u_y,u_z) ! used to be
conv1n


      include 'SIZE'

      include 'DXYZ'

      include 'INPUT'

      include 'GEOM'

      include 'SOLN'

      include 'TSTEP'


      parameter (lt=lx1*ly1*lz1*lelt)

      real u_x (size(vx,1),size(vx,2),size(vx,3),size(vx,4))

      real u_y (size(vx,1),size(vx,2),size(vx,3),size(vx,4))

      real u_z (size(vx,1),size(vx,2),size(vx,3),size(vx,4))


      real u_xx (size(vx,1),size(vx,2),size(vx,3),size(vx,4))

      real u_xy (size(vx,1),size(vx,2),size(vx,3),size(vx,4))

      real u_xz (size(vx,1),size(vx,2),size(vx,3),size(vx,4))


      real u_yx (size(vy,1),size(vy,2),size(vy,3),size(vy,4))

      real u_yy (size(vy,1),size(vy,2),size(vy,3),size(vy,4))

      real u_yz (size(vy,1),size(vy,2),size(vy,3),size(vy,4))


      real u_zx (size(vz,1),size(vz,2),size(vz,3),size(vz,4))

      real u_zy (size(vz,1),size(vz,2),size(vz,3),size(vz,4))

      real u_zz (size(vz,1),size(vz,2),size(vz,3),size(vz,4))


      real conv_x (size(vx,1),size(vx,2),size(vx,3),size(vx,4))

      real conv_y (size(vy,1),size(vy,2),size(vy,3),size(vy,4))

      real conv_z (size(vz,1),size(vz,2),size(vz,3),size(vz,4))


      *call gradm1(u_xx,u_xy,u_xz,u_x)*

*      call gradm1(u_yx,u_yy,u_yz,u_y)*

*      call gradm1(u_zx,u_zy,u_zz,u_z)*

*
*

*      conv_x = vx*u_xx + vy*u_yx + vz*u_zx*

*      conv_y = vx*u_xy + vy*u_yy + vz*u_zy*

*      conv_z = vx*u_xz + vy*u_yz + vz*u_zz*


      return

      end

However, I'm not sure I've done it properly or not (arrays declaration for
instance). Especially, I had no idea wheter I have to use the gradm1
subroutine or its weak counterpart wgradm1. Finally in perturb.f, I've
slightly change the advap subroutine the following way:

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 conv_adj(ta1,ta2,ta3,tb1,tb2,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

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

c

      endif

c

      return

      end


Any help to spot potential erros or to speed up these operation would be
highly appreciated.

Best regards,
-- 
Jean-Christophe
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20110505/294ce4cc/attachment.html>