[Nek5000-users] Adjoint perturbation
nek5000-users at lists.mcs.anl.gov
nek5000-users at lists.mcs.anl.gov
Tue Jun 21 13:47:59 CDT 2011
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 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> 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
>
More information about the Nek5000-users
mailing list