Ah alright, I forgot to include 'ADJOINT'.<br>Fair enough then, I'll try again tomorrow morning.<div><br></div><div>Thans a lots.<br><br><div class="gmail_quote">On 21 June 2011 20:47, <span dir="ltr"><<a href="mailto:nek5000-users@lists.mcs.anl.gov">nek5000-users@lists.mcs.anl.gov</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">Hi Jean-Christophe,<br>
<br>
<br>
It looks like the only place in the code where ifadj shows up is in perturb.f:107<br>
<br>
if (ifnav.and.(.not.ifchar).and.( ifadj))call advabp_adjoint<br>
<br>
so you can set ifadj=.true. somewhere in userchk or usrdat -- just make that you also add<br>
<br>
include 'ADJOINT'<br>
<br>
there which is not part of 'TOTAL'<br>
<br>
Best,<br>
Aleks<div><div></div><div class="h5"><br>
<br>
<br>
<br>
On Tue, 21 Jun 2011, <a href="mailto:nek5000-users@lists.mcs.anl.gov" target="_blank">nek5000-users@lists.mcs.anl.<u></u>gov</a> wrote:<br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Hi Nek's,<br>
<br>
I have seen the adjoint perturbation mode has been added to the current.<br>
Where/how do I have to specify ifadj = .true. if I want to use it? Is there<br>
some parameters in the .rea file like for direct perturbation mode?<br>
<br>
Sincerely,<br>
<br>
On 9 June 2011 11:50, Jean-Christophe Loiseau <<a href="mailto:loiseau.jc@gmail.com" target="_blank">loiseau.jc@gmail.com</a>> wrote:<br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Hi Neks,<br>
<br>
I was just wondering if the adjoint linear solver has already been released<br>
in the current repo or not?<br>
<br>
Regards,<br>
<br>
<br>
On 6 May 2011 15:49, <<a href="mailto:nek5000-users@lists.mcs.anl.gov" target="_blank">nek5000-users@lists.mcs.anl.<u></u>gov</a>> wrote:<br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<br>
<br>
Dear Jean-Christophe,<br>
<br>
My apologies for the delay in getting back to you - we've been<br>
saturated with proposal writing on our end.<br>
<br>
Your approach looks reasonable -- though I would avoid at all<br>
costs the f90 type declarations and stick with the std f77-based<br>
approach to be consistent with the code.<br>
<br>
I'll look into your proposed approach. I should point out also<br>
that we've recently had an adjoint solver developed by some colleagues<br>
that is ready to go into the svn repo. I've just<br>
been delayed in getting it in because of travel and proposals.<br>
My hope is that this will be in the repo in about 10-15 days.<br>
<br>
Regarding the proper choice of gradm1 vs wgradm1, the question<br>
ultimately comes down to the degree of differentiability of the<br>
input --- functions in nek are only C0 and hence one-time differentiable.<br>
This applies also to the test functions. Thus, the std. practice<br>
is to move a derivative from the solution to the test functions (i.e.,<br>
to use the weak derivative) whenever the solution is to be differentiated<br>
twice.<br>
<br>
Regards,<br>
<br>
Paul<br>
<br>
<br>
<br>
<br>
<br>
On Thu, 5 May 2011, <a href="mailto:nek5000-users@lists.mcs.anl.gov" target="_blank">nek5000-users@lists.mcs.anl.<u></u>gov</a> wrote:<br>
<br>
Hi Nek's,<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<br>
<br>
I am interested into the linear adjoint perturbation. The equations read:<br>
<br>
-du/dt = (U.Grad) u - transpose(Grad(U)) u + 1/Re Laplacian(u) - grad(p)<br>
<br>
div(u) = 0<br>
<br>
Assuming t' = T-t, one can rewrite: du/dt' = (U.Grad) u -<br>
transpose(Grad(U))<br>
u + 1/Re Laplacian(u) - grad(p)<br>
Only small modifications to the perturbation mode are necessary, mainly<br>
to<br>
evaluate the tranpose gradient terms. Here is what I've done to evaluate<br>
these:<br>
<br>
subroutine conv_adj(conv_x,conv_y,conv_z,<u></u>u_x,u_y,u_z) ! used to be<br>
conv1n<br>
<br>
<br>
include 'SIZE'<br>
<br>
include 'DXYZ'<br>
<br>
include 'INPUT'<br>
<br>
include 'GEOM'<br>
<br>
include 'SOLN'<br>
<br>
include 'TSTEP'<br>
<br>
<br>
parameter (lt=lx1*ly1*lz1*lelt)<br>
<br>
real u_x (size(vx,1),size(vx,2),size(<u></u>vx,3),size(vx,4))<br>
<br>
real u_y (size(vx,1),size(vx,2),size(<u></u>vx,3),size(vx,4))<br>
<br>
real u_z (size(vx,1),size(vx,2),size(<u></u>vx,3),size(vx,4))<br>
<br>
<br>
real u_xx (size(vx,1),size(vx,2),size(<u></u>vx,3),size(vx,4))<br>
<br>
real u_xy (size(vx,1),size(vx,2),size(<u></u>vx,3),size(vx,4))<br>
<br>
real u_xz (size(vx,1),size(vx,2),size(<u></u>vx,3),size(vx,4))<br>
<br>
<br>
real u_yx (size(vy,1),size(vy,2),size(<u></u>vy,3),size(vy,4))<br>
<br>
real u_yy (size(vy,1),size(vy,2),size(<u></u>vy,3),size(vy,4))<br>
<br>
real u_yz (size(vy,1),size(vy,2),size(<u></u>vy,3),size(vy,4))<br>
<br>
<br>
real u_zx (size(vz,1),size(vz,2),size(<u></u>vz,3),size(vz,4))<br>
<br>
real u_zy (size(vz,1),size(vz,2),size(<u></u>vz,3),size(vz,4))<br>
<br>
real u_zz (size(vz,1),size(vz,2),size(<u></u>vz,3),size(vz,4))<br>
<br>
<br>
real conv_x (size(vx,1),size(vx,2),size(<u></u>vx,3),size(vx,4))<br>
<br>
real conv_y (size(vy,1),size(vy,2),size(<u></u>vy,3),size(vy,4))<br>
<br>
real conv_z (size(vz,1),size(vz,2),size(<u></u>vz,3),size(vz,4))<br>
<br>
<br>
*call gradm1(u_xx,u_xy,u_xz,u_x)*<br>
<br>
* call gradm1(u_yx,u_yy,u_yz,u_y)*<br>
<br>
* call gradm1(u_zx,u_zy,u_zz,u_z)*<br>
<br>
*<br>
*<br>
<br>
* conv_x = vx*u_xx + vy*u_yx + vz*u_zx*<br>
<br>
* conv_y = vx*u_xy + vy*u_yy + vz*u_zy*<br>
<br>
* conv_z = vx*u_xz + vy*u_yz + vz*u_zz*<br>
<br>
<br>
return<br>
<br>
end<br>
<br>
However, I'm not sure I've done it properly or not (arrays declaration<br>
for<br>
instance). Especially, I had no idea wheter I have to use the gradm1<br>
subroutine or its weak counterpart wgradm1. Finally in perturb.f, I've<br>
slightly change the advap subroutine the following way:<br>
<br>
subroutine advabp<br>
<br>
C<br>
<br>
C Eulerian scheme, add convection term to forcing function<br>
<br>
C at current time step.<br>
<br>
C<br>
<br>
include 'SIZE'<br>
<br>
include 'INPUT'<br>
<br>
include 'SOLN'<br>
<br>
include 'MASS'<br>
<br>
include 'TSTEP'<br>
<br>
C<br>
<br>
COMMON /SCRNS/ TA1 (LX1*LY1*LZ1*LELV)<br>
<br>
$ , TA2 (LX1*LY1*LZ1*LELV)<br>
<br>
$ , TA3 (LX1*LY1*LZ1*LELV)<br>
<br>
$ , TB1 (LX1*LY1*LZ1*LELV)<br>
<br>
$ , TB2 (LX1*LY1*LZ1*LELV)<br>
<br>
$ , TB3 (LX1*LY1*LZ1*LELV)<br>
<br>
C<br>
<br>
ntot1 = nx1*ny1*nz1*nelv<br>
<br>
ntot2 = nx2*ny2*nz2*nelv<br>
<br>
c<br>
<br>
if (if3d) then<br>
<br>
call opcopy (tb1,tb2,tb3,vx,vy,vz) ! Save<br>
velocity<br>
<br>
call opcopy (vx,vy,vz,vxp(1,jp),vyp(1,jp),<u></u>vzp(1,jp)) ! U <-- du<br>
<br>
* call conv_adj(ta1,ta2,ta3,tb1,tb2,<u></u>tb3)*<br>
<br>
* call opcopy (vx,vy,vz,tb1,tb2,tb3) ! Restore velocity*<br>
<br>
*c*<br>
<br>
* do i=1,ntot1*<br>
<br>
* tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,<u></u>ifield)*<br>
<br>
* bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)*<br>
<br>
* bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)*<br>
<br>
* bfzp(i,jp) = bfzp(i,jp)-tmp*ta3(i)*<br>
<br>
* enddo*<br>
<br>
c<br>
<br>
call convop (ta1,vxp(1,jp)) ! U.grad dU<br>
<br>
call convop (ta2,vyp(1,jp))<br>
<br>
call convop (ta3,vzp(1,jp))<br>
<br>
c<br>
<br>
do i=1,ntot1<br>
<br>
tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,<u></u>ifield)<br>
<br>
*bfxp(i,jp) = bfxp(i,jp) + tmp*ta1(i)*<br>
<br>
* bfyp(i,jp) = bfyp(i,jp) + tmp*ta2(i)*<br>
<br>
* bfzp(i,jp) = bfzp(i,jp) + tmp*ta3(i)*<br>
<br>
enddo<br>
<br>
c<br>
<br>
else ! 2D<br>
<br>
c<br>
<br>
call opcopy (tb1,tb2,tb3,vx,vy,vz) ! Save<br>
velocity<br>
<br>
call opcopy (vx,vy,vz,vxp(1,jp),vyp(1,jp),<u></u>vzp(1,jp)) ! U <-- dU<br>
<br>
call convop (ta1,tb1) ! du.grad U<br>
<br>
call convop (ta2,tb2)<br>
<br>
call opcopy (vx,vy,vz,tb1,tb2,tb3) ! Restore velocity<br>
<br>
c<br>
<br>
do i=1,ntot1<br>
<br>
tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,<u></u>ifield)<br>
<br>
bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)<br>
<br>
bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)<br>
<br>
enddo<br>
<br>
c<br>
<br>
call convop (ta1,vxp(1,jp)) ! U.grad dU<br>
<br>
call convop (ta2,vyp(1,jp))<br>
<br>
c<br>
<br>
do i=1,ntot1<br>
<br>
tmp = bm1(i,1,1,1)*vtrans(i,1,1,1,<u></u>ifield)<br>
<br>
bfxp(i,jp) = bfxp(i,jp)-tmp*ta1(i)<br>
<br>
bfyp(i,jp) = bfyp(i,jp)-tmp*ta2(i)<br>
<br>
enddo<br>
<br>
c<br>
<br>
endif<br>
<br>
c<br>
<br>
return<br>
<br>
end<br>
<br>
<br>
Any help to spot potential erros or to speed up these operation would be<br>
highly appreciated.<br>
<br>
Best regards,<br>
--<br>
Jean-Christophe<br>
<br>
______________________________<u></u>_________________<br>
</blockquote>
Nek5000-users mailing list<br>
<a href="mailto:Nek5000-users@lists.mcs.anl.gov" target="_blank">Nek5000-users@lists.mcs.anl.<u></u>gov</a><br>
<a href="https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users" target="_blank">https://lists.mcs.anl.gov/<u></u>mailman/listinfo/nek5000-users</a><br>
<br>
</blockquote>
<br>
<br>
<br>
--<br>
Jean-Christophe<br>
<br>
</blockquote>
<br>
<br>
<br>
-- <br>
Jean-Christophe<br>
<br>
</blockquote>
______________________________<u></u>_________________<br>
Nek5000-users mailing list<br>
<a href="mailto:Nek5000-users@lists.mcs.anl.gov" target="_blank">Nek5000-users@lists.mcs.anl.<u></u>gov</a><br>
<a href="https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users" target="_blank">https://lists.mcs.anl.gov/<u></u>mailman/listinfo/nek5000-users</a><br>
</div></div></blockquote></div><br><br clear="all"><br>-- <br>Jean-Christophe<br>
</div>