cyclic boundary conditions - 2D or 3D?

49 views
Skip to first unread message

Philipp Schlatter

unread,
Mar 25, 2025, 6:51:03 AM3/25/25
to nek5000, Mammadli, Isa

Hi,

I had a quick question. Do the cyclic boundary conditions, as implemented right now in Nek5000, also work for a completely arbitrary rotation in 3D, or only in 2D (i.e. around an axis)? We got confused because the main routine, rotate_cyc in navier1.f:

02648       nface = 2*ndim
02649       do e=1,nelfld(ifield)
02650       do f=1,nface
02651 
02652          if(cbc(f,e,ifield) .eq. 'P  '.or.cbc(f,e,ifield).eq.'p  ')then
02653 
02654             call facind2 (js1,jf1,jskip1,js2,jf2,jskip2,f)
02655             if (idir.eq.1) then
02656               k=0
02657               do j2=js2,jf2,jskip2
02658               do j1=js1,jf1,jskip1
02659                k=k+1
02660 
02661                dotprod = unx(k,1,f,e)*ym1(j1,j2,1,e)
02662      $                  -uny(k,1,f,e)*xm1(j1,j2,1,e)
02663                ifxy = .false.
02664                if (abs(unz(k,1,f,e)).lt.0.0001) ifxy = .true.
02665 
02666                cost =  unx(k,1,f,e)
02667                sint =  uny(k,1,f,e)
02668                rnor = ( r1(j1,j2,1,e)*cost + r2(j1,j2,1,e)*sint )
02669                rtn1 = (-r1(j1,j2,1,e)*sint + r2(j1,j2,1,e)*cost )
02670 
02671                if (ifxy.and.dotprod .ge. 0.0) then 
02672                   r1(j1,j2,1,e) = rnor
02673                   r2(j1,j2,1,e) = rtn1
02674                elseif (ifxy) then
02675                   r1(j1,j2,1,e) =-rnor
02676                   r2(j1,j2,1,e) =-rtn1
02677                endif
02678               enddo
02679               enddo

seems to only do a 2D rotation. Or am I missing something? Essentially, we would like to map two completely arbitrary planes with each other.

Thanks,
philipp


Fischer, Paul

unread,
Mar 25, 2025, 8:00:54 AM3/25/25
to Philipp Schlatter, nek5000, Mammadli, Isa
Hi Philipp,

I think you could modify the setup to handle the arbitrary case, given that the math on every timestep would be the same.  It's just that you'd need to be careful in the set up.  One reason not to go for the general case at the time is that you then also need to worry about users running with element faces that are not topologically congruent, e.g., one face "rotated" w.r.t.  the other.   The most common use case seemed to be people who wanted to run just a piece of pie in an axisymmetric geometry.  Even that case is fraught with challenges because of the pole singularity.

hth,
Paul


From: nek...@googlegroups.com <nek...@googlegroups.com> on behalf of Philipp Schlatter <psch...@gmail.com>
Sent: Tuesday, March 25, 2025 5:50 AM
To: nek5000 <nek...@googlegroups.com>
Cc: Mammadli, Isa <isa.ma...@fau.de>
Subject: [nek5000] cyclic boundary conditions - 2D or 3D?
 
--
You received this message because you are subscribed to the Google Groups "Nek5000" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nek5000+u...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/nek5000/e4b4dcd5-8d83-46c2-b85b-2c4536cf9ec2%40gmail.com.

Philipp Schlatter

unread,
Mar 25, 2025, 9:50:29 AM3/25/25
to Fischer, Paul, nek5000, Mammadli, Isa

Thanks Paul. Actually, the latest version of Nek has a slightly updated rotate_cyc routine, which always does the 2D rotation, but permits a potential error in z (I think). Let me think a bit more about our helical pipes then :-)


Best regards,
Philipp

Reply all
Reply to author
Forward
0 new messages