Spring 2019
Jewett

Issues unique to the C language - that we will encounter in class - are discussed here.

Compiling

• We will generally use a Makefile when compiling.  But in general, to compile a code file with the Intel compiler icc:
  
       icc   name-of-your-program.c
  
  
• As we'll discuss in class, there are times when you need to employ additional checking of your code, to help track down possible errors (bugs).  In such cases, we will change the OPTIONS in our Makefile to enable additional code-checking when compiling; some of those compiler options are described here.  For example, we could do this:
  
      icc  -g -debug extended -traceback  name-of-your-program.c
  
  
• The Portland Group C compiler is available via the module command on Stampede.
  
       pgcc  name-of-your-program.c
  
  
• The GNU C compiler (free on all platforms)f:
  
       gcc  name-of-your-program.c
  
  

1-Boundary conditions and array layout

In C, 1-D arrays are usually defined/declared with statements like float s[nx].
This means that there are nx array elements, and we address them as s[0], s[1], ..., s[nx-1].

In some cases this is fine.  However, for arrays used in our integration, we will encounter issues connected to boundary conditions.

Imagine we have a domain that is nx grid points in size.  And, assume we use an array named s1 to hold the conditions at the present time, and array s2 to hold conditions in the future - at the next time step.  For some array index j, determining s2[j] - the future value of array s2 at index point j - requires knowing s1[j], the current value at point j.  However, most integration methods also require knowing the neighboring grid values: s1[j-1] and s1[j+1], for example - the grid points to the left and right of point j.

Why do we care?  Well -

• away from the boundaries - far from j=0 or j=(nx-1) - we can examine points s1[j-1] and s1[j+1] without problems.
• at the left boundary, s1[j=0], though, we would need the j-1 point:  s1[-1], which is off the left end of an array s1[nx].
  
• at the right boundary, s1[j=nx-1], we would need the j+1 point: s1[nx], which is off the right end of an array s1[nx].

Most codes deal with this one of two ways. 

1. One approach breaks down the integration procedure with checks to see if we are near a boundary.  For example,
   
       for ( j=0; j<nx; j++ ) {
           if (j == 0) then {
               do something special for integration at the left-boundary to get s2[j]
          } else if (j == (nx-1)) {
               do something special for integration at the right-boundary to get s2[j]
          } else {
               this is an "interior" grid point: do the usual integration, e.g. s2[j] = some function of s1[j-1], s1[j], and s1[j+1]
          }
       }
   
   This gets pretty messy, particularly in 2-D and 3-D; there are boundaries and corner points to consider.
   
   
2. Another approach, which we will follow: deal with the boundaries before you do the integration.

Because C does not allow variably-dimensioned boundaries (e.g. s1[-1:nx]), we have a few tricks to do.

• we need nx grid points, plus points on either side (in this case, the 'either side' boundary width is equal to 1)
• thus, here, we need nx+2 points
• we define new array starting and ending points to define the real, physical part of the array; call these I1 and I2
• the off-the-edge grid points are now I1-1 and I2+1
• the code now looks like this:

        bc( this subroutine sets the boundary points )
        for ( j=I1; j<=I2; j++ ) {
            do the same procedure for all j = I1, I1+1,...,I2-1,I2
        enddo

In other words, we set the grid points we need "off the ends" of our arrays before we do the integration.

We call those 'extra' grid points "ghost zones" or "ghost grid points" - they are there to help us solve things quickly and neatly (there are also some considerations for ghost or "halo" points for parallel computing we can talk about in class).

How we define and use arrays with boundary points when programming in C:

You will see something like the following in the program 1 source code:

#define NX 25
#define BC_WIDTH 1
#define I1 BC_WIDTH
#define I2 I1+NX-1
#define NXDIM NX+2*BC_WIDTH

float s[NXDIM] ...

• here, the number of extra-points-off-the-ends-of-the-array, BC_WIDTH, is = 1 (so we use points j-1, j, j+1).
• we need one 'extra grid point' on each end of the array, so the total array size, NXDIM, is NX + 2*BC_WIDTH
• our array is thus defined s1[NXDIM], which really means s1[nx+2] in our case.
• we use those NXDIM points in this way:
  
• s[0] is the left 'extra grid point' -- ghost point -- we refer to this as s[I1-1]
  
• s[I1] is the first physical or real grid point; I1=1 here; in general, I1 = BC_WIDTH
• s[I2] is the last physical or real grid point; I2=NXDIM-2 here.
• s[NXDIM-1] is the last (right side) 'extra grid point' -- ghost point -- we refer to this as =s[I2+1]
  

So you will see plenty of places in our C codes where extra array points are used, more than the simple [nx)]
In the demo code, we define arrays float s1[NXDIM], s2[NXDIM].  Only I1,...,I2 are physical; the others are 'helper' points.

Click on the following icon to see an image summarizing the above: