componentupdate.c

a full 3D simulation of electromagnetic waves with efficient absorbing boundary a full 3D simulation

/* Copyright notice:
   radarFDTD - A free 3D-FDTD simulation for EM-waves with GPML-ABCs ;-)
   Copyright (C) 2000 Carsten Aulbert
   ( 
   I used the following program as a 'manual', so in my opinion this needs to be quoted:
   ToyFDTD, version 1.02 
   The if-I-can-do-it-you-can-do-it FDTD! 
   Copyright (C) 1998,1999 Laurie E. Miller, Paul Hayes, Matthew O'Keefe 
   )

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

/* all inner updates come here */

/* Update the hx values: */

void UpdateInnerHx(void)
{
  int i,j,k;
  PRECISION nextY, nextZ; /* the upper neighbour-values in y and z direction */
  
  for(i=0; i<(nx); i++) 
    for(j=0; j<(ny); j++)
      for(k=0; k<(nz); k++)
	{
	 
	  if (j < ny-1)
	    nextY = ez[i][j+1][k];
	  else
	    nextY = pmlEZXTopXZ[i][0][k] + pmlEZYTopXZ[i][0][k];
	  
	  if (k < nz-1)
	    nextZ = ey[i][j][k+1];
	  else
	    nextZ = pmlEYXTopXY[i][j+pmlWidth][0] + pmlEYZTopXY[i][j+pmlWidth][0];
	  hx[i][j][k] = hx[i][j][k] + (materialConstants[material[i][j][k]][DTMUDZ]*(nextZ - ey[i][j][k]) - materialConstants[material[i][j][k]][DTMUDY]*(nextY - ez[i][j][k]));
	}
}


 /* Update the hy values: */
void UpdateInnerHy(void)
{
  int i,j,k;
  PRECISION nextX, nextZ; /*  the upper neighbour-values in y and z direction */

  for(i=0; i<(nx); i++)
    for(j=0; j<(ny); j++)
      for(k=0; k<(nz); k++)
	{
	  if (i < nx-1)
	    nextX = ez[i+1][j][k];
	  else
	    nextX = pmlEZXTopYZ[0][j+pmlWidth][k+pmlWidth] + pmlEZYTopYZ[0][j+pmlWidth][k+pmlWidth];
	  
	  if (k < nz-1)
	    nextZ = ex[i][j][k+1];
	  else
	    nextZ = pmlEXYTopXY[i][j+pmlWidth][0] + pmlEXZTopXY[i][j+pmlWidth][0];
	  
	  hy[i][j][k] = hy[i][j][k] + (materialConstants[material[i][j][k]][DTMUDX]*(nextX - ez[i][j][k]) - materialConstants[material[i][j][k]][DTMUDZ]*(nextZ - ex[i][j][k]));
	}
}

/* Update the hz values: */
void UpdateInnerHz(void)
{
  int i,j,k;
  PRECISION nextX, nextY;
  
  for(i=0; i<(nx); i++)
    for(j=0; j<(ny); j++)
      for(k=0; k<(nz); k++)
	{
	  if (i < nx-1)
	    nextX = ey[i+1][j][k];
	  else
	    nextX = pmlEYXTopYZ[0][j+pmlWidth][k+pmlWidth] +  pmlEYZTopYZ[0][j+pmlWidth][k+pmlWidth];
	  
	  if (j < ny-1)
	    nextY = ex[i][j+1][k];
	  else
	    nextY = pmlEXYTopXZ[i][0][k] + pmlEXZTopXZ[i][0][k];
	  
	  hz[i][j][k] = hz[i][j][k] + (materialConstants[material[i][j][k]][DTMUDY]*(nextY - ex[i][j][k]) - materialConstants[material[i][j][k]][DTMUDX]*(nextX - ey[i][j][k]));
	}
}

/* Update the ex values: */
void UpdateInnerEx(void)
{
  int i,j,k;
  PRECISION lowerY, lowerZ;
  
  for(i=0; i<(nx); i++)
    for(j=0; j<(ny); j++) 
      for(k=0; k<(nz); k++) 
	{
	  if (j > 0)
	    lowerY = hz[i][j-1][k];
	  else
	    lowerY = pmlHZXBottomXZ[i][pmlWidth-1][k] + pmlHZYBottomXZ[i][pmlWidth-1][k];
	  
	  if (k > 0)
	    lowerZ = hy[i][j][k-1];
	  else
	    lowerZ = pmlHYXBottomXY[i][j+pmlWidth][pmlWidth-1] + pmlHYZBottomXY[i][j+pmlWidth][pmlWidth-1];
	  
	  switch(material[i][j][k])
	    {
	    case 0: 
	      ex[i][j][k] = ex[i][j][k] + materialConstants[0][DTEPSDY] * (hz[i][j][k] - lowerY) - materialConstants[0][DTEPSDZ] * (hy[i][j][k] - lowerZ);
	      break;
	    case 1:
	      ex[i][j][k] = 0;
	      break;
	    default:
	      ex[i][j][k] = materialConstants[material[i][j][k]][DTSIGEPS]*ex[i][j][k] + materialConstants[material[i][j][k]][DTEPSDY]*(hz[i][j][k] - lowerY) - materialConstants[material[i][j][k]][DTEPSDZ]*(hy[i][j][k] - lowerZ);
	    }
	}
}

/* Update the ey values: */
  
void UpdateInnerEy(void)
{
  int i,j,k;
  PRECISION lowerX, lowerZ;
  
  for(i=0; i<(nx); i++) 
    for(j=0; j<(ny); j++)
      for(k=0; k<(nz); k++) 
	{
	  if (i > 0)
	    lowerX = hz[i-1][j][k];
	  else
	    lowerX = pmlHZXBottomYZ[pmlWidth-1][j+pmlWidth][k+pmlWidth] + pmlHZYBottomYZ[pmlWidth-1][j+pmlWidth][k+pmlWidth];
	  
	  if (k > 0)
	    lowerZ = hx[i][j][k-1]; 
	  else
	    lowerZ = pmlHXYBottomXY[i][j+pmlWidth][pmlWidth-1] +  pmlHXZBottomXY[i][j+pmlWidth][pmlWidth-1];
	  
	  switch(material[i][j][k])
	    {
	    case 0: 
	      ey[i][j][k] = ey[i][j][k] + materialConstants[0][DTEPSDZ] * (hx[i][j][k] - lowerZ) - materialConstants[0][DTEPSDX] * (hz[i][j][k] - lowerX);
	      break;
	    case 1:
	      ey[i][j][k] = 0;
	      break;
	    default:
	      ey[i][j][k] = materialConstants[material[i][j][k]][DTSIGEPS]*ey[i][j][k] + materialConstants[material[i][j][k]][DTEPSDZ]*(hx[i][j][k] - lowerZ) - materialConstants[material[i][j][k]][DTEPSDX]*(hz[i][j][k] - lowerX);
	    }
	}
}
 
/* Update the ez values: */
void UpdateInnerEz(void)
{
  int i,j,k;
  PRECISION lowerX, lowerY;

  for(i=0; i<(nx); i++) 
    for(j=0; j<(ny); j++)
      for(k=0; k<(nz); k++)
	{
	  if (i > 0)
	    lowerX = hy[i-1][j][k];
	  else
	    lowerX = pmlHYXBottomYZ[pmlWidth-1][j+pmlWidth][k+pmlWidth] + pmlHYZBottomYZ[pmlWidth-1][j+pmlWidth][k+pmlWidth];
	  
	  if (j > 0)
	    lowerY = hx[i][j-1][k];
	  else
	    lowerY = pmlHXYBottomXZ[i][pmlWidth-1][k] + pmlHXZBottomXZ[i][pmlWidth-1][k];

	  switch(material[i][j][k])
	    {
	    case 0: 
	      ez[i][j][k] = ez[i][j][k] + materialConstants[0][DTEPSDX] * (hy[i][j][k] - lowerX) - materialConstants[0][DTEPSDY] * (hx[i][j][k] - lowerY);
	      break;
	    case 1:
	      ez[i][j][k] = 0;
	      break;
	    default:
	      ez[i][j][k] = materialConstants[material[i][j][k]][DTSIGEPS]*ez[i][j][k] + materialConstants[material[i][j][k]][DTEPSDX]*(hy[i][j][k] - lowerX) - materialConstants[material[i][j][k]][DTEPSDY]*(hx[i][j][k] - lowerY);
	    }
	}
}

/* THIS IS THE PML SECTION */

void UpdatePMLHx(void)
{
  int i,j,k;
  PRECISION mudt = MU_0 / dt;
  PRECISION topY1, topY2, topZ1, topZ2; /* components which may be in different memory blocks than the one worked on */
  PRECISION commonFactor;
  
  /* speed up variables */
  int ny_plus_2pmlWidth = ny + 2*pmlWidth;
  int nz_plus_2pmlWidth = nz + 2*pmlWidth;

  int ny_plus_pmlWidth  = ny + pmlWidth;

  int pmlWidth_minus_1  = pmlWidth-1;

  int nz_minus_1        = nz-1;
  
  /* start with the YZ planes */

  /* there are no problems at any boundary ! */
  for(i=0; i<pmlWidth; i++)
    for(j=0; j<ny_plus_2pmlWidth; j++)
      for(k=0; k<nz_plus_2pmlWidth; k++)
	{ /* bottom */
	    
	    commonFactor =  mudt + 0.5 * pmlSigmaYStarBottomYZ[i][j][k];
	    pmlHXYBottomYZ[i][j][k] = pmlHXYBottomYZ[i][j][k] * 
	       (mudt - 0.5 * pmlSigmaYStarBottomYZ[i][j][k]) / commonFactor -
	      (pmlEZXBottomYZ[i][j+1][k] - pmlEZXBottomYZ[i][j][k] + 
	       pmlEZYBottomYZ[i][j+1][k] - pmlEZYBottomYZ[i][j][k]) / 
	      (dy * pmlSYBottomYZ[i][j][k] * commonFactor);
	    
	    commonFactor = mudt + 0.5 * pmlSigmaZStarBottomYZ[i][j][k];
	    pmlHXZBottomYZ[i][j][k] = pmlHXZBottomYZ[i][j][k] * 
	       (mudt - 0.5 * pmlSigmaZStarBottomYZ[i][j][k]) / commonFactor +
	      (pmlEYZBottomYZ[i][j][k+1] - pmlEYZBottomYZ[i][j][k] + 
	       pmlEYXBottomYZ[i][j][k+1] - pmlEYXBottomYZ[i][j][k]) /
	      (dz * pmlSZBottomYZ[i][j][k] * commonFactor);
	}
  /* top, there are no problems at any boundary */
  for(i=0; i<pmlWidth; i++)
    for(j=0; j<ny_plus_2pmlWidth; j++)
      for(k=0; k<nz_plus_2pmlWidth; k++)
	{
	    commonFactor =  mudt + 0.5 * pmlSigmaYStarBottomYZ[i][j][k];
	    pmlHXYTopYZ[i][j][k] = pmlHXYTopYZ[i][j][k] * 
	       (mudt - 0.5 * pmlSigmaYStarTopYZ[i][j][k]) / commonFactor -
	      (pmlEZXTopYZ[i][j+1][k] - pmlEZXTopYZ[i][j][k] + 
	       pmlEZYTopYZ[i][j+1][k] - pmlEZYTopYZ[i][j][k]) / 
	      (dy * pmlSYTopYZ[i][j][k] * commonFactor);
	    
	    commonFactor = mudt + 0.5 * pmlSigmaZStarTopYZ[i][j][k];
	    pmlHXZTopYZ[i][j][k] = pmlHXZTopYZ[i][j][k] * 
	       (mudt - 0.5 * pmlSigmaZStarTopYZ[i][j][k]) / commonFactor +
	      (pmlEYZTopYZ[i][j][k+1] - pmlEYZTopYZ[i][j][k] + 
	       pmlEYXTopYZ[i][j][k+1] - pmlEYXTopYZ[i][j][k]) /
	      (dz * pmlSZTopYZ[i][j][k] * commonFactor);

	}

  /* now the XY planes */

  for(i=0; i<nx; i++)
    for(j=0; j<ny_plus_2pmlWidth; j++)
      for(k=0; k<pmlWidth; k++)
	{ /* bottom */
	  /* there are three boundary-problems:
	     0           <= y < pmlWidth:     XZ - plane (bottom)
	     pmlWidth    <= y < pmlWidth+ny:  SimSpace
	     pmlWidth+ny <= y:                XZ - plane (top)
	  */

	  if (k < pmlWidth_minus_1)
	    {
	      topZ1 = pmlEYZBottomXY[i][j][k+1];
	      topZ2 = pmlEYXBottomXY[i][j][k+1];
	    }
	  else
	    {
	      if (j < pmlWidth)
		{
		  topZ1 = pmlEYZBottomXZ[i][j][0];
		  topZ2 = pmlEYXBottomXZ[i][j][0];
		}
	      else if (j < ny_plus_pmlWidth)
		{
		  topZ1 = ey[i][j-pmlWidth][0];
		  topZ2 = 0.0;
		}
	      else
		{
		  topZ1 = pmlEYZTopXZ[i][j-ny_plus_pmlWidth][0];
		  topZ2 = pmlEYXTopXZ[i][j-ny_plus_pmlWidth][0];
		}
	    }

	  commonFactor =  mudt + 0.5 * pmlSigmaYStarBottomXY[i][j][k];
	  pmlHXYBottomXY[i][j][k] = pmlHXYBottomXY[i][j][k] * 
	    (mudt - 0.5 * pmlSigmaYStarBottomXY[i][j][k]) / commonFactor -
	    (pmlEZXBottomXY[i][j+1][k] - pmlEZXBottomXY[i][j][k] + 
	     pmlEZYBottomXY[i][j+1][k] - pmlEZYBottomXY[i][j][k]) / 
	    (dy * pmlSYBottomXY[i][j][k] * commonFactor);
	  
	  commonFactor = mudt + 0.5 * pmlSigmaZStarBottomXY[i][j][k];
	  pmlHXZBottomXY[i][j][k] = pmlHXZBottomXY[i][j][k] * 
	    (mudt - 0.5 * pmlSigmaZStarBottomXY[i][j][k]) / commonFactor +
	    (topZ1 - pmlEYZBottomXY[i][j][k] + 
	     topZ2 - pmlEYXBottomXY[i][j][k]) /
	    (dz * pmlSZBottomXY[i][j][k] * commonFactor);
	  
	}
  
  for(i=0; i<nx; i++)
    for(j=0; j<ny_plus_2pmlWidth; j++)
      for(k=0; k<pmlWidth; k++)
	{ 
	  /* top */
	  /* there are no problems at this boundary */

	  commonFactor =  mudt + 0.5 * pmlSigmaYStarTopXY[i][j][k];
	  pmlHXYTopXY[i][j][k] = pmlHXYTopXY[i][j][k] * 
	    (mudt - 0.5 * pmlSigmaYStarTopXY[i][j][k]) / commonFactor -
	    (pmlEZXTopXY[i][j+1][k] - pmlEZXTopXY[i][j][k] + 
	     pmlEZYTopXY[i][j+1][k] - pmlEZYTopXY[i][j][k]) / 
	    (dy * pmlSYTopXY[i][j][k] * commonFactor);
	  
	  commonFactor = mudt + 0.5 * pmlSigmaZStarTopXY[i][j][k];
	  pmlHXZTopXY[i][j][k] = pmlHXZTopXY[i][j][k] * 
	    (mudt - 0.5 * pmlSigmaZStarTopXY[i][j][k]) / commonFactor +
	    (pmlEYZTopXY[i][j][k+1] - pmlEYZTopXY[i][j][k] + 
	     pmlEYXTopXY[i][j][k+1] - pmlEYXTopXY[i][j][k]) /
	    (dz * pmlSZTopXY[i][j][k] * commonFactor);	  
	} 
  
  /* and finally the XZ planes */
  for(i=0; i<nx; i++)
    for(j=0; j<pmlWidth; j++)
      for(k=0; k<nz; k++)
	{ /* bottom */
	  /* there are boundaries at:
	     j = pmlWidth-1 for all i,k      SimSpace is neighbour
	     k = ny_plus_pmlWidth-1 for all i,j   TopXY is neighbour 
	  */

	  if (j < pmlWidth_minus_1)
	    {
	  	topY1 = pmlEZXBottomXZ[i][j+1][k];
		topY2 = pmlEZYBottomXZ[i][j+1][k];
	    }
	  else
	    {
	      topY1 = ez[i][0][k];
	      topY2 = 0.0;
	    }
	  
	  if (k < nz_minus_1)
	    {
	      topZ1 = pmlEYZBottomXZ[i][j][k+1];
	      topZ2 = pmlEYXBottomXZ[i][j][k+1];
	    }
	  else
	    {
	      topZ1 = pmlEYZTopXY[i][j][0];
	      topZ2 = pmlEYXTopXY[i][j][0];
	    }
	  
	  commonFactor =  mudt + 0.5 * pmlSigmaYStarBottomXZ[i][j][k];
	  pmlHXYBottomXZ[i][j][k] = pmlHXYBottomXZ[i][j][k] * 
	    (mudt - 0.5 * pmlSigmaYStarBottomXZ[i][j][k]) / commonFactor -
	    (topY1 - pmlEZXBottomXZ[i][j][k] + 
	     topY2 - pmlEZYBottomXZ[i][j][k]) /