[Nek5000-users] Adjoint perturbation

Tue Jun 21 10:45:36 CDT 2011

Hi Nek's,

I have seen the adjoint perturbation mode has been added to the current.
Where/how do I have to specify ifadj = .true. if I want to use it? Is there
some parameters in the .rea file like for direct perturbation mode?

Sincerely,

On 9 June 2011 11:50, Jean-Christophe Loiseau <loiseau.jc at gmail.com> wrote:

> Hi Neks,
>
> I was just wondering if the adjoint linear solver has already been released
> in the current repo or not?
>
> Regards,
>
>
> On 6 May 2011 15:49, <nek5000-users at lists.mcs.anl.gov> wrote:
>
>>
>>
>> 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
>>>
>>>  _______________________________________________
>> Nek5000-users mailing list
>> Nek5000-users at lists.mcs.anl.gov
>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>>
>
>
>
> --
> Jean-Christophe
>

-- 
Jean-Christophe
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20110621/ae014e30/attachment.html>