Handling Boundary Conditions
============================

Author: Tom Goodale
Date: 16 Dec 2002
Status: Fixed typos, explicit spec of arg lists, clarified some points.

This discussion is about a generic method for people to share boundary
conditions, to have run-time selection of boundary conditions, and to have
symmetry bcs interact cleanly with other bcs without people having to
hard-code things like

if(bc = "radiative")
  call radiative_bc
elsif
  ...
else
  ...
fi

#ifdef CARTOON2D
if(symm = "cartoon")
  call cartoon
fi
#endif

which is ugly and means you need to know about every possible bc when you 
write your code.

We have two sorts of boundary conditions, in general, symmetry and 
physical, and each of these can be local, i.e. based upon some
point-wise operation from neighbouring points, or non-local,
i.e. based upon some data which is not in the neighbourhood of the
boundary point.

Symmetry BCs apply some symmetry condition, such as periodic or
rotational, to whichever variables are having boundary conditions
applied, and are, in general, non-local, 'though things like the
Cartoon2D boundary can be implemented as a local boundary condition.
In principle many different symmetry conditions could be active at
once.  Physics codes should be independent of which symmetries there
are, as these are properties of the grid.

Physical BCs apply a specific boundary condition to a variable.  Each
variable only has one physical boundary condition at any boundary point,
'though, in principle, a different one could be applied at each boundary
point.  In practice the most common case is treating the whole boundary
the same or the whole of particular faces the same.  Physical boundary
conditions are, in general, local boundary conditions, though there could
also be a physical boundary condition which involves solving an elliptic
equation across the boundary.

To summarise the difference between the two:  there is at most one 
physical boundary condition at any point, whilst there may be many 
symmetry boundary conditions at a point.  An evolution scheme only knows 
about the physical boundary conditions which must be applied, not the 
symmetry boundary conditions.

So I'd like to propose the following scheme.

We have an implementation Boundary which provides the following aliased 
functions:

Boundary_RegisterPhysicalBC(const cGH *cctkGH, 
                            CCTK_FPOINTER function_pointer, 
                            const char * bc_name)

where function pointer is a pointer to a function which returns an int
and takes arguments

(const cGH *, int size, int *indices, int *table_handles)

where indices are the (size) indices of the GVs requiring the bc, and
table_handles is an array of (size) associated tables containing specific
boundary condition information for the corresponding GV.  These arrays are
sorted first on table handle and then on variable index so that variables
with precisely the same bc are grouped together, and within this grouping
variables' groups are similarly grouped together.

Boundary_SelectVarForBC(const cGH *cctkGH, 
                        int table_handle, 
                        const char *var_name, 
                        const char *bc_name)

Boundary_SelectVarForBCI(const cGH *cctkGH, 
                        int table_handle, 
                        int var_index, 
                        const char *bc_name)

A thorn can request that boundary condition bc_name gets called on a
grid variable.  The grid variable can be specified either by its
variable name or index.

Boundary_SelectGroupForBC(const cGH *cctkGH, 
                          int table_handle, 
                          const char *group_name, 
                          const char *bc_name)

Boundary_SelectGroupForBCI(const cGH *cctkGH, 
                           int table_handle, 
                           int group_index, 
                           const char *bc_name)

For convenience, an entire group of grid variables can be selected at
once.  Of course each variable within the group must use the same
table handle.

int Boundary_SelectedGVs(const cGH *cctkGH, 
                         int array_size,
                         int *indices,
                         int *tables,
                         const char *bc_name)

This function places a list of (array_size) grid variable indices, sorted
in the manner expected by the function implementing the boundary
condition, which are selected for 'bc_name', into the array 'indices'.  
The corresponding (array_size) table handles are placed into the array
'tables'.  It returns the total number of grid variables selected for
boundary condition 'bc_name'.  Thus to find out how big array_size is the
user may call the routine first with an array_size of 0.
If bc_name is passed as NULL, all selected variables (for any boundary
condition) are returned.

The implementation Boundary also creates two schedule groups

ApplyBCs

BoundaryConditions in ApplyBCs BEFORE Boundary_ClearSelection

It schedules two functions

Boundary_ApplyPhysicalBCs IN BoundaryConditions

Boundary_ClearSelection IN ApplyBCs

in the schedule groups.

Then each thorn which provides a physical boundary condition registers a
routine for its bc with the Boundary_RegisterPhysicalBC routine.

Each thorn providing a non-local physical bc additionally schedules a
routine in the BoundaryConditions schedule group.  It should register NULL
as its function in the function registry.  (Implementation note:  for
efficiency this should probably be turned, internally, into registration
of a function which does nothing, thus removing an if from the main
program flow which would be redundant in the most common path.)

Each thorn providing a symmetry bc schedules a routine in the 
BoundaryConditions schedule group.  Note that consistency of the symmetry 
conditions will be dealt with separately by the upcoming revised symmetry 
and coordinate stuff.

When a thorn needs to apply a boundary condition to a function it calls 
Boundary_Select*ForBC* with the appropriate info.  Then, at the right place 
for BCs to be applied, it schedules the ApplyBCs group.

The routines providing symmetry bcs and non-local physical bcs then get the
list of marked GVs and apply the bc to them.

The Boundary_ApplyPhysicalBCs routine calls the registered routines for
the local bcs.

The Boundary_ClearSelection clears the selection list.

(At the moment the GHs in the argument list are ignored, this is there
so that the interface doesn't change in 4.1.)

This system may be revised in 4.1 if certain modifications to the 
scheduler and parameter systems make a more efficient, simpler, or easier 
system possible.  See previous version of this spec and mails on the 
mailing list for some such alternatives.