Fix FV operators in double null configurations #2626
Conversation
Used mesh::firstY() and mesh::lastY() rather than mesh::firstY(xind) and mesh::lastY(xind). Those only check the global Y index, and so don't check for boundaries in the middle of the Y domain. This causes a problem for double null tokamak simulations.
The versions of these functions without X index don't actually give any useful information regarding mesh topology. They give the correct result in many cases but not all. Instead the versions with X index argument should be used.
|
clang-tidy review says "All clean, LGTM! 👍" |
There was a problem hiding this comment.
When using firstY(ix) and lastY(ix), is it also necessary to check periodicY(ix) the way D4DY4_Index() does?
Lines 363 to 366 in 3bf597b
FV::Div_par does something similar
BOUT-dev/include/bout/fv_ops.hxx
Line 220 in 3bf597b
There was something in the comments for firstY(ix)/lastY(ix) about just being the first/last in the communicator?
src/mesh/fv_ops.cxx
Outdated
| @@ -325,7 +325,7 @@ Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f) { | |||
| result(i, j, k) += flux / (coord->J(i, j, k) * coord->dy(i, j, k)); | |||
| result(i, j + 1, k) -= flux / (coord->J(i, j + 1, k) * coord->dy(i, j + 1, k)); | |||
|
|
|||
| if(j == mesh->ystart && (!mesh->firstY())) { | |||
| if(j == mesh->ystart && (!mesh->firstY(i))) { | |||
There was a problem hiding this comment.
Why does this function not need any check for the upper boundary (mesh->lastY(i)), when it does for the lower boundary?
There was a problem hiding this comment.
Thanks @johnomotani ! I've rewritten this function to be more consistent in handling the boundaries
Thanks to @johnomotani! The periodicY(x) value should be checked when using firstY / lastY
|
clang-tidy review says "All clean, LGTM! 👍" |
The
Mesh::firstY()andMesh::lastY()functions can be used in slab and single-null tokamak geometry to determine if a cell is at a Y boundary. In double null geometries this fails because the upper target boundaries are in the middle of the Y domain.The versions of these functions taking an X index use the communicators to determine if a point is first or last on its flux surface.
This fix replaces calls to no-argument
firstY()andlastY()with calls the versions taking an X index. It also deprecates the no-argument functions as they are likely to cause similar confusion in future.