[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