Hi Nek's,<br><br>I've been trying to implement the selective frequency damping (Akervik E., Brandt L., Henningson D.S., Hoepffner J., Marxen O., Schlatter P. 2006. <i>Steady solutions of the Navier-Stokes equations by selective frequency damping.</i> Physics of fluids <b>18</b>) which enables to compute stationnary solution. Here is what I have written in my .usr file:<br>
<br><div style="margin-left: 40px;">In userchk:<br><br> common /vxold/ vxold(size(vx,1),size(vx,2),size(vx,3),size(vx,4))<br> common /vyold/ vyold(size(vy,1),size(vy,2),size(vy,3),size(vy,4))<br> common /vzold/ vzold(size(vz,1),size(vz,2),size(vz,3),size(vz,4))<br>
<br> common /qxold/ qxold(size(vx,1),size(vx,2),size(vx,3),size(vx,4))<br> common /qyold/ qyold(size(vy,1),size(vy,2),size(vy,3),size(vy,4))<br> common /qzold/ qzold(size(vz,1),size(vz,2),size(vz,3),size(vz,4))<br>
<br> common /for_x/ f_x(size(vx,1),size(vx,2),size(vx,3),size(vx,4))<br> common /for_y/ f_y(size(vy,1),size(vy,2),size(vy,3),size(vy,4))<br> common /for_z/ f_z(size(vz,1),size(vz,2),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.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,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,vxold,vyold,vzold)<br> call opcopy(vxold,vyold,vzold,vx,vy,vz)<br> endif<br><br>and in userf:<br><br> common /for_x/ f_x(size(vx,1),size(vx,2),size(vx,3),size(vx,4))<br>
common /for_y/ f_y(size(vx,1),size(vx,2),size(vx,3),size(vx,4))<br> common /for_z/ f_z(size(vx,1),size(vx,2),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></div><br>I have tried to compute with such routines the base flow for a 2D lid-driven cavity flow at Re = 8500, however depending on the number of spectral elements I use, I sometimes get a stationnary solution that looks slightly as the one I'm looking for (but some differences near the lid) and sometimes I get one which is actually not stationnary. I applied the same algorithm for a finite-difference code and encountered no problem, so I was wondering then if I missed something specific to spectral elements or Nek 5000 (such as a mass matrix multplication for instance) or to F77 (which I'm not at ease with, see my definition of the common for instance)?<br>
<br>Sincerely yours,<br><br>-- <br>Jean-Christophe<br>