#350: Autogenerating cctk_Loop.h --------------------------+------------------------------------------------- Reporter: eschnett | Owner: Type: enhancement | Status: review Priority: minor | Milestone: Component: Cactus | Version: Resolution: | Keywords: --------------------------+-------------------------------------------------
Comment (by hinder):
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.