- changed status to open
- removed comment
Autogenerating cctk_Loop.h
Issue #350
closed
I wrote a perl script to auto-generate the Cactus header file Cactus/src/include/cctk_Loop.h.The script writes the code to stdout. I attach it to this ticket, as well as the output it generates.
I suggest to keep this perl script next to its output in the include directory, and to run it manually whenever required.
Keyword:
Comments (18)
-
reporter
-
reporter
- removed comment
- ping*
-
- removed comment
This is interesting with the scheduler changes in mind (because it would make writing those loops easier). Do you think it should be applied as is, or would you change something with that already in mind? Is there documentation to these new defines?
-
- removed comment
Having an easy-to-use looping facility in the flesh is a very good idea. However, this still leaves the task of finding the boundary sizes to the user. There are many functions around which do this, and I think we should provide one standard version. e.g. thorn Boundary could provide the required iblo, jblo, ibhi, jbhi variables for the interior and boundary, and could in fact provide macros similar to the ones in cctk_loop.h but which didn't require the user to find these sizes themselves.
The generated code looks OK (though very complicated, I assume this is unavoidable). Looks OK to apply.
-
reporter
- removed comment
I attach my current version of these files. There are generic loops CCTK_LOOP3, constructs with pre-defined sizes CCTK_ENDLOOP3_ALL, CCTK_LOOP3_INTERIOR, CCTK_LOOP3_BOUNDARIES, and also strided versions of each (CCTK_LOOP3STR etc.). The strided versions operate only on every nth grid point.
-
- removed comment
Carpet does not compile with this version:
/home/knarf/Cactus_CIGR/arrangements/Carpet/Carpet/src/Recompose.cc(2000): error: expected "while" } CCTK_ENDLOOP3_ALL(CarpetClassifyPoints); ^
-
- changed status to open
- removed comment
-
reporter
- changed status to open
- removed comment
This has the error above corrected.
It also introduces new macros CCTK_LOOP_INT, CCTK_LOOP_BND, and CCTK_LOOP_INTBND that do not require specifying the boundary widths -- these are obtained via an aliased function from thorn CoordBase.
-
- removed comment
Very nice! Getting this into the flesh means we can use it in lots of places, which should tidy up a number of things. I have a few comments/questions:
1. The lc_ prefix presumably comes from LoopControl; could this be CCTK_ instead? Since this is auto-generated code, such a change is probably easy to make in the Perl script. Is there a rule that users should not declare variables which start with CCTK_? Ideally we would use namespaces, but C does not support these.
2. I don't understand how OpenMP works with these loops.
3. The boundary loops go over all 26 boundary zones (6 faces, 12 edges and 8 corners). Is it possible to also provide a variable for the boundary condition code to find out which of these zones it is on? This is the "normal" provided by GenericFD. It is essentially an N-component array of face * dir. This could be called something like CCTK_NORMAL as a standard name, and declared by the looping macro, rather than requiring the user to provide such an array. This allows the boundary condition code to select the appropriate finite differencing operator. Thinking ahead, it is also useful to know how far away from the boundary you are in each direction. You can then reduce the stencil size point by point as you approach the boundary. The sign of the normal tells you whether you are on the lower or upper face. We could also use the magnitude of the components of the normal to represent the number of points less than the boundary size which are available. There might be better representations. I would like to be able to do something equivalent to this:
switch(CCTK_NORMAL[0]) { case 0: dfdx = D4(f); break;
case 1: case -1: dfdx = D2(f); break;
case 2: dfdx = D2Minus(f); break;
case -2: dfdx = D2Plus(f); break; } If we declare the normal inside the macro, rather than asking the user to provide an array, then we can add this later.
4. As far as I can tell, there are LOOP, LOOP_INTERIOR, LOOP_BOUNDARIES, LOOP_INTBOUNDARIES, LOOP_ALL, LOOP_INT, LOOP_BND, LOOP_INTBND, all defined as well for stride and dimension. 1. LOOP is the most basic, and requires you to specify the min and max coordinates of the grid. Could we (later) add some bounds checking code to ensure that the user doesn't loop outside the grid? 2. LOOP_INTERIOR loops over the interior, requiring you to tell it the boundary sizes. 3. LOOP_BOUNDARIES loops over all 26 boundary zones, requiring you to tell it the boundary sizes. 4. LOOP_INTBOUNDARIES: what does this do? Is this related to not looping over ghost zones? 5. LOOP_ALL loops over all points. 6. LOOP_INT is the same as LOOP_INTERIOR, but the boundary sizes are determined from the aliased function. 7. LOOP_BND and LOOP_INTBND similarly. The macro names seem a little arbitrary; it's not clear to the user that the difference between LOOP_BND and LOOP_BOUNDARIES relates to specifying additional information. Ideally we could just use LOOP_BOUNDARIES overloaded on the number of arguments, but I think the C preprocessor does not support this except via horrible tricks with variadic macros. Since the LOOP_BND case will be the most common, we could rename this to LOOP_BOUNDARIES (or LOOP_BOUNDARY? Should it be considered as a set of points or a list of zones?) and rename LOOP_BOUNDARIES as LOOP_BOUNDARIES_FROM_INDICES to indicate that you are going to provide the indices yourself. Similarly for the other macros.
5. I think that the BND macro loops over all non-interprocessor boundary zones; i.e. it loops over both physical and symmetry boundaries. In GenericFD, we only allow loops over the physical boundaries, determined by an additional array identifying whether the boundaries are symmetry boundaries or not, determined by an aliased function from SymBase. What do we want to do here? Do we want to allow the implementation of symmetry boundary conditions using this interface? Or should we just exclude symmetry boundaries from these loops? I think the latter - we are really targetting application code, and they should be concerned with physical boundaries only.
-
- removed comment
6. Optional: It would be convenient to be able to run the Perl script without having to remember to redirect to a specific file. i.e. the Perl script would write to the file.
-
reporter
- removed comment
1. Sure. 2. There is no OpenMP support. LoopControl can overload these macros, adding OpenMP support. We could add simple OpenMP support to this script as well. 3. For the boundaries, the user specifies three integers i,j,k for the current position, and three integers ni,nj,nk for the normal. This is currently only for INTBOUNDARIES; we can add this for BOUNDARIES as well. 4. INTBOUNDARIES loops over "internal boundaries", i.e. over those boundaries that are not also ghost or symmetry points. This allows taking tangential derivatives. I did not want to change the behaviour of the existing LOOP_INTERIOR and LOOP_BOUNDARIES. 5. Thanks -- this would not be intended, I'll look int this. 6. Can do.
-
- removed comment
2. I think it would be good to have OpenMP support in this, as then we can say that Cactus supports OpenMP out of the box, without adding the complexity of LoopControl, which is really designed to solve a much harder problem. 3. You mean variable names for each? The macro then declares variables named in this way for these things. What is the reason for letting the user choose these variable names, rather than providing a standard, such as CCTK_NORMAL? Maybe your way is better. I think it should be consistent for all the macros. The current implementation does not do everything that I said above; you cannot determine how far you are from the boundary. What do you think about this? 4. The usual case is to loop over "internal boundaries", so I think these should be called just "boundaries". Isn't it quite unusual to want to loop over ghost points and symmetry points? That version should have the less common name. Maybe we can find a way around the backwards compatibility issue. This code has never been released (it was committed after the release in November), and I suspect might not have been used. "internal boundaries" sounds like it loops over all ghost zones, as those are internal to the domain.
-
reporter
- removed comment
2. Ok.
3. If we use CCTK_NORMAL for the normal, then we should also use CCTK_INDEX for the loop index... Joke aside, we can store the distance from the outer boundary in the normal. (The "number of points less than the boundary size which are available" would be less useful -- it would be large at the boundary, and small in the interior, thus it wouldn't indicate what stencil you can use.)
4. The issue isn't looping over physics vs. looping over symmetry zones. The issue is that there are three boundary specifiers (physical, symmetry, ghost), and in three dimension, and boundary can have up to three specifies attach to it (corners!). Yes, a boundary can be all of a physical, a symmetry, and a ghost boundary. The issue is in handling these. BOUNDARIES loops over all boundaries that are physical boundaries (irrespective of whether they are also a symmetry or a ghost boundary), INTBOUNDARIES loops over all boundaries that are only physical boundaries. The former does not require applying symmetry conditions and synchronisation afterwards, the latter allows taking tangential derivatives.
-
- removed comment
Replying to [comment:13 eschnett]:
3. If we use CCTK_NORMAL for the normal, then we should also use CCTK_INDEX for the loop index...
Yes. Why not? Is the reason only due to brevity? Or is there another reason?
Joke aside, we can store the distance from the outer boundary in the normal. (The "number of points less than the boundary size which are available" would be less useful -- it would be large at the boundary, and small in the interior, thus it wouldn't indicate what stencil you can use.)
By "distance from the outer boundary", I assume you mean the distance from the outermost point in the boundary. But which one? Lower or upper? I thought about both options, and the most symmetric seemed to be the one I proposed. This would be of the form [-2, -1, 0, 0, ..., 0, 1, 2], assuming a boundary size of 2. You are correct that this array does not contain exactly the information about the stencil that you can use, but if you assume that the interior corresponds to a stencil of the same size as the boundary, then this is the missing information. Maybe this is too complicated. For your suggestion, I believe we would need two variables, one for the lower boundary and one for the upper boundary. This would then look like this: nx_lo = [0, 1, 2, ...]; nx_hi = [,. 2, 1, 0], then we could do
if (nx_lo[0] >= 2 && nx_hi[0] >= 2) { dfdx = D4(f); } else if (nx_lo[0] == 1 || nx_hi[0] == 1) { dfdx = D2(f); } else if (nx_lo[0] == 0) { dfdx = D2Plus(f); } else if (nx_hi[0] == 0) { dfdx = D2Minus(f); } else assert(0);
4. The issue isn't looping over physics vs. looping over symmetry zones. The issue is that there are three boundary specifiers (physical, symmetry, ghost), and in three dimension, and boundary can have up to three specifies attach to it (corners!). Yes, a boundary can be all of a physical, a symmetry, and a ghost boundary. The issue is in handling these. BOUNDARIES loops over all boundaries that are physical boundaries (irrespective of whether they are also a symmetry or a ghost boundary), INTBOUNDARIES loops over all boundaries that are only physical boundaries. The former does not require applying symmetry conditions and synchronisation afterwards, the latter allows taking tangential derivatives.
So BOUNDARIES would be useful if you were only taking derivatives normal to the boundary, as it would allow you to skip expensive interprocessor synchronisation and symmetry boundary conditions, and INTBOUNDARIES would be useful if you also needed to take tangential derivatives? This seems OK, but I wonder if we can come up with clearer names?
-
reporter
- removed comment
I suggest to set up the normal according to the distance from the boundary, plus one:
-1 -2 -3 0 0 0 0 0 +3 +2 +1 This is symmetric, and immediately useful for selecting a stencil. There is no need to determine and subtract the actual boundary width, nor to ensure that the boundary with is being set correctly in the parameter file.
-
reporter
- removed comment
Yes, brevity is the reason for not using CCTK_INDEX. This would require
CCTK_GFINDEX3D(cctkGH, CCTK_INDEX[0], CCTK_INDEX[1], CCTK_INDEX[2]) or (in Fortran)
gxx(CCTK_INDEX(1), CCTK_INDEX(2), CCTK_INDEX(3)) which may lead people to not use these macros (or make them add their own boilerplate to introduce shorter names).
-
reporter
- changed status to resolved
- removed comment
Applied.
-
- edited description
- changed status to closed
- Log in to comment
