componentupdate.c

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

	    (dy * pmlSYBottomXZ[i][j][k] * commonFactor);
	  
	  commonFactor = mudt + 0.5 * pmlSigmaZStarBottomXZ[i][j][k];
	  pmlHXZBottomXZ[i][j][k] = pmlHXZBottomXZ[i][j][k] * 
	    (mudt - 0.5 * pmlSigmaZStarBottomXZ[i][j][k]) / commonFactor +
	    (topZ1 - pmlEYZBottomXZ[i][j][k] + 
	     topZ2 - pmlEYXBottomXZ[i][j][k]) /
	    (dz * pmlSZBottomXZ[i][j][k] * commonFactor);

	}
  
  for(i=0; i<nx; i++)
    for(j=0; j<pmlWidth; j++)
      for(k=0; k<nz; k++)
	{  /* top */
	  /* Here you have only one boundary:
	     k = nz_minus_1 for all i,j          TopXY is neighbour */
	  if (k < nz_minus_1)
	    {
	      topZ1 = pmlEYZTopXZ[i][j][k+1];
	      topZ2 = pmlEYXTopXZ[i][j][k+1];
	    }
	  else
	    {
	      topZ1 = pmlEYZTopXY[i][j+ny_plus_pmlWidth][0];
	      topZ2 = pmlEYXTopXY[i][j+ny_plus_pmlWidth][0];
	    }
	  
	  commonFactor =  mudt + 0.5 * pmlSigmaYStarTopXZ[i][j][k];
	  pmlHXYTopXZ[i][j][k] = pmlHXYTopXZ[i][j][k] * 
	    (mudt - 0.5 * pmlSigmaYStarTopXZ[i][j][k]) / commonFactor -
	    (pmlEZXTopXZ[i][j+1][k] - pmlEZXTopXZ[i][j][k] + 
	     pmlEZYTopXZ[i][j+1][k] - pmlEZYTopXZ[i][j][k]) / 
	    (dy * pmlSYTopXZ[i][j][k] * commonFactor);
	  
	  commonFactor = mudt + 0.5 * pmlSigmaZStarTopXZ[i][j][k];
	  pmlHXZTopXZ[i][j][k] = pmlHXZTopXZ[i][j][k] * 
	    (mudt - 0.5 * pmlSigmaZStarTopXZ[i][j][k]) / commonFactor +
	    (topZ1 - pmlEYZTopXZ[i][j][k] + 
	     topZ2 - pmlEYXTopXZ[i][j][k]) /
	    (dz * pmlSZTopXZ[i][j][k] * commonFactor);
	} 
  /* that was easy, wasn磘 it? */
}

void UpdatePMLHy(void)
{
  int i,j,k;
  PRECISION mudt = MU_0 / dt;
  PRECISION topX1, topX2, 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 nz_plus_pmlWidth  = nz + pmlWidth;

  int pmlWidth_minus_1  = pmlWidth-1;

  int nx_minus_1        = nx-1;
  int nz_minus_1        = nz-1;
  
  
  /* starting with YZ - plane (bottom) */

  for(i=0; i<pmlWidth; i++)
    for(j=0; j<ny_plus_2pmlWidth; j++)
      for(k=0; k<nz_plus_2pmlWidth; k++)
	{
	  if (i < pmlWidth_minus_1)
	    {
	      topX1 = pmlEZXBottomYZ[i+1][j][k];
	      topX2 = pmlEZYBottomYZ[i+1][j][k];
	    }
	  else
	    {
	      if (k < pmlWidth)
		{
		  topX1 = pmlEZXBottomXY[0][j][k];
		  topX2 = pmlEZYBottomXY[0][j][k];
		}
	      else if (k < nz_plus_pmlWidth)
		{
		  if (j < pmlWidth)
		    {
		      topX1 = pmlEZXBottomXZ[0][j][k-pmlWidth];
		      topX2 = pmlEZYBottomXZ[0][j][k-pmlWidth];
		    }
		  else if (j < ny_plus_pmlWidth)
		    {
		      topX1 = ez[0][j-pmlWidth][k-pmlWidth];
		      topX2 = 0.0;
		    }
		  else
		    {
		      topX1 = pmlEZXTopXZ[0][j-ny_plus_pmlWidth][k-pmlWidth];
		      topX2 = pmlEZYTopXZ[0][j-ny_plus_pmlWidth][k-pmlWidth];
		    }
		}
	      else 
		{
		  topX1 = pmlEZXTopXY[0][j][k-nz_plus_pmlWidth];
		  topX2 = pmlEZYTopXY[0][j][k-nz_plus_pmlWidth];
		}
	    }
	 
	  commonFactor = mudt + 0.5 * pmlSigmaZStarBottomYZ[i][j][k];
	  pmlHYZBottomYZ[i][j][k] = pmlHYZBottomYZ[i][j][k] *
	    (mudt - 0.5 * pmlSigmaZStarBottomYZ[i][j][k]) / commonFactor -
	    (pmlEXYBottomYZ[i][j][k+1] - pmlEXYBottomYZ[i][j][k] +
	     pmlEXZBottomYZ[i][j][k+1] - pmlEXZBottomYZ[i][j][k]) /
	    (dz * pmlSZBottomYZ[i][j][k] * commonFactor);
	  
	  commonFactor = mudt + 0.5 * pmlSigmaXStarBottomYZ[i][j][k];
	  pmlHYXBottomYZ[i][j][k] = pmlHYXBottomYZ[i][j][k] *
	     (mudt - 0.5 * pmlSigmaXStarBottomYZ[i][j][k])  / commonFactor +
	    (topX1 - pmlEZXBottomYZ[i][j][k] +
	     topX2 - pmlEZYBottomYZ[i][j][k]) /
	    (dx * pmlSXBottomYZ[i][j][k] * commonFactor);
	}

  /* and the top part of it */

  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 * pmlSigmaZStarTopYZ[i][j][k];
	  pmlHYZTopYZ[i][j][k] = pmlHYZTopYZ[i][j][k] *
	    (mudt - 0.5 * pmlSigmaZStarTopYZ[i][j][k]) / commonFactor -
	    (pmlEXYTopYZ[i][j][k+1] - pmlEXYTopYZ[i][j][k] +
	     pmlEXZTopYZ[i][j][k+1] - pmlEXZTopYZ[i][j][k]) /
	    (dz * pmlSZTopYZ[i][j][k] * commonFactor);
	  
	  commonFactor = mudt + 0.5 * pmlSigmaXStarTopYZ[i][j][k];
	  pmlHYXTopYZ[i][j][k] = pmlHYXTopYZ[i][j][k] *
	     (mudt - 0.5 * pmlSigmaXStarTopYZ[i][j][k])  / commonFactor +
	    (pmlEZXTopYZ[i+1][j][k] - pmlEZXTopYZ[i][j][k] +
	     pmlEZYTopYZ[i+1][j][k] - pmlEZYTopYZ[i][j][k]) /
	    (dx * pmlSXTopYZ[i][j][k] * commonFactor);
	}

  /* starting with bottomXY plane */

  for(i=0; i<nx; i++)
    for(j=0; j<ny_plus_2pmlWidth; j++)
      for(k=0; k<pmlWidth; k++)
	{
	  if (k < pmlWidth_minus_1)
	    {
	      topZ1 = pmlEXYBottomXY[i][j][k+1];
	      topZ2 = pmlEXZBottomXY[i][j][k+1];
	    }
	  else
	    {
	      if (j < pmlWidth)
		{
		  topZ1 = pmlEXYBottomXZ[i][j][0];
		  topZ2 = pmlEXZBottomXZ[i][j][0];
		}
	      else if (j < ny_plus_pmlWidth)
		{
		  topZ1 = ex[i][j-pmlWidth][0];
		  topZ2 = 0.0;
		}
	      else
		{
		  topZ1 = pmlEXYTopXZ[i][j-ny_plus_pmlWidth][0];
		  topZ2 = pmlEXZTopXZ[i][j-ny_plus_pmlWidth][0];
		}
	    }

	  if (i <  nx_minus_1)
	    {
	      topX1 = pmlEZXBottomXY[i+1][j][k];
	      topX2 = pmlEZYBottomXY[i+1][j][k];
	    }
	  else
	    {
	      topX1 = pmlEZXTopYZ[0][j][k];
	      topX2 = pmlEZYTopYZ[0][j][k];
	    }

	  commonFactor = mudt + 0.5 * pmlSigmaZStarBottomXY[i][j][k];
	  pmlHYZBottomXY[i][j][k] = pmlHYZBottomXY[i][j][k] *
	    (mudt - 0.5 * pmlSigmaZStarBottomXY[i][j][k]) / commonFactor -
	    (topZ1 - pmlEXYBottomXY[i][j][k] +
	     topZ2 - pmlEXZBottomXY[i][j][k]) /
	    (dz * pmlSZBottomXY[i][j][k] * commonFactor);
	  
	  commonFactor = mudt + 0.5 * pmlSigmaXStarBottomXY[i][j][k];
	  pmlHYXBottomXY[i][j][k] = pmlHYXBottomXY[i][j][k] *
	     (mudt - 0.5 * pmlSigmaXStarBottomXY[i][j][k])  / commonFactor +
	    (topX1 - pmlEZXBottomXY[i][j][k] +
	     topX2 - pmlEZYBottomXY[i][j][k]) /
	    (dx * pmlSXBottomXY[i][j][k] * commonFactor);
	  
	}

  /* topXY - plane */

  for(i=0; i<nx; i++)
    for(j=0; j<ny_plus_2pmlWidth; j++)
      for(k=0; k<pmlWidth; k++)
	{
	  if (i <  nx_minus_1)
	    {
	      topX1 = pmlEZXTopXY[i+1][j][k];
	      topX2 = pmlEZYTopXY[i+1][j][k];
	    }
	  else
	    {
	      topX1 = pmlEZXTopYZ[0][j][k+nz_plus_pmlWidth];
	      topX2 = pmlEZYTopYZ[0][j][k+nz_plus_pmlWidth];
	    }
	
	  commonFactor = mudt + 0.5 * pmlSigmaZStarTopXY[i][j][k];
	  pmlHYZTopXY[i][j][k] = pmlHYZTopXY[i][j][k] *
	    (mudt - 0.5 * pmlSigmaZStarTopXY[i][j][k]) / commonFactor -
	    (pmlEXYTopXY[i][j][k+1] - pmlEXYTopXY[i][j][k] +
	     pmlEXZTopXY[i][j][k+1] - pmlEXZTopXY[i][j][k]) /
	    (dz * pmlSZTopXY[i][j][k] * commonFactor);
	  
	  commonFactor = mudt + 0.5 * pmlSigmaXStarTopXY[i][j][k];
	  pmlHYXTopXY[i][j][k] = pmlHYXTopXY[i][j][k] *
	     (mudt - 0.5 * pmlSigmaXStarTopXY[i][j][k])  / commonFactor +
	    (topX1 - pmlEZXTopXY[i][j][k] +
	     topX2 - pmlEZYTopXY[i][j][k]) /
	    (dx * pmlSXTopXY[i][j][k] * commonFactor);
	}

  /* finally the XZ - plane (bottom) */
  
  for(i=0; i<nx; i++)
    for(j=0; j<pmlWidth; j++)
      for(k=0; k<nz; k++)
	{
	  if (i <  nx_minus_1)
	    {
	      topX1 = pmlEZXBottomXZ[i+1][j][k];
	      topX2 = pmlEZYBottomXZ[i+1][j][k];
	    }
	  else 
	    {
	      topX1 = pmlEZXTopYZ[0][j][k+pmlWidth];
	      topX2 = pmlEZYTopYZ[0][j][k+pmlWidth];
	    }
	  
	  if (k < nz_minus_1)
	    {
	      topZ1 = pmlEXYBottomXZ[i][j][k+1];
	      topZ2 = pmlEXZBottomXZ[i][j][k+1];
	    }
	  else
	    {
	      topZ1 = pmlEXYTopXY[i][j][0];
	      topZ2 = pmlEXZTopXY[i][j][0];
	    }

	  commonFactor = mudt + 0.5 * pmlSigmaZStarBottomXZ[i][j][k];
	  pmlHYZBottomXZ[i][j][k] = pmlHYZBottomXZ[i][j][k] *
	    (mudt - 0.5 * pmlSigmaZStarBottomXZ[i][j][k]) / commonFactor -
	    (topZ1 - pmlEXYBottomXZ[i][j][k] +
	     topZ2 - pmlEXZBottomXZ[i][j][k]) /
	    (dz * pmlSZBottomXZ[i][j][k] * commonFactor);
	  
	  commonFactor = mudt + 0.5 * pmlSigmaXStarBottomXZ[i][j][k];
	  pmlHYXBottomXZ[i][j][k] = pmlHYXBottomXZ[i][j][k] *
	     (mudt - 0.5 * pmlSigmaXStarBottomXZ[i][j][k])  / commonFactor +
	    (topX1 - pmlEZXBottomXZ[i][j][k] +
	     topX2 - pmlEZYBottomXZ[i][j][k]) /
	    (dx * pmlSXBottomXZ[i][j][k] * commonFactor);
	}

  /* and the top part */

  for(i=0; i<nx; i++)
    for(j=0; j<pmlWidth; j++)
      for(k=0; k<nz; k++)
	{
	  if (i <  nx_minus_1)
	    {
	      topX1 = pmlEZXTopXZ[i+1][j][k];
	      topX2 = pmlEZYTopXZ[i+1][j][k];
	    }
	  else
	    {
	      topX1 = pmlEZXTopYZ[0][j+ny_plus_pmlWidth][k+pmlWidth];
	      topX2 = pmlEZYTopYZ[0][j+ny_plus_pmlWidth][k+pmlWidth];
	    }
	  if (k < nz_minus_1)
	    {
	      topZ1 = pmlEXYTopXZ[i][j][k+1];
	      topZ2 = pmlEXZTopXZ[i][j][k+1];
	    }
	  else
	    {
	      topZ1 = pmlEXYTopXY[i][j+ny_plus_pmlWidth][0];
	      topZ2 = pmlEXZTopXY[i][j+ny_plus_pmlWidth][0];
	    }
	  
	  commonFactor = mudt + 0.5 * pmlSigmaZStarTopXZ[i][j][k];
	  pmlHYZTopXZ[i][j][k] = pmlHYZTopXZ[i][j][k] *
	    (mudt - 0.5 * pmlSigmaZStarTopXZ[i][j][k]) / commonFactor -
	    (topZ1 - pmlEXYTopXZ[i][j][k] +
	     topZ2 - pmlEXZTopXZ[i][j][k]) /
	    (dz * pmlSZTopXZ[i][j][k] * commonFactor);
	  
	  commonFactor = mudt + 0.5 * pmlSigmaXStarTopXZ[i][j][k];
	  pmlHYXTopXZ[i][j][k] = pmlHYXTopXZ[i][j][k] *
	     (mudt - 0.5 * pmlSigmaXStarTopXZ[i][j][k])  / commonFactor +
	    (topX1 - pmlEZXTopXZ[i][j][k] +
	     topX2 - pmlEZYTopXZ[i][j][k]) / /* !!!!! */
	    (dx * pmlSXTopXZ[i][j][k] * commonFactor);

	}
}

void UpdatePMLHz(void)
{
  int i,j,k;
  PRECISION mudt = MU_0 / dt;
  PRECISION topX1, topX2, topY1, topY2; /* 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 nz_plus_pmlWidth  = nz + pmlWidth;

  int pmlWidth_minus_1  = pmlWidth-1;

  int nx_minus_1        = nx-1;
  
  /* starting with YZ - plane (bottom) */
  
  for(i=0; i<pmlWidth; i++)
    for(j=0; j<ny_plus_2pmlWidth; j++)
      for(k=0; k<nz_plus_2pmlWidth; k++)
	{
	
	  if (i < pmlWidth_minus_1)
	    {
	      topX1 = pmlEYZBottomYZ[i+1][j][k];
	      topX2 = pmlEYXBottomYZ[i+1][j][k];
	    }
	  else
	    {
	      if (k < pmlWidth)
		{
		  topX1 = pmlEYZBottomXY[0][j][k];
		  topX2 = pmlEYXBottomXY[0][j][k];
		}
	      else if (k < nz_plus_pmlWidth)
		{
		  if (j < pmlWidth)
		    {
		      topX1 = pmlEYZBottomXZ[0][j][k-pmlWidth];
		      topX2 = pmlEYXBottomXZ[0][j][k-pmlWidth];
		    }
		  else if (j < ny_plus_pmlWidth)
		    {
		      topX1 = ey[0][j-pmlWidth][k-pmlWidth];
		      topX2 = 0.0;
		    }
		  else
		    {
		      topX1 = pmlEYZTopXZ[0][j-ny_plus_pmlWidth][k-pmlWidth];
		      topX2 = pmlEYXTopXZ[0][j-ny_plus_pmlWidth][k-pmlWidth];
		    }
		}
	      else
		{