It indeed solved my problems.<br>Thanks a lot.<br><br><div class="gmail_quote">On 13 September 2011 16:56, <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>
A couple of comments.<br>
<br>
1). I would use declaration of arrays vxold,...,qxold,...,f_x,...<br>
similar to declaration of vx, vy,... in<br>
<br>
nek5_svn/trunk/nek/SOLN :<br>
<br>
<br>
COMMON /VPTSOL/<br>
...<br>
$ , VX (LX1,LY1,LZ1,LELV)<br>
...<br>
<br>
<br>
2). To set an array on velocity mesh to zero, please use<br>
<br>
ntot = nx1*ny1*nz1*nelv<br>
call rzero(vxold,ntot)<br>
<br>
instead of<br>
<br>
vxold = 0.0D+00<br>
<br>
<br>
3). You do need local element number for setting the forcing in userf, namely<br>
<br>
e = gllel(eg)<br>
ffx = f_x(ix,iy,iz,e)<br>
<br>
Otherwise it will work correctly only on one core.<br>
<br>
<br>
Let me know if it helps.<br>
Best,<br>
Aleks<div><div></div><div class="h5"><br>
<br>
<br>
<br>
<br>
On Tue, 13 Sep 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've been trying to implement the selective frequency damping (Akervik E.,<br>
Brandt L., Henningson D.S., Hoepffner J., Marxen O., Schlatter P. 2006. *Steady<br>
solutions of the Navier-Stokes equations by selective frequency<br>
damping.*Physics of fluids<br>
*18*) which enables to compute stationnary solution. Here is what I have<br>
written in my .usr file:<br>
<br>
In userchk:<br>
<br>
common /vxold/ vxold(size(vx,1),size(vx,2),<u></u>size(vx,3),size(vx,4))<br>
common /vyold/ vyold(size(vy,1),size(vy,2),<u></u>size(vy,3),size(vy,4))<br>
common /vzold/ vzold(size(vz,1),size(vz,2),<u></u>size(vz,3),size(vz,4))<br>
<br>
common /qxold/ qxold(size(vx,1),size(vx,2),<u></u>size(vx,3),size(vx,4))<br>
common /qyold/ qyold(size(vy,1),size(vy,2),<u></u>size(vy,3),size(vy,4))<br>
common /qzold/ qzold(size(vz,1),size(vz,2),<u></u>size(vz,3),size(vz,4))<br>
<br>
common /for_x/ f_x(size(vx,1),size(vx,2),<u></u>size(vx,3),size(vx,4))<br>
common /for_y/ f_y(size(vy,1),size(vy,2),<u></u>size(vy,3),size(vy,4))<br>
common /for_z/ f_z(size(vz,1),size(vz,2),<u></u>size(vz,3),size(vz,4))<br>
<br>
real chi, omega_c<br>
integer step_filtr<br>
<br>
<br>
chi = .25<br>
omega_c = .1<br>
step_filtr = 020000<br>
<br>
if(istep.EQ.step_filtr.AND.<u></u>nid.EQ.0) write(*,*) 'Filtrage'<br>
<br>
if(istep.EQ.0) then<br>
vxold = 0.0D+00<br>
vyold = 0.0D+00<br>
vzold = 0.0D+00<br>
endif<br>
<br>
if(istep.LT.step_filtr) then<br>
qxold = 0.0D+00<br>
qyold = 0.0D+00<br>
qzold = 0.0D+00<br>
<br>
f_x = 0.0D+00<br>
f_y = 0.0D+00<br>
f_z = 0.0D+00<br>
<br>
call opcopy(vxold,vyold,vzold,vx,<u></u>vy,vz)<br>
endif<br>
<br>
if(istep.GE.step_filtr) then<br>
qxold = qxold + dt*omega_c * (vxold - qxold)<br>
qyold = qyold + dt*omega_c * (vyold - qyold)<br>
qzold = qzold + dt*omega_c * (vzold - qzold)<br>
<br>
f_x = -chi * (vxold - qxold)<br>
f_y = -chi * (vyold - qyold)<br>
f_z = -chi * (vzold - qzold)<br>
<br>
call opcopy(vxxold,vyyold,vzzold,<u></u>vxold,vyold,vzold)<br>
call opcopy(vxold,vyold,vzold,vx,<u></u>vy,vz)<br>
endif<br>
<br>
and in userf:<br>
<br>
common /for_x/ f_x(size(vx,1),size(vx,2),<u></u>size(vx,3),size(vx,4))<br>
common /for_y/ f_y(size(vx,1),size(vx,2),<u></u>size(vx,3),size(vx,4))<br>
common /for_z/ f_z(size(vx,1),size(vx,2),<u></u>size(vx,3),size(vx,4))<br>
<br>
<br>
integer e,f,eg<br>
c e = gllel(eg)<br>
<br>
<br>
c Note: this is an acceleration term, NOT a force!<br>
c Thus, ffx will subsequently be multiplied by rho(x,t).<br>
<br>
<br>
ffx = f_x(ix,iy,iz,eg)<br>
ffy = f_y(ix,iy,iz,eg)<br>
ffz = f_z(ix,iy,iz,eg)<br>
<br>
<br>
I have tried to compute with such routines the base flow for a 2D lid-driven<br>
cavity flow at Re = 8500, however depending on the number of spectral<br>
elements I use, I sometimes get a stationnary solution that looks slightly<br>
as the one I'm looking for (but some differences near the lid) and sometimes<br>
I get one which is actually not stationnary. I applied the same algorithm<br>
for a finite-difference code and encountered no problem, so I was wondering<br>
then if I missed something specific to spectral elements or Nek 5000 (such<br>
as a mass matrix multplication for instance) or to F77 (which I'm not at<br>
ease with, see my definition of the common for instance)?<br>
<br>
Sincerely yours,<br>
<br>
-- <br>
Jean-Christophe<br>
<br>
</blockquote></div></div>
______________________________<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>
</blockquote></div><br><br clear="all"><br>-- <br>Jean-Christophe<br>