[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