abc2d.c

for applying a second-order ABC to a 1D grid

/* package for applying a second-order ABC to a 1D grid. */
#include <stdlib.h>

#include "fdtd-general.h"  // has the allocation macro

#define Ez(M,N) ez[(M)*size_y +(N)]

/* Define macros for arrays which store the previous values of the
 * fields.  For each one of these arrays the three arguments are as
 * follows:
 *
 *   first argument: spatial displacement from the boundary
 *   second argument: displacement back in time
 *   third argument: distance from either bottom or left side of grid
 *                   (depending on orientation)
 */
#define EzLeft(M,Q,N)   ezLeft[  (N)*6 + (Q)*3 + (M)]
#define EzRight(M,Q,N)  ezRight[ (N)*6 + (Q)*3 + (M)]
#define EzTop(N,Q,M)    ezTop[   (M)*6 + (Q)*3 + (N)]
#define EzBottom(N,Q,M) ezBottom[(M)*6 + (Q)*3 + (N)]

/* global variables not visible outside of this package */
static int size_x, size_y; 
static double *ez, s;
static double *ezLeft, *ezRight, *ezTop, *ezBottom;

/* initialization function */
void abc2d_init(double *the_ez,
		int the_size_x, int the_size_y,
		double the_sc_local)
{
  /* record arguments in local static variables */
  ez = the_ez;
  size_x = the_size_x;
  size_y = the_size_y;
  s = the_sc_local;

  /* allocate memory for ABC arrays */
  ALLOC_1D(ezLeft,size_y*6,double);
  ALLOC_1D(ezRight,size_y*6,double);
  ALLOC_1D(ezTop,size_x*6,double);
  ALLOC_1D(ezBottom,size_x*6,double);

  return;
} /* end abc_init() */

/* function which applies ABC -- called once per time step */
void abc2d()
{
  int m, n;

  /* ABC at left side of grid */
  for (n=0; n<size_y; n++) {
    Ez(0,n) = -((1./s-2.+s)*(Ez(2,n)+EzLeft(0,0,n))
     + 2.*(s-1./s)*(EzLeft(0,1,n)+EzLeft(2,1,n)
		    -Ez(1,n)-EzLeft(1,0,n))
     - 4.*(1./s+s)*(EzLeft(1,1,n)))/(1./s+2.+s)-EzLeft(2,0,n);

    /* memorize old fields */ 
    for (m=0; m<3; m++) {
      EzLeft(m,0,n)=EzLeft(m,1,n);
      EzLeft(m,1,n)=Ez(m,n);
    }
  }
  
  /* ABC at right side of grid */
  for (n=0; n<size_y; n++) {
    Ez(size_x-1,n) = -((1./s-2.+s)*(Ez(size_x-3,n)+EzRight(0,0,n))
     + 2.*(s-1./s)*(EzRight(0,1,n)+EzRight(2,1,n)
		    -Ez(size_x-2,n)-EzRight(1,0,n))
     -4.*(1./s+s)*(EzRight(1,1,n)))/(1./s+2.+s)-EzRight(2,0,n);
  
    /* memorize old fields */
    for (m=0; m<3; m++) {
      EzRight(m,0,n)=EzRight(m,1,n);
      EzRight(m,1,n)=Ez(size_x-1-m,n);
    }
  }

  /* ABC at bottom of grid */
  for (m=0; m<size_x; m++) {
    Ez(m,0) = -((1./s-2.+s)*(Ez(m,2)+EzBottom(0,0,m))
     + 2.*(s-1./s)*(EzBottom(0,1,m)+EzBottom(2,1,m)
		    -Ez(m,1)-EzBottom(1,0,m))
     - 4.*(1./s+s)*(EzBottom(1,1,m)))/(1./s+2.+s)-EzBottom(2,0,m);

    /* memorize old fields */ 
    for (n=0; n<3; n++) {
      EzBottom(n,0,m)=EzBottom(n,1,m);
      EzBottom(n,1,m)=Ez(m,n);
    }
  }
  
  /* ABC at top of grid */
  for (m=0; m<size_x; m++) {
    Ez(m,size_y-1) = -((1./s-2.+s)*(Ez(m,size_y-3)+EzTop(0,0,m))
     + 2.*(s-1./s)*(EzTop(0,1,m)+EzTop(2,1,m)
		    -Ez(m,size_y-2)-EzTop(1,0,m))
     -4.*(1./s+s)*(EzTop(1,1,m)))/(1./s+2.+s)-EzTop(2,0,m);
  
    /* memorize old fields */
    for (n=0; n<3; n++) {
      EzTop(n,0,m)=EzTop(n,1,m);
      EzTop(n,1,m)=Ez(m,size_y-1-n);
    }
  }

  return;
} /* end abc() */