maxwell_eps.c

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

     n_last = md->last_dim_size / 2;
     rank = (n3 == 1) ? (n2 == 1 ? 1 : 2) : 3;

     for (i = 0; i < n_other; ++i)
	  for (j = 0; j < n_last; ++j)
     {
	  int eps_index = i * n_last + j;
	  int i2, j2, k2;
	  switch (rank) {
	      case 2: i2 = i; j2 = j; k2 = 0; break;
	      case 3: i2 = i / n2; j2 = i % n2; k2 = j; break;
	      default: i2 = j; j2 = k2 = 0;  break;
	  }

#  else /* HAVE_MPI */

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

     /* For a real->complex transform, the last dimension is cut in
	half.  For a 2d transform, this is taken into account in local_ny
	already, but for a 3d transform we must compute the new n3: */
     if (n3 > 1)
	  local_n3 = md->last_dim_size / 2;
     else
	  local_n3 = 1;
     
     /* 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 < local_n3; ++k)
     {
#         define i2 i
	  int j2 = j + local_y_start;
#         define k2 k
	  int eps_index = ((j * n1 + i) * local_n3 + k);

#  endif  /* HAVE_MPI */

#endif /* not SCALAR_COMPLEX */

     {
	  int mi, mj, mk;
#ifdef WITH_HERMITIAN_EPSILON
	  symmetric_matrix eps_mean = {0,0,0,{0,0},{0,0},{0,0}}, 
			   eps_inv_mean = {0,0,0,{0,0},{0,0},{0,0}},
			   eps_mean_inv;
#else
	  symmetric_matrix eps_mean = {0,0,0,0,0,0}, 
			   eps_inv_mean = {0,0,0,0,0,0}, eps_mean_inv;
#endif
	  real norm_len;
	  real norm0, norm1, norm2;
	  short means_different_p, diag_eps_p;
	  
	  for (mi = 0; mi < mesh_size[0]; ++mi)
	       for (mj = 0; mj < mesh_size[1]; ++mj)
		    for (mk = 0; mk < mesh_size[2]; ++mk) {
			 real r[3];
			 symmetric_matrix eps, eps_inv;
			 r[0] = i2 * s1 + (mi - mesh_center[0]) * m1;
			 r[1] = j2 * s2 + (mj - mesh_center[1]) * m2;
			 r[2] = k2 * s3 + (mk - mesh_center[2]) * m3;
			 epsilon(&eps, &eps_inv, r, epsilon_data);
			 eps_mean.m00 += eps.m00;
			 eps_mean.m11 += eps.m11;
			 eps_mean.m22 += eps.m22;
			 eps_inv_mean.m00 += eps_inv.m00;
			 eps_inv_mean.m11 += eps_inv.m11;
			 eps_inv_mean.m22 += eps_inv.m22;
#ifdef WITH_HERMITIAN_EPSILON
			 CACCUMULATE_SUM(eps_mean.m01, eps.m01);
			 CACCUMULATE_SUM(eps_mean.m02, eps.m02);
			 CACCUMULATE_SUM(eps_mean.m12, eps.m12);
			 CACCUMULATE_SUM(eps_inv_mean.m01, eps_inv.m01);
			 CACCUMULATE_SUM(eps_inv_mean.m02, eps_inv.m02);
			 CACCUMULATE_SUM(eps_inv_mean.m12, eps_inv.m12);
#else
			 eps_mean.m01 += eps.m01;
			 eps_mean.m02 += eps.m02;
			 eps_mean.m12 += eps.m12;
			 eps_inv_mean.m01 += eps_inv.m01;
			 eps_inv_mean.m02 += eps_inv.m02;
			 eps_inv_mean.m12 += eps_inv.m12;
#endif
		    }
     
	  diag_eps_p = DIAG_SYMMETRIC_MATRIX(eps_mean);
	  if (diag_eps_p) { /* handle the common case of diagonal matrices: */
	       eps_mean_inv.m00 = mesh_prod / eps_mean.m00;
	       eps_mean_inv.m11 = mesh_prod / eps_mean.m11;
	       eps_mean_inv.m22 = mesh_prod / eps_mean.m22;
#ifdef WITH_HERMITIAN_EPSILON
	       CASSIGN_ZERO(eps_mean_inv.m01);
	       CASSIGN_ZERO(eps_mean_inv.m02);
	       CASSIGN_ZERO(eps_mean_inv.m12);
#else
	       eps_mean_inv.m01 = eps_mean_inv.m02 = eps_mean_inv.m12 = 0.0;
#endif
	       eps_inv_mean.m00 *= mesh_prod_inv;
	       eps_inv_mean.m11 *= mesh_prod_inv;
	       eps_inv_mean.m22 *= mesh_prod_inv;

	       means_different_p = 
		    fabs(eps_mean_inv.m00 - eps_inv_mean.m00) > SMALL ||
		    fabs(eps_mean_inv.m11 - eps_inv_mean.m11) > SMALL ||
		    fabs(eps_mean_inv.m22 - eps_inv_mean.m22) > SMALL;
	  }
	  else {
	       eps_inv_mean.m00 *= mesh_prod_inv;
	       eps_inv_mean.m11 *= mesh_prod_inv;
	       eps_inv_mean.m22 *= mesh_prod_inv;
	       eps_mean.m00 *= mesh_prod_inv;
	       eps_mean.m11 *= mesh_prod_inv;
	       eps_mean.m22 *= mesh_prod_inv;
#ifdef WITH_HERMITIAN_EPSILON
	       eps_mean.m01.re *= mesh_prod_inv;
	       eps_mean.m01.im *= mesh_prod_inv;
	       eps_mean.m02.re *= mesh_prod_inv;
	       eps_mean.m02.im *= mesh_prod_inv;
	       eps_mean.m12.re *= mesh_prod_inv;
	       eps_mean.m12.im *= mesh_prod_inv;
	       eps_inv_mean.m01.re *= mesh_prod_inv;
	       eps_inv_mean.m01.im *= mesh_prod_inv;
	       eps_inv_mean.m02.re *= mesh_prod_inv;
	       eps_inv_mean.m02.im *= mesh_prod_inv;
	       eps_inv_mean.m12.re *= mesh_prod_inv;
	       eps_inv_mean.m12.im *= mesh_prod_inv;
#else
	       eps_mean.m01 *= mesh_prod_inv;
	       eps_mean.m02 *= mesh_prod_inv;
	       eps_mean.m12 *= mesh_prod_inv;
	       eps_inv_mean.m01 *= mesh_prod_inv;
	       eps_inv_mean.m02 *= mesh_prod_inv;
	       eps_inv_mean.m12 *= mesh_prod_inv;
#endif
	       maxwell_sym_matrix_invert(&eps_mean_inv, &eps_mean);
	       
	       means_different_p = 
		    fabs(eps_mean_inv.m00 - eps_inv_mean.m00) > SMALL ||
		    fabs(eps_mean_inv.m11 - eps_inv_mean.m11) > SMALL ||
		    fabs(eps_mean_inv.m22 - eps_inv_mean.m22) > SMALL;
#ifdef WITH_HERMITIAN_EPSILON
	       means_different_p = means_different_p ||
		    fabs(eps_mean_inv.m01.re - eps_inv_mean.m01.re) > SMALL ||
		    fabs(eps_mean_inv.m02.re - eps_inv_mean.m02.re) > SMALL ||
		    fabs(eps_mean_inv.m12.re - eps_inv_mean.m12.re) > SMALL ||
		    fabs(eps_mean_inv.m01.im - eps_inv_mean.m01.im) > SMALL ||
		    fabs(eps_mean_inv.m02.im - eps_inv_mean.m02.im) > SMALL ||
		    fabs(eps_mean_inv.m12.im - eps_inv_mean.m12.im) > SMALL;
#else
	       means_different_p = means_different_p ||
		    fabs(eps_mean_inv.m01 - eps_inv_mean.m01) > SMALL ||
		    fabs(eps_mean_inv.m02 - eps_inv_mean.m02) > SMALL ||
		    fabs(eps_mean_inv.m12 - eps_inv_mean.m12) > SMALL;
#endif
	  }

	  /* if the two averaging methods yielded different results,
	     which usually happens if epsilon is not constant, then
	     we need to find the normal vector to the dielectric interface: */
	  if (means_different_p) {
	       real moment0 = 0, moment1 = 0, moment2 = 0;

	       for (mi = 0; mi < size_moment_mesh; ++mi) {
		    real r[3], eps_trace;
		    symmetric_matrix eps, eps_inv;
		    r[0] = i2 * s1 + moment_mesh[mi][0];
		    r[1] = j2 * s2 + moment_mesh[mi][1];
		    r[2] = k2 * s3 + moment_mesh[mi][2];
		    epsilon(&eps, &eps_inv, r, epsilon_data);
		    eps_trace = eps.m00 + eps.m11 + eps.m22;
		    eps_trace *= moment_mesh_weights[mi];
		    moment0 += eps_trace * moment_mesh[mi][0];
		    moment1 += eps_trace * moment_mesh[mi][1];
		    moment2 += eps_trace * moment_mesh[mi][2];
	       }
	       /* need to convert moment from lattice to cartesian coords: */
	       norm0 = R[0][0]*moment0 + R[1][0]*moment1 + R[2][0]*moment2;
	       norm1 = R[0][1]*moment0 + R[1][1]*moment1 + R[2][1]*moment2;
	       norm2 = R[0][2]*moment0 + R[1][2]*moment1 + R[2][2]*moment2;
	  
	       norm_len = sqrt(norm0*norm0 + norm1*norm1 + norm2*norm2);
	  }
	  
	  if (means_different_p && norm_len > SMALL) {
	       real x0, x1, x2;

	       norm_len = 1.0/norm_len;
	       norm0 *= norm_len;
	       norm1 *= norm_len;
	       norm2 *= norm_len;

	       /* Compute the effective inverse dielectric tensor.
		  We define this as:
                     1/2 ( {eps_inv_mean, P} + {eps_mean_inv, 1-P} )
	          where P is the projection matrix onto the normal direction
                  (P = norm ^ norm), and {a,b} is the anti-commutator ab+ba.
		   = 1/2 {eps_inv_mean - eps_mean_inv, P} + eps_mean_inv
		   = 1/2 (n_i conj(x_j) + x_i n_j) + (eps_mean_inv)_ij
		  where n_k is the kth component of the normal vector and 
		     x_i = (eps_inv_mean - eps_mean_inv)_ik n_k  
		  Note the implied summations (Einstein notation).

		  Note that the resulting matrix is symmetric, and we get just
		  eps_inv_mean if eps_inv_mean == eps_mean_inv, as desired.

	          Note that P is idempotent, so for scalar epsilon this
	          is just eps_inv_mean * P + eps_mean_inv * (1-P)
                        = (1/eps_inv_mean * P + eps_mean * (1-P)) ^ (-1),
	          which corresponds to the expression in the Meade paper. */
	       
	       x0 = (eps_inv_mean.m00 - eps_mean_inv.m00) * norm0;
	       x1 = (eps_inv_mean.m11 - eps_mean_inv.m11) * norm1;
	       x2 = (eps_inv_mean.m22 - eps_mean_inv.m22) * norm2;
	       if (diag_eps_p) {
#ifdef WITH_HERMITIAN_EPSILON
		    md->eps_inv[eps_index].m01.re = 0.5*(x0*norm1 + x1*norm0);
		    md->eps_inv[eps_index].m01.im = 0.0;
		    md->eps_inv[eps_index].m02.re = 0.5*(x0*norm2 + x2*norm0);
		    md->eps_inv[eps_index].m02.im = 0.0;
		    md->eps_inv[eps_index].m12.re = 0.5*(x1*norm2 + x2*norm1);
		    md->eps_inv[eps_index].m12.im = 0.0;
#else
		    md->eps_inv[eps_index].m01 = 0.5*(x0*norm1 + x1*norm0);
		    md->eps_inv[eps_index].m02 = 0.5*(x0*norm2 + x2*norm0);
		    md->eps_inv[eps_index].m12 = 0.5*(x1*norm2 + x2*norm1);
#endif
	       }
	       else {
#ifdef WITH_HERMITIAN_EPSILON
		    real x0i, x1i, x2i;
		    x0 += ((eps_inv_mean.m01.re - eps_mean_inv.m01.re)*norm1 + 
			   (eps_inv_mean.m02.re - eps_mean_inv.m02.re)*norm2);
		    x1 += ((eps_inv_mean.m01.re - eps_mean_inv.m01.re)*norm0 + 
			   (eps_inv_mean.m12.re - eps_mean_inv.m12.re)*norm2);
		    x2 += ((eps_inv_mean.m02.re - eps_mean_inv.m02.re)*norm0 +
			   (eps_inv_mean.m12.re - eps_mean_inv.m12.re)*norm1);
		    x0i = ((eps_inv_mean.m01.im - eps_mean_inv.m01.im)*norm1 + 
			   (eps_inv_mean.m02.im - eps_mean_inv.m02.im)*norm2);
		    x1i = (-(eps_inv_mean.m01.im - eps_mean_inv.m01.im)*norm0+ 
			   (eps_inv_mean.m12.im - eps_mean_inv.m12.im)*norm2);
		    x2i = -((eps_inv_mean.m02.im - eps_mean_inv.m02.im)*norm0 +
			    (eps_inv_mean.m12.im - eps_mean_inv.m12.im)*norm1);

		    md->eps_inv[eps_index].m01.re = (0.5*(x0*norm1 + x1*norm0) 
						     + eps_mean_inv.m01.re);
		    md->eps_inv[eps_index].m02.re = (0.5*(x0*norm2 + x2*norm0) 
						     + eps_mean_inv.m02.re);
		    md->eps_inv[eps_index].m12.re = (0.5*(x1*norm2 + x2*norm1) 
						     + eps_mean_inv.m12.re);
		    md->eps_inv[eps_index].m01.im = (0.5*(x0i*norm1-x1i*norm0) 
						     + eps_mean_inv.m01.im);
		    md->eps_inv[eps_index].m02.im = (0.5*(x0i*norm2-x2i*norm0) 
						     + eps_mean_inv.m02.im);
		    md->eps_inv[eps_index].m12.im = (0.5*(x1i*norm2-x2i*norm1) 
						     + eps_mean_inv.m12.im);
#else
		    x0 += ((eps_inv_mean.m01 - eps_mean_inv.m01) * norm1 + 
			   (eps_inv_mean.m02 - eps_mean_inv.m02) * norm2);
		    x1 += ((eps_inv_mean.m01 - eps_mean_inv.m01) * norm0 + 
			   (eps_inv_mean.m12 - eps_mean_inv.m12) * norm2);
		    x2 += ((eps_inv_mean.m02 - eps_mean_inv.m02) * norm0 +
			   (eps_inv_mean.m12 - eps_mean_inv.m12) * norm1);

		    md->eps_inv[eps_index].m01 = (0.5*(x0*norm1 + x1*norm0) 
						  + eps_mean_inv.m01);
		    md->eps_inv[eps_index].m02 = (0.5*(x0*norm2 + x2*norm0) 
						  + eps_mean_inv.m02);
		    md->eps_inv[eps_index].m12 = (0.5*(x1*norm2 + x2*norm1) 
						  + eps_mean_inv.m12);
#endif
	       }
	       md->eps_inv[eps_index].m00 = x0*norm0 + eps_mean_inv.m00;
	       md->eps_inv[eps_index].m11 = x1*norm1 + eps_mean_inv.m11;
	       md->eps_inv[eps_index].m22 = x2*norm2 + eps_mean_inv.m22;
	  }
	  else { /* undetermined normal vector and/or constant eps */
	       md->eps_inv[eps_index] = eps_mean_inv;
	  }
	  
	  eps_inv_total += (md->eps_inv[eps_index].m00 + 
			    md->eps_inv[eps_index].m11 + 
			    md->eps_inv[eps_index].m22);
     }}  /* end of loop body */

     mpi_allreduce_1(&eps_inv_total, real, SCALAR_MPI_TYPE,
		     MPI_SUM, MPI_COMM_WORLD);
     n1 = md->fft_output_size;
     mpi_allreduce_1(&n1, int, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
     md->eps_inv_mean = eps_inv_total / (3 * n1);
}