[Nek5000-users] Adjoint perturbation

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

Dear Jean-Christophe,

My apologies for the delay in getting back to you - we've been
saturated with proposal writing on our end.

Your approach looks reasonable -- though I would avoid at all
costs the f90 type declarations and stick with the std f77-based
approach to be consistent with the code.

I'll look into your proposed approach.  I should point out also
that we've recently had an adjoint solver developed by some 
colleagues that is ready to go into the svn repo.  I've just
been delayed in getting it in because of travel and proposals.
My hope is that this will be in the repo in about 10-15 days.

Regarding the proper choice of gradm1 vs wgradm1, the question
ultimately comes down to the degree of differentiability of the
input --- functions in nek are only C0 and hence one-time differentiable.
This applies also to the test functions.  Thus, the std. practice
is to move a derivative from the solution to the test functions (i.e.,
to use the weak derivative) whenever the solution is to be differentiated
twice.

Regards,

Paul




On Thu, 5 May 2011, nek5000-users at lists.mcs.anl.gov wrote:

> 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
>