How to set a axisymmetric boundary condition

[steven]

Hi, guys

I am currently conducting a simulation of axisymmetric boundary condition. I generated my mesh with gmsh, in which I set the axisymmetric boundary with ID 4. In the .usr file, I also set a mapping from ID to BC in subroutuine userdat2. If the ID is 4, it will be given the BC 'A  '.

The code is as follows:

-------------------------------------------------------------------------------------------
      subroutine usrdat2()  ! This routine to modify mesh coordinates

c      implicit none

      include 'SIZE'
      include 'TOTAL'

c     Define loop variables
      integer iel, ifc, id

c     Loop over all elements
      do iel=1,nelv
c        Loop over all faces of each element
         do ifc=1,2*ndim

            ! Get the Physical ID defined in Gmsh for this face
            ! 'boundaryID' is the array where Nek5000 stores the IDs from the .re2 file
            id = boundaryID(ifc,iel)

            if (id.eq.1) then
               ! If ID is 1, set as Velocity boundary (v) - e.g., Piston
               cbc(ifc,iel,1) = 'v  '

            elseif (id.eq.2) then
               ! If ID is 2, set as Wall boundary (W)
               cbc(ifc,iel,1) = 'W  '

            elseif (id.eq.3) then
               ! If ID is 3, set as Outflow boundary (O)
               cbc(ifc,iel,1) = 'O  '

            elseif (id.eq.4) then
               cbc(ifc,iel,1) = 'A  '

            endif

           enddo
      enddo

      return
      end

-----------------------------------------------------------------------------------

And I have also set the Axisymmetry = yes in the .par file. The mesh is generated in the zone where y >= 0.

However, the calculation failed all the time, giving the message below:

-------------------------------------------------------------
Reading sjet2d.ma2                                                                                                    
 done :: partioning   0.30059     sec

 reading mesh
   readp_re2_mesh:pack/cr/unpack : 0.50E-05 0.13E-03 0.13E-04
   readp_re2_mesh:byte_read_mpi  : 0.77E-01
 reading boundary faces                  228 for ifield   1
   readp_re2_bc:pack/cr/unpack : 0.17E-05 0.41E-04 0.28E-05
   readp_re2_bc:byte_read_mpi  : 0.40E-01
 done :: read .re2 file    0.12     sec

 nelgt/nelgv/lelt:        1795        1795         266
 lx1/lx2/lx3/lxd:            6           4           6           9

 setup mesh topology
   Right-handed check complete for        1795 elements. OK.
 Checking axisymmetry BCs:
IFAXIS Error: BC must be
IFAXIS Error: BC must be

IFAXIS Error: BC must be

IFAXIS Error: BC must be
IFAXIS Error: BC must be
IFAXIS Error: BC must be x
IFAXIS Error: BC must be y
IFAXIS Error: BC must be z
IFAXIS Error: BC must be {
IFAXIS Error: BC must be <
IFAXIS Error: BC must be =
IFAXIS Error: BC must be >
IFAXIS Error: BC must be ?
IFAXIS Error: BC must be @
IFAXIS Error: BC must be A
IFAXIS Error: BC must be B
IFAXIS Error: BC must be C

.................
-------------------------------------------------------------------------------------------

Does anyone know what is going on with it? I have been trapped here for one day.

Any help will be appreciated. Thank you guys.

[steven]

  Due to the presence of some errors in the NEK5000 code, the error code output I provided earlier contained garbled characters. This is the error code output after I completed the debugging.
------------------------------------------------------------
 Checking axisymmetry BCs:
IFAXIS Error: BC must be  'A  ' for Face 1, Elem      265.
IFAXIS Error: BC must be  'A  ' for Face 1, Elem      266.
IFAXIS Error: BC must be  'A  ' for Face 1, Elem      267.
IFAXIS Error: BC must be  'A  ' for Face 1, Elem      268.
IFAXIS Error: BC must be  'A  ' for Face 1, Elem      269.
IFAXIS Error: BC must be  'A  ' for Face 1, Elem      270.
IFAXIS Error: BC must be  'A  ' for Face 1, Elem      271.
IFAXIS Error: BC must be  'A  ' for Face 1, Elem      272.
IFAXIS Error: BC must be  'A  ' for Face 1, Elem      273.
IFAXIS Error: BC must be ' A  ' for Face 1, Elem      274.
IFAXIS Error: BC must be  'A  ' for Face 1, Elem      275.
IFAXIS Error: BC must be  'A  ' for Face 1, Elem      376.
IFAXIS Error: BC must be  'A  ' for Face 1, Elem      377.

......

--------------------------------------------------------------------------------------------------

  Guys, I thought I have set all the axisymmetric boundaries correctly. But there must be something wrong. Please help me. Thanks a lot!!

[Fischer, Paul]

For axisymmtric, only ifc=1 is allowed to have the A bc --- that is, the elements on the axis boundary must be oriented such that "face 1" is on the axis.

Hth,

Paul

---

From: nek...@googlegroups.com <nek...@googlegroups.com> on behalf of steven <sijunc...@gmail.com>
Sent: Sunday, January 18, 2026 3:33 AM
To: Nek5000 <nek...@googlegroups.com>
Subject: [nek5000] How to set a axisymmetric boundary condition

 

--
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/fc801146-08f3-449f-ab23-c7d9a347b6f8n%40googlegroups.com.

[Bibhas Kumar]

Hi,
Apart from the BC ID, making physical groups in Gmsh, you would also need to take care of face IDs of each element. All quad elements are assigned face IDs 1-4, Nek requires the face ID of the BC elements that touch the axis of symmetry to be 1. This can be controlled in Gmsh. Show me your mesh so I can see if I can help.
Best,
Bibhas

[steven]

Dear Paul,

Thanks a lot for your advice.

But I thouhght the "face 1" of my mesh was on the axis, as the picture shows.

Have I grasped your point correctly?

Thanks again.

Steven  

 

[steven]

Dear Bibhas,

Thank you very much for your advice.

This is how I defined the physical groups in my .geo file.

As you said, should I define the axisymmetric lines first so that the ID of it can be 1?

Thanks again.

Steven

 

[Bibhas Kumar]

Thanks for providing the image. FACE IDs are different from Boundary IDs. Physical Curve tags the BC but not the faces. It is inherently controlled by the Tranfinite Surface commands when properly defined with the ordered corners.
 

When you are transfiniting the surfaces which lie on the 'A  ' BC (the bottom three (2, 3, and 5)), make sure the order of the corners is such that the points connecting the line which lie on the 'A  ' BC (lines 5, 9, and 14) are the first two corners. (Corners are points).

hth,

Bibhas

[steven]

Thank you for your advice, Bibhas.

  I have followed your instructions to the letter, but the program is still throwing errors. What is more, the elements triggering the chk_aixs error are not the ones located at y=0 at all. I strongly suspect there is a minor flaw in the element numbering mapping at this point in the source code. I wonder if you have encountered the same issue.

Thank you again.
Steven  

[Bibhas Kumar]

Yes, I have encountered a similar problem before; this issue is a little trickier than it seems (please see Prof. Fischer's previous replies). To resolve this, I created a range of geometries to get to the core of the issue, starting with a box, as my advisor suggested, and we eventually resolved this issue. I would encourage you to do the same.

I am attaching a .geo and .msh file that is very similar to your domain, and it worked for IFAXIS T and BC 'A  '. You can give it a try and play with dimensions in the .geo to make it similar to your domain.

hth,
Bibhas