[Nek5000-users] Adjoint perturbation

nek5000-users at lists.mcs.anl.gov nek5000-users at lists.mcs.anl.gov
Tue Jun 21 17:41:21 CDT 2011


Ah alright, I forgot to include 'ADJOINT'.
Fair enough then, I'll try again tomorrow morning.

Thans a lots.

On 21 June 2011 20:47, <nek5000-users at lists.mcs.anl.gov> wrote:

> Hi Jean-Christophe,
>
>
> It looks like the only place in the code where ifadj shows up is in
> perturb.f:107
>
>      if (ifnav.and.(.not.ifchar).and.(     ifadj))call advabp_adjoint
>
> so you can set ifadj=.true. somewhere in userchk or usrdat -- just make
> that you also add
>
>      include 'ADJOINT'
>
> there which is not part of 'TOTAL'
>
> Best,
> Aleks
>
>
>
>
> On Tue, 21 Jun 2011, nek5000-users at lists.mcs.anl.**gov<nek5000-users at lists.mcs.anl.gov>wrote:
>
>  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<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<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 <Nek5000-users at lists.mcs.anl.gov>
>>>> https://lists.mcs.anl.gov/**mailman/listinfo/nek5000-users<https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users>
>>>>
>>>>
>>>
>>>
>>> --
>>> Jean-Christophe
>>>
>>>
>>
>>
>> --
>> Jean-Christophe
>>
>>  ______________________________**_________________
> Nek5000-users mailing list
> Nek5000-users at lists.mcs.anl.**gov <Nek5000-users at lists.mcs.anl.gov>
> https://lists.mcs.anl.gov/**mailman/listinfo/nek5000-users<https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users>
>



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


More information about the Nek5000-users mailing list