[Nek5000-users] Selective frequency damping
nek5000-users at lists.mcs.anl.gov
nek5000-users at lists.mcs.anl.gov
Tue Sep 13 09:35:51 CDT 2011
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20110913/42d20bf2/attachment.html>
More information about the Nek5000-users
mailing list