[Nek5000-users] Selective frequency damping

[nek5000-users at lists.mcs.anl.gov]

Hi Jean-Christophe,

A couple of comments.

1). I would use declaration of arrays vxold,...,qxold,...,f_x,...
similar to declaration of vx, vy,... in

nek5_svn/trunk/nek/SOLN :


    COMMON /VPTSOL/
...
      $               , VX     (LX1,LY1,LZ1,LELV)
...


2). To set an array on velocity mesh to zero, please use

       ntot = nx1*ny1*nz1*nelv
       call rzero(vxold,ntot)

instead of

       vxold = 0.0D+00


3). You do need local element number for setting the forcing in userf, 
namely

       e   = gllel(eg)
       ffx = f_x(ix,iy,iz,e)

Otherwise it will work correctly only on one core.


Let me know if it helps.
Best,
Aleks




On Tue, 13 Sep 2011, nek5000-users at lists.mcs.anl.gov wrote:

> Hi Nek's,
>
> I've been trying to implement the selective frequency damping (Akervik E.,
> Brandt L., Henningson D.S., Hoepffner J., Marxen O., Schlatter P. 2006. *Steady
> solutions of the Navier-Stokes equations by selective frequency
> damping.*Physics of fluids
> *18*) which enables to compute stationnary solution. Here is what I have
> written in my .usr file:
>
> In userchk:
>
>       common /vxold/ vxold(size(vx,1),size(vx,2),size(vx,3),size(vx,4))
>       common /vyold/ vyold(size(vy,1),size(vy,2),size(vy,3),size(vy,4))
>       common /vzold/ vzold(size(vz,1),size(vz,2),size(vz,3),size(vz,4))
>
>       common /qxold/ qxold(size(vx,1),size(vx,2),size(vx,3),size(vx,4))
>       common /qyold/ qyold(size(vy,1),size(vy,2),size(vy,3),size(vy,4))
>       common /qzold/ qzold(size(vz,1),size(vz,2),size(vz,3),size(vz,4))
>
>       common /for_x/ f_x(size(vx,1),size(vx,2),size(vx,3),size(vx,4))
>       common /for_y/ f_y(size(vy,1),size(vy,2),size(vy,3),size(vy,4))
>       common /for_z/ f_z(size(vz,1),size(vz,2),size(vz,3),size(vz,4))
>
>       real chi, omega_c
>       integer step_filtr
>
>
>       chi = .25
>       omega_c = .1
>       step_filtr = 020000
>
>       if(istep.EQ.step_filtr.AND.nid.EQ.0) write(*,*) 'Filtrage'
>
>       if(istep.EQ.0) then
>       vxold = 0.0D+00
>       vyold = 0.0D+00
>       vzold = 0.0D+00
>       endif
>
>       if(istep.LT.step_filtr) then
>       qxold = 0.0D+00
>       qyold = 0.0D+00
>       qzold = 0.0D+00
>
>       f_x = 0.0D+00
>       f_y = 0.0D+00
>       f_z = 0.0D+00
>
>       call opcopy(vxold,vyold,vzold,vx,vy,vz)
>       endif
>
>       if(istep.GE.step_filtr) then
>       qxold = qxold + dt*omega_c * (vxold - qxold)
>       qyold = qyold + dt*omega_c * (vyold - qyold)
>       qzold = qzold + dt*omega_c * (vzold - qzold)
>
>       f_x = -chi * (vxold - qxold)
>       f_y = -chi * (vyold - qyold)
>       f_z = -chi * (vzold - qzold)
>
>       call opcopy(vxxold,vyyold,vzzold,vxold,vyold,vzold)
>       call opcopy(vxold,vyold,vzold,vx,vy,vz)
>       endif
>
> and in userf:
>
>       common /for_x/ f_x(size(vx,1),size(vx,2),size(vx,3),size(vx,4))
>       common /for_y/ f_y(size(vx,1),size(vx,2),size(vx,3),size(vx,4))
>       common /for_z/ f_z(size(vx,1),size(vx,2),size(vx,3),size(vx,4))
>
>
>      integer e,f,eg
> c     e = gllel(eg)
>
>
> c     Note: this is an acceleration term, NOT a force!
> c     Thus, ffx will subsequently be multiplied by rho(x,t).
>
>
>     ffx = f_x(ix,iy,iz,eg)
>     ffy = f_y(ix,iy,iz,eg)
>     ffz = f_z(ix,iy,iz,eg)
>
>
> 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)?
>
> Sincerely yours,
>
> -- 
> Jean-Christophe
>