[Nek5000-users] Adjoint perturbation
nek5000-users at lists.mcs.anl.gov
nek5000-users at lists.mcs.anl.gov
Thu May 5 06:11:41 CDT 2011
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>
More information about the Nek5000-users
mailing list