[Nek5000-users] Stream function

nek5000-users at lists.mcs.anl.gov nek5000-users at lists.mcs.anl.gov
Thu Jul 12 03:37:29 CDT 2018


Ok, so can you please just try the attached code in userchck? It should 
work (at least does for us) in 2D. Essentially, the steps are:
1. copmute z-vorticity in a new array vorticity
2. construct the rhs by multiplying with -bm1
3. setup the Helmholtz problem and tolerances
4. solve the Helmholtz problem and storing the results in a new array phi
5. write out everything in two files.

In this code you do not use the temperature at all, which I guess is 
fine for a 2D case (which I guess is what you look at).

best regards,
Philipp

On 2018-07-10 22:33, Philipp Schlatter wrote:
> Dear Hannah,
> sorry that I only see the post now. Are you still having problems with 
> the stream function, or does it work now?
> Best regards,
> Philipp
> 
> On 2018-07-02 16:40, nek5000-users at lists.mcs.anl.gov wrote:
>> Hi JC,
>>
>> Yes, the mv_wall example computes the streamfunction, but it applies 
>> tmask, tmult that I could not find where they are defined in the 
>> file.  I was experimenting with different masks.
>>
>> Do you know how postnek compute the streamfunction?
>>
>> I think it is computed in the subroutine strfct(psi) in postnek2.f, 
>> but I am not sure what happens in the routine?
>>
>> Thanks,
>>
>> Hannah
>> ________________________________________
>> From: Nek5000-users [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: Friday, June 29, 2018 3:47 AM
>> To: nek5000-users at lists.mcs.anl.gov
>> Subject: Re: [Nek5000-users] Stream function
>>
>> Hi Hannah,
>>
>> I tried to implement the streamfunction as well (for a 2D cylinder flow)
>> and it would appear that, doing naïvely, the tmask is set to zero
>> everywhere. Hence, when rhs * tmask is computed in the hmholtz
>> subroutine, rhs is set to zero everywhere and it then makes sense that
>> the returned streamfunction is zero everywhere as well.
>>
>> Not sure how to do fix this properly though. Any help would be
>> appreciated indeed :)
>>
>>
>> ++
>> JC
>>
>>
>> On 06/28/2018 09:41 PM, nek5000-users at lists.mcs.anl.gov wrote:
>>> Hi Philipp,
>>>
>>> Thanks very much for your comments.  Yes, the delta* works only for 
>>> the parallel plate.
>>>
>>> The streamfunction computation with the Blasius example is still not 
>>> working. I suspect I am not using the correct mask (vmult /tmult, 
>>> imsh as well).
>>>
>>> Using the Blasius example, v1mask, vmult (imsh = 1), I dumped the 
>>> streamfuntion to the temperature field. However, it only contains 
>>> negative values and does not look physical (I pasted userchk below).
>>>
>>> Another Nek Example, “mv_wall”, computes the streamfunction from the 
>>> Possison equation using vorticity. It uses tmask, tmult, and imsh = 2.
>>>
>>> I tried implementing what was done in the “mv_wall” example in the 
>>> Blasius example, but the computed streamfunction is zero everywhere.
>>>
>>> Any ideas what I might be doing wrong?
>>>
>>> Any advice would be appreciated!
>>>
>>> Thanks again,
>>>
>>> Hannah
>>>
>>> c-----------------------------------------------------------------------
>>>         subroutine userchk
>>>         include 'SIZE'
>>>         include 'TOTAL'
>>>
>>>         common /myblas/ ub(lx1,ly1,lz1,lelt),vb(lx1,ly1,lz1,lelt)
>>>
>>>         common /blasiusr/ u(lx1,ly1,lz1,lelv),v(lx1,ly1,lz1,lelv)
>>>         common /blasiusg/ bin(lx1,ly1,lz1,lelv)
>>>         common /scrns/    dum(lx1*ly1*lz1*lelv)
>>>
>>>         common /xtream/ psi(lx1*ly1*lz1*lelt)
>>>        $              , rhs(lx1*ly1*lz1*lelt)
>>>        $              , h1 (lx1*ly1*lz1*lelt)
>>>        $              , h2 (lx1*ly1*lz1*lelt)
>>>
>>>         parameter (lt=lx1*ly1*lz1*lelt)
>>>         common /scrns/ vort(lt,3),w1(lt),w2(lt)
>>>
>>>         integer bin
>>>         real normu(50),normv(50)
>>>
>>>         real xmax,xmin,dx
>>>         save xmax,xmin,dx
>>>
>>>         n = nx1*ny1*nz1*nelv
>>>         ifheat=.false.
>>>         call comp_vort3(vort , w1, w2, vx, vy, vz)
>>>         call copy (t,vort,n)   ! Vorticity --> T
>>>         ifto = .true.          ! Dump vorticity as T
>>> c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
>>> c     Compute streamfunction from vorticity
>>>         call col3 (rhs,t,bm1,n) ! B*omega
>>>         call rone (h1,n)
>>>         call rzero(h2,n)
>>>
>>>         ifield = 1
>>>         tol    = param(22)
>>>         imsh   = 1
>>>         call hmholtz('psi ',psi,rhs,h1,h2,v1mask,vmult,imsh,
>>>        $        tol,200,1)
>>>         call copy (t,psi,n)    ! Psi --> T (for output)
>>> c. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
>>>
>>>         nfldt = 1 + npscal
>>>
>>>         if (istep.eq.0) call outpost2(vx,vy,vz,pr,t,nfldt,'   ')
>>>
>>>         do i=1,n
>>>            dum(i) = 1.-vx(i,1,1,1)
>>>         enddo
>>>         tot_deficit = glsc2(dum,bm1,n)
>>>         xmin        = glmin(xm1,n)
>>>         xmax        = glmax(xm1,n)
>>>         delta_star  = tot_deficit/(xmax-xmin)
>>>         tmin        = glmin(t,n)
>>>         tmax        = glmax(t,n)
>>>         vortmin      = glmin(vort,n)
>>>         vortmax      = glmax(vort,n)
>>>         if (nid.eq.0) write(6,1) istep,time,delta_star,vortmin,vortmax
>>>         if (nid.eq.0) write(6,1) istep,tmin,tmax
>>>       1 format(i8,1p4e13.5,'  blasius: delta*')
>>>
>>>         return
>>>         end
>>> c-----------------------------------------------------------------------
>>>
>>>
>>> ________________________________________
>>> From: Nek5000-users [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: Tuesday, June 26, 2018 1:22 PM
>>> To: nek5000-users at lists.mcs.anl.gov
>>> Subject: Re: [Nek5000-users] Stream function
>>>
>>> Hi,
>>> To correctly compute the stream function, the rhs of the Poisson
>>> equation needs to be the (negative) vorticity. Therefore, the line
>>>
>>>       call col2(rhs,bm1,n)
>>>
>>> is wrong, and should be something like
>>>
>>>       call col3(rhs, vorticity, bm1, n)
>>>
>>> where vorticity should be the computed vorticity using comp_vort.
>>>
>>> I guess your calculation of delta* only works for parallel boundary 
>>> layers.
>>>
>>> Philipp
>>>
>>>
>>>
>>> On 2018-06-26 18:35, nek5000-users at lists.mcs.anl.gov wrote:
>>>> Hi Phillip,
>>>>
>>>> I am using the “Blasius” example test outputting the streamfunction 
>>>> as a
>>>> scalar.
>>>>
>>>> I included my code below. I am not sure what I am doing wrong.
>>>>
>>>> In userchk, I added “compute streamfunction from vorticity” from the
>>>> psi_omega example, but it does not compute the psi properly… (max_psi =
>>>> min_psi = 0 in every step)
>>>>
>>>> c----------------------------------------------------------------------- 
>>>>
>>>>
>>>>          subroutine userchk
>>>>
>>>>          include 'SIZE'
>>>>
>>>>          include 'TOTAL'
>>>>
>>>>          common /myblas/ ub(lx1,ly1,lz1,lelt),vb(lx1,ly1,lz1,lelt)
>>>>
>>>>          common /blasiusr/ u(lx1,ly1,lz1,lelv),v(lx1,ly1,lz1,lelv)
>>>>
>>>>          common /blasiusg/ bin(lx1,ly1,lz1,lelv)
>>>>
>>>>          common /scrns/    dum(lx1*ly1*lz1*lelv)
>>>>
>>>>          common /xtream/ psi(lx1*ly1*lz1*lelt)
>>>>
>>>>         $              , rhs(lx1*ly1*lz1*lelt)
>>>>
>>>>         $              , h1 (lx1*ly1*lz1*lelt)
>>>>
>>>>         $              , h2 (lx1*ly1*lz1*lelt)
>>>>
>>>>          integer bin
>>>>
>>>>          real normu(50),normv(50)
>>>>
>>>>          real xmax,xmin,dx
>>>>
>>>>          save xmax,xmin,dx
>>>>
>>>>          n = nx1*ny1*nz1*nelv
>>>>
>>>> c. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
>>>>
>>>> c     Compute streamfunction from vorticity
>>>>
>>>>          call col2(rhs,bm1,n)
>>>>
>>>>          call rone (h1 ,n)
>>>>
>>>>          call rzero(h2 ,n)
>>>>
>>>>          tol = param(22)
>>>>
>>>>          imsh   = 1
>>>>
>>>>          call hmholtz('psi ',psi,rhs,h1,h2,v1mask,vmult,imsh,tol,200,1)
>>>>
>>>> c. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
>>>>
>>>>          if (istep.eq.0) call outpost(vx,vy,vz,pr,psi,'   ')
>>>>
>>>> c     Compute delta*
>>>>
>>>>          do i=1,n
>>>>
>>>>             dum(i) = 1.-vx(i,1,1,1)
>>>>
>>>>          enddo
>>>>
>>>>          tot_deficit = glsc2(dum,bm1,n)
>>>>
>>>>          xmin        = glmin(xm1,n)
>>>>
>>>>          xmax        = glmax(xm1,n)
>>>>
>>>>          delta_star  = tot_deficit/(xmax-xmin)
>>>>
>>>>          psimin      = glmin(psi,n)
>>>>
>>>>          psimax      = glmax(psi,n)
>>>>
>>>>          if (nid.eq.0) write(6,1) istep,time,delta_star,psimin,psimax
>>>>
>>>>        1 format(i8,1p4e13.5,'  blasius: delta*')
>>>>
>>>>          return
>>>>
>>>>          end
>>>>
>>>>
>>>>
>>>> _______________________________________________
>>>> Nek5000-users mailing list
>>>> Nek5000-users at lists.mcs.anl.gov
>>>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>>>>
>>> _______________________________________________
>>> Nek5000-users mailing list
>>> Nek5000-users at lists.mcs.anl.gov
>>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>>> _______________________________________________
>>> Nek5000-users mailing list
>>> Nek5000-users at lists.mcs.anl.gov
>>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>>
>> _______________________________________________
>> Nek5000-users mailing list
>> Nek5000-users at lists.mcs.anl.gov
>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>> _______________________________________________
>> Nek5000-users mailing list
>> Nek5000-users at lists.mcs.anl.gov
>> https://lists.mcs.anl.gov/mailman/listinfo/nek5000-users
>>
-------------- next part --------------
c-----------------------------------------------------------------------
C
C  USER SPECIFIED ROUTINES:
C
C     - boundary conditions
C     - initial conditions
C     - variable properties
C     - local acceleration for fluid (a)
C     - forcing function for passive scalar (q)
C     - general purpose routine for checking errors etc.
C
c-----------------------------------------------------------------------
      subroutine uservp (ix,iy,iz,eg)

      return
      end
c-----------------------------------------------------------------------
      subroutine userf  (ix,iy,iz,ieg)

      return
      end
c-----------------------------------------------------------------------
      subroutine userq  (ix,iy,iz,eg)


      return
      end
c-----------------------------------------------------------------------
      subroutine userchk

      include 'SIZE'
      include 'INPUT'           ! param
      include 'SOLN'            ! vx, vy, vz, pr, t, v2mask, vmult
      include 'MASS'            ! bm1
      

      real dummy1(lx1,ly1,lz1,lelv)
      real dummy2(lx1,ly1,lz1,lelv)
      real vorticity(lx1,ly1,lz1,lelv)
      real phi(lx1,ly1,lz1,lelv)
      real h1(lx1,ly1,lz1,lelv)
      real h2(lx1,ly1,lz1,lelv)
      real rhs(lx1,ly1,lz1,lelv)

      integer n, maxit
      real toler, p22
      
      n = nelv*lx1*ly1*lz1
c
c     Compute streamfunction phi for a 2D velocity field with
c     all Dirichlet boundary conditions (i.e. phi=0 along all boundaries)
c
c     Solve the Poisson equation \/^2 phi = -w
c     with w the vorticity.
c
c     initialize the vectors to zero
      call rzero(vorticity,n)
c     compute vorticity
      call comp_vort3(vorticity,dummy1,dummy2,vx,vy,vz)                                                    
      call chsign(vorticity,n)
c     multiply with the mass matrix
      call col3(rhs,vorticity,bm1,n)
c     setup Poisson problem: h2=0, h1=1
      call rzero(h2,n)
      call rone(h1,n)
c     set solver tolerance
      toler = 1e-14
      p22 = param(22)
      param(22) = toler
      maxit = 10000      
c     Solve Poisson equation:    
      call hmholtz('PHI ',phi,rhs,h1,h2,v2mask,vmult,1,
     $     toler,maxit,1)
c     reset tolerance
      param(22) = p22

c     output streamfunction phi
      call outpost(rhs,phi,vz,pr,t,'stf')
c     output normal velocity field
      call outpost(vx,vy,vz,pr,t,'   ')
      
      call exitt()

      return
      end
c-----------------------------------------------------------------------
      subroutine userbc (ix,iy,iz,iside,ieg)

      return
      end
c-----------------------------------------------------------------------
      subroutine useric (ix,iy,iz,ieg)

      return
      end
c-----------------------------------------------------------------------
      subroutine usrdat

      return
      end
c-----------------------------------------------------------------------
      subroutine usrdat2
c
      return
      end
c-----------------------------------------------------------------------
      subroutine usrdat3
c
      return
      end
c-----------------------------------------------------------------------


More information about the Nek5000-users mailing list