Hi Neks,<div><br></div><div>I was just wondering if the adjoint linear solver has already been released in the current repo or not?</div><div><br></div><div>Regards,<br><br><div class="gmail_quote">On 6 May 2011 15:49, <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;"><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 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<div><div></div><div class="h5"><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.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>
<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 - transpose(Grad(U))<br>
u + 1/Re Laplacian(u) - grad(p)<br>
Only small modifications to the perturbation mode are necessary, mainly 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_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(vx,3),size(vx,4))<br>
<br>
real u_y (size(vx,1),size(vx,2),size(vx,3),size(vx,4))<br>
<br>
real u_z (size(vx,1),size(vx,2),size(vx,3),size(vx,4))<br>
<br>
<br>
real u_xx (size(vx,1),size(vx,2),size(vx,3),size(vx,4))<br>
<br>
real u_xy (size(vx,1),size(vx,2),size(vx,3),size(vx,4))<br>
<br>
real u_xz (size(vx,1),size(vx,2),size(vx,3),size(vx,4))<br>
<br>
<br>
real u_yx (size(vy,1),size(vy,2),size(vy,3),size(vy,4))<br>
<br>
real u_yy (size(vy,1),size(vy,2),size(vy,3),size(vy,4))<br>
<br>
real u_yz (size(vy,1),size(vy,2),size(vy,3),size(vy,4))<br>
<br>
<br>
real u_zx (size(vz,1),size(vz,2),size(vz,3),size(vz,4))<br>
<br>
real u_zy (size(vz,1),size(vz,2),size(vz,3),size(vz,4))<br>
<br>
real u_zz (size(vz,1),size(vz,2),size(vz,3),size(vz,4))<br>
<br>
<br>
real conv_x (size(vx,1),size(vx,2),size(vx,3),size(vx,4))<br>
<br>
real conv_y (size(vy,1),size(vy,2),size(vy,3),size(vy,4))<br>
<br>
real conv_z (size(vz,1),size(vz,2),size(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 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),vzp(1,jp)) ! U <-- du<br>
<br>
* call conv_adj(ta1,ta2,ta3,tb1,tb2,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,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,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),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,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,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>
</blockquote></div></div>
_______________________________________________<br>
Nek5000-users mailing list<br>
<a href="mailto:Nek5000-users@lists.mcs.anl.gov" target="_blank">Nek5000-users@lists.mcs.anl.gov</a><br>
<a href="https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users" target="_blank">https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users</a><br>
</blockquote></div><br><br clear="all"><br>-- <br>Jean-Christophe<br>
</div>