maxwell_eps.c

MIT开发出来的计算光子晶体的软件

	  
	  /* first, the vertices of an icosohedron: */
	  x0 = sqrt(0.5 - sqrt(0.05));
	  y0 = sqrt(0.5 + sqrt(0.05));
	  z0 = 0;
	  if (num_sq_pts == 72)
	       w = 125 / 10080.0;
	  else
	       w = 1 / 12.0;
	  for (i = 0; i < 2; ++i) {
	       x0 = -x0;
	       for (j = 0; j < 2; ++j) {
		    y0 = -y0;
		    for (k = 0; k < 3; ++k) {
			 SHIFT3(x0,y0,z0);
			 x[n] = x0; y[n] = y0; z[n] = z0; weight[n++] = w;
		    }
	       }
	  }
	  
	  if (num_sq_pts == 72) {
	       /* it would be nice, for completeness, to have more
                  digits here: */
	       real coords[3][5] = {
		    { -0.151108275, 0.315838353, 0.346307112, -0.101808787, -0.409228403 },
		    { 0.155240600, 0.257049387, 0.666277790,  0.817386065, 0.501547712 },
		    { 0.976251323, 0.913330032, 0.660412970,  0.567022920, 0.762221757 }
	       };
	       
	       w = 143 / 10080.0;
	       for (l = 0; l < 5; ++l) {
		    x0 = coords[0][l]; y0 = coords[1][l]; z0 = coords[2][l];
		    for (i = 0; i < 3; ++i) {
			 double dummy = x0;
			 x0 = z0;
			 z0 = -y0;
			 y0 = -dummy;
			 for (j = 0; j < 3; ++j) {
			      SHIFT3(x0,y0,z0);
			      x[n] = x0; y[n] = y0; z[n] = z0; weight[n++] = w;
			 }
			 y0 = -y0;
			 z0 = -z0;
			 x[n] = x0; y[n] = y0; z[n] = z0; weight[n++] = w;
		    }
	       }
	  }
     }
     else
	  mpi_die("spherical_quadrature_points: passed unknown # points!");

     CHECK(n == num_sq_pts,
	   "bug in spherical_quadrature_points: wrong number of points!");
}

#define K_PI 3.141592653589793238462643383279502884197
#define SMALL 1.0e-6
#define MAX2(a,b) ((a) > (b) ? (a) : (b))
#define MIN2(a,b) ((a) < (b) ? (a) : (b))

#define NUM_SQ_PTS 50  /* number of spherical quadrature points to use */
#define MAX_MOMENT_MESH NUM_SQ_PTS /* max # of moment-mesh vectors */
#define MOMENT_MESH_R 0.5

/* A function to set up the mesh given the grid dimensions, mesh size,
   and lattice/reciprocal vectors.  (Any mesh sizes < 1 are taken to
   be 1.)  The returned values are:

   mesh_center: the center of the mesh, relative to the integer
                mesh coordinates; e.g. the mesh_center for a 3x3x3
                mesh is the point (1,1,1). 
   mesh_prod: the product of the mesh sizes.
   moment_mesh: an array of size_moment_mesh vectors, in lattice
                coordinates, of a "spherically-symmetric" mesh of
                points centered on the origin, designed to be
		used for averaging the first moment of epsilon at 
		a grid point (for finding the local surface normal).
   moment_mesh_weights: an array of size_moment_mesh weights to multiply
                        the integrand values by.  */
static void get_mesh(int nx, int ny, int nz, const int mesh_size[3],
		     real R[3][3], real G[3][3],
		     real mesh_center[3], int *mesh_prod,
		     real moment_mesh[MAX_MOMENT_MESH][3],
		     real moment_mesh_weights[MAX_MOMENT_MESH],
		     int *size_moment_mesh)
{
     int i,j;
     real min_diam = 1e20;
     real mesh_total[3] = { 0, 0, 0 };
     int rank = nz > 1 ? 3 : (ny > 1 ? 2 : 1);
     real weight_sum = 0.0;
	  
     *mesh_prod = 1;
     for (i = 0; i < 3; ++i) {
	  int ms = MAX2(mesh_size[i], 1);
	  mesh_center[i] = (ms - 1) * 0.5;
	  *mesh_prod *= ms;
     }

     /* Set the mesh to compute the "dipole moment" of epsilon over,
	in Cartesian coordinates (used to compute the normal vector
	at surfaces).  Ideally, a spherically-symmetrical distribution
	of points on the radius MOMENT_MESH_R sphere. */
     if (rank == 1) {
	  /* one-dimensional: just two points (really, normal
	     vector is always along x): */
	  *size_moment_mesh = 2;
	  moment_mesh[0][0] = -MOMENT_MESH_R;
	  moment_mesh[0][1] = 0.0;
	  moment_mesh[0][2] = 0.0;
	  moment_mesh[1][0] = MOMENT_MESH_R;
	  moment_mesh[1][1] = 0.0;
	  moment_mesh[1][2] = 0.0;
	  moment_mesh_weights[0] = moment_mesh_weights[1] = 0.5;
     }
     else if (rank == 2) {
	  /* two-dimensional: 12 points at 30-degree intervals: */
	  *size_moment_mesh = 12;
	  for (i = 0; i < *size_moment_mesh; ++i) {
	       moment_mesh[i][0] =
		    cos(2*i * K_PI / *size_moment_mesh) * MOMENT_MESH_R;
	       moment_mesh[i][1] =
		    sin(2*i * K_PI / *size_moment_mesh) * MOMENT_MESH_R;
	       moment_mesh[i][2] = 0.0;
	       moment_mesh_weights[i] = 1.0 / *size_moment_mesh;
	  }
     }
     else {
	  real x[NUM_SQ_PTS], y[NUM_SQ_PTS], z[NUM_SQ_PTS];

	  *size_moment_mesh = NUM_SQ_PTS;
	  spherical_quadrature_points(x,y,z, moment_mesh_weights, NUM_SQ_PTS);
	  for (i = 0; i < *size_moment_mesh; ++i) {
	       moment_mesh[i][0] = x[i] * MOMENT_MESH_R;
	       moment_mesh[i][1] = y[i] * MOMENT_MESH_R;
	       moment_mesh[i][2] = z[i] * MOMENT_MESH_R;
	  }
     }

     CHECK(*size_moment_mesh <= MAX_MOMENT_MESH, "yikes, moment mesh too big");

     for (i = 0; i < *size_moment_mesh; ++i)
	  weight_sum += moment_mesh_weights[i];
     CHECK(fabs(weight_sum - 1.0) < SMALL, "bug, incorrect moment weights");

     /* scale the moment-mesh vectors so that the sphere has a
	diameter of 2*MOMENT_MESH_R times the diameter of the
	smallest grid direction: */

     /* first, find the minimum distance between grid points along the
	lattice directions (should we use the maximum instead?): */
     for (i = 0; i < rank; ++i) {
	  real ri = R[i][0] * R[i][0] + R[i][1] * R[i][1] + R[i][2] * R[i][2];
	  ri = sqrt(ri) / (i == 0 ? nx : (i == 1 ? ny : nz));
	  min_diam = MIN2(min_diam, ri);
     }
     
     /* scale moment_mesh by this diameter: */
     for (i = 0; i < *size_moment_mesh; ++i) {
	  real len = 0;
	  for (j = 0; j < 3; ++j) {
	       moment_mesh[i][j] *= min_diam;
	       len += moment_mesh[i][j] * moment_mesh[i][j];
	       mesh_total[j] += moment_mesh[i][j];
	  }
	  CHECK(fabs(len - min_diam*min_diam*(MOMENT_MESH_R*MOMENT_MESH_R)) 
		< SMALL,
		"bug in get_mesh: moment_mesh vector is wrong length");
     }
     CHECK(fabs(mesh_total[0]) + fabs(mesh_total[1]) + fabs(mesh_total[2])
	   < SMALL, "bug in get_mesh: moment_mesh does not average to zero");

     /* Now, convert the moment_mesh vectors to lattice/grid coordinates;
        to do this, we multiply by the G matrix (inverse of R transposed) */
     for (i = 0; i < *size_moment_mesh; ++i) {
	  real v[3];
	  for (j = 0; j < 3; ++j)    /* make a copy of moment_mesh[i] */
	       v[j] = moment_mesh[i][j];
	  for (j = 0; j < 3; ++j)
	       moment_mesh[i][j] = G[j][0]*v[0] + G[j][1]*v[1] + G[j][2]*v[2];
     }
}

/**************************************************************************/

/* The following function initializes the dielectric tensor md->eps_inv,
   using the dielectric function epsilon(&eps, &eps_inv, r, epsilon_data).

   epsilon is averaged over a rectangular mesh spanning the space between
   grid points; the size of the mesh is given by mesh_size.

   R[0..2] are the spatial lattice vectors, and are used to convert
   the discretization grid into spatial coordinates (with the origin
   at the (0,0,0) grid element.

   In most places, the dielectric tensor is equal to eps_inv, but at
   dielectric interfaces it varies depending upon the polarization of
   the field (for faster convergence).  In particular, it depends upon
   the direction of the field relative to the surface normal vector,
   so we must compute the latter.  The surface normal is approximated
   by the "dipole moment" of the dielectric function over a spherical
   mesh.

   Implementation note: md->eps_inv is chosen to have dimensions matching
   the output of the FFT.  Thus, its dimensions depend upon whether we are
   doing a real or complex and serial or parallel FFT. */

void set_maxwell_dielectric(maxwell_data *md,
			    const int mesh_size[3],
			    real R[3][3], real G[3][3],
			    maxwell_dielectric_function epsilon,
			    void *epsilon_data)
{
     real s1, s2, s3, m1, m2, m3;  /* grid/mesh steps */
     real mesh_center[3];
     real moment_mesh[MAX_MOMENT_MESH][3];
     real moment_mesh_weights[MAX_MOMENT_MESH];
     real eps_inv_total = 0.0;
     int i, j, k;
     int mesh_prod;
     real mesh_prod_inv;
     int size_moment_mesh = 0;
     int n1, n2, n3;
#ifdef HAVE_MPI
     int local_n2, local_y_start, local_n3;
#endif
#ifndef SCALAR_COMPLEX
     int n_other, n_last, rank;
#endif

     n1 = md->nx; n2 = md->ny; n3 = md->nz;

     get_mesh(n1, n2, n3, mesh_size, R, G, 
	      mesh_center, &mesh_prod, moment_mesh, moment_mesh_weights,
	      &size_moment_mesh);
     mesh_prod_inv = 1.0 / mesh_prod;

#ifdef DEBUG
     mpi_one_printf("Using moment mesh (%d):\n", size_moment_mesh);
     for (i = 0; i < size_moment_mesh; ++i)
	  mpi_one_printf("   (%g, %g, %g) (%g)\n",
		 moment_mesh[i][0], moment_mesh[i][1], moment_mesh[i][2],
		 moment_mesh_weights[i]);
#endif

     s1 = 1.0 / n1;
     s2 = 1.0 / n2;
     s3 = 1.0 / n3;
     m1 = s1 / MAX2(1, mesh_size[0] - 1);
     m2 = s2 / MAX2(1, mesh_size[1] - 1);
     m3 = s3 / MAX2(1, mesh_size[2] - 1);

     /* Here we have different loops over the coordinates, depending
	upon whether we are using complex or real and serial or
        parallel transforms.  Each loop must define, in its body,
        variables (i2,j2,k2) describing the coordinate of the current
        point, and eps_index describing the corresponding index in 
	the array md->eps_inv[]. */

#ifdef SCALAR_COMPLEX

#  ifndef HAVE_MPI
     
     for (i = 0; i < n1; ++i)
	  for (j = 0; j < n2; ++j)
	       for (k = 0; k < n3; ++k)
     {
#         define i2 i
#         define j2 j
#         define k2 k
	  int eps_index = ((i * n2 + j) * n3 + k);

#  else /* HAVE_MPI */

     local_n2 = md->local_ny;
     local_y_start = md->local_y_start;

     /* first two dimensions are transposed in MPI output: */
     for (j = 0; j < local_n2; ++j)
          for (i = 0; i < n1; ++i)
	       for (k = 0; k < n3; ++k)
     {
#         define i2 i
	  int j2 = j + local_y_start;
#         define k2 k
	  int eps_index = ((j * n1 + i) * n3 + k);

#  endif /* HAVE_MPI */

#else /* not SCALAR_COMPLEX */

#  ifndef HAVE_MPI

     n_other = md->other_dims;