[Nek5000-users] dissipation

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

Hi Ilias,


How about:


      subroutine dissipation(sij2)


      include 'SIZE'

      include 'TOTAL'


      real sij2(lx1*ly1*lz1,lelt)


      common /scrns/         sij (lx1*ly1*lz1,6,lelv)


      parameter (lr=lx1*ly1*lz1)

      common /scruz/         ur(lr),us(lr),ut(lr)

     $                     , vr(lr),vs(lr),vt(lr)

     $                     , wr(lr),ws(lr),wt(lr)


      integer e


      nij = 3+3*(ndim-2)


      call comp_sij(sij,nij,vx,vy,vz,ur,us,ut,vr,vs,vt,wr,ws,wt)


      nxyz = nx1*ny1*nz1

      n    = nx1*ny1*nz1*nelv


      call rzero(sij2,n)


      do e=1,nelv

      do j=1,nij

      do i=1,nxyz

         sij2(i,e)=sij2(i,e)+sij(i,j,e)**2

      enddo

      enddo

      enddo


      return

      end

?

This should function correctly.

Best, Paul


________________________________
From: nek5000-users-bounces at lists.mcs.anl.gov [nek5000-users-bounces at lists.mcs.anl.gov] on behalf of nek5000-users at lists.mcs.anl.gov [nek5000-users at lists.mcs.anl.gov]
Sent: Sunday, November 08, 2015 4:37 PM
To: nek5000-users at lists.mcs.anl.gov
Subject: [Nek5000-users] dissipation

Dear neks,
I am trying to user comp_sij to compute dissipation in every point of the domain.
The data storage is organised according to
      real sij(lx1*ly1*lz1,nij,lelv)
so the index of a component of tensor
precedes number of element, and
the is no continuous block of lx1*ly1*lz1,lelv points for an element of tensor.
And so, the global functions like col3 of glsc can not be applied.
At least I don't understand how.

To bypass this, I copied comp_sij to my usr file and changed the order of
last indices in sij
      real sij(lx1*ly1*lz1,lelv,nij)


After that I use the followng code to compute viscous dissipation

       nij = 3
       if (if3d) nij=6

      call il_comp_sij(sij,nij,vx,vy,vz,ur,us,ut,vr,vs,vt,wr,ws,wt)

      call rzero(il_dissip,lv)

      do j=1,nij
      call col3(sij2,sij(1,1,1,1,j),sij(1,1,1,1,j),nxyz)
      il_dissip=il_dissip+sij2
       if (if3d .and. j.ge.4) then
      il_dissip=il_dissip+sij2
       end if
      if ( .not. if3d .and. j.eq.3 ) then
      il_dissip=il_dissip+sij2
      end if
      il_dissip=2*param(2)*il_dissip
      enddo

Is is a fair way to compute dissipation?

Thanks, Ilias


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/nek5000-users/attachments/20151109/a3aa7052/attachment.html>