[Nek5000-users] Adjoint perturbation

Thu Jun 9 04:50:08 CDT 2011

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20110609/846a1aa6/attachment.html>