Feature/cyclic bc

[IsaMammadli]

In spectral element method, Direct Stiffness Sumation (DSS) is used to enforces continuity across element boundary nodes. Neko (as well as Nek5000) achieves this with gather-scatter operation (gs%op(..,.., GS_OP_ADD). Since the GLL nodes at periodic faces are also neighbors they are subject to GS operation. When the periodic faces do not align with each other, this summation can no longer be performed in Cartesian basis (see figure below).

In Nek5000 this is handled with cyclic boundary conditions, which temporarily rotates the vectors at periodic faces to an intermediate basis to apply DSS correctly. This occurs as following: [ rotate forward ] → [ gather ] → [ scatter ] → [ rotate backward ]

For bent channels, toroidal and helical pipes we rotate the vectors into Frenet-Serret (a.k.a TNB) basis which is a moving reference frame for differentiable curves. In this proposed approach, I inserted rotate_cyc( ux, uy, uz, idir, coef) subroutine in relevant locations of the code, as in the following diagram:

• idir=1 is for Cartesian → TNB rotation and idir=0 is TNB → Cartesian.
• The 5 lines of rotate-gs calls are analogous to opdssum subroutine in Nek5000, but kept explicit intentionally since gs%op calls are not always consecutive in Neko (e.g. in fused_cpld_cg.F90 there are device_event_sync and blstz%apply in between gs%op calls.)
• A new subroutine coef_generate_cyclic_bc computes rotation angle as part of coef%init. It uses surface normals which are also generated in coef.f90
• Here cyc_msk array is also generated. It is a mask array that helps to iterate GLL nodes at periodic faces and components of rotation matrix.
• CPU and device implementation of rotate_cyc have two versions for rank 1 and rank 4 arrays.
• Cyclic BC is enabled in .case file manually.

There can be different ways to handle these steps, so I am open to suggestions and I list here some of them for your comments.

1. Do we want to combine rotate - gs process to make it more compact like opdssum ?
2. Currently rotate_cyc only handles 2D rotations, which are sufficient for toroids and helical pipes. This is the same as in Nek5000 because the rotation is geometry dependent and difficult to generalize. However, it would require minor changes to adjust it for 3D case, if any. (Unfortunately we dont have an example for that at this point, unless we orient pipes weirdly.)
3. I think keeping cyclic BC as a manual switch in the user file is a good idea, because the user must be aware of the requirement, its implementation details and if current implementation suffices without any modifications. I have a check_cyclic subroutine that hints about this. We can provide it in the documentation and instruct how modifications can be done.
4. I don't actually calculate the rotation angle but rather the components of the 2x2 rotation matrix R11 and R12. These two are sufficient as R11=R22=cos(theta) and R12=-R21=-sin(theta). So we save 2 arrays [cos(theta) and sin(theta)] instead of 1 angle array [theta], but avoid calculation of cosine and sine functions every time rotate_cyc is called. This is the original reason why I followed this approach. However, in GPU loading 2 arrays from global memory is gonna be slower than sine-cosine calculation so I think of going back to angles. Let me know if my reasoning is faulty.
5. I assigned one angle value for every GLL node on periodic faces. This keeps it more general, if the periodic faces are deformed for some reason (e.g. a wavy surface).
6. opr_device_rotate_cyc GPU implementation does not follow same approach as other subroutines in opr_device.f90 i.e. it does not copy arrays into shared memory before using them. We load and store them once so I believe shared memory is redundant.
7. Consequently, it is more similar to subroutines in math.f90 and mathops.f90, so we may consider putting it into another file.
8. Once we clarify these points, we can easily extend implementations to HIP and OpenCL.

Thanks!

[timfelle]

Great idea moving into the Frenet-Serret basis. But I don't quite see where you do that in the code, where do you estimate the center curve, and do you use bezier, b-splines or something completely different?

For a helical pipe you should also need the full 3D rotation, since your bi-normal would not be parallel to the z-axis, as your tangent would be pointing slightly up. (assuming a helicoid wound arround the z axis.)

[timofeymukha]

Very nice!

What one kind of wishes when one sees the code is to bake in the rotation into a new type of gs call for vector fields. So that one does not have to think about that "manually". However, that would need passing coef to gs, which currently will not work :P. So, this is not on you Isa, but rather the old can of works about the SEM type structure.

Isa, my understanding is that right now on the user side this is transparent, right? In the sense that as long as the needed transformation is supported, there is nothing else that needs to be done? I guess we need to put some description in the docs, and also have a warning when we detect that the transformation is not trivial, i.e. some rotation is applied, so that the user knows. Maybe an explicit user switch is good to raise awareness even further, as you write. I am not 100% sure.

[IsaMammadli]

@timofeymukha yes, once 'cyclic' in case file is set to true, nothing else needs to be done for the geometries mentioned above. Once PR goes through I can prepare the documentation. Because there are other small details e.g. only coupled CG solvers can be used and velocity projection must be switched off .

@timfelle very relevant points! This picture shows that we use surface normals of periodic faces and GLL nodes (x,y) (left) to derive T and N of TNB basis (right). This setup assumes that curvature is centered at (x,y)=(0,0) and it is wound around z axis. I will share some text with you for further explanation.