maxwell_op.c

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

   space and corresponds to the output from maxwell_compute_d_from_H,
   above. */
void maxwell_compute_e_from_d(maxwell_data *d,
			      scalar_complex *dfield,
			      int cur_num_bands)
{
     int i, b;

     CHECK(d, "null maxwell data pointer!");
     CHECK(dfield, "null field input/output data!");

     for (i = 0; i < d->fft_output_size; ++i) {
	  symmetric_matrix eps_inv = d->eps_inv[i];
	  for (b = 0; b < cur_num_bands; ++b) {
	       int ib = 3 * (i * cur_num_bands + b);
	       assign_symmatrix_vector(&dfield[ib], eps_inv, &dfield[ib]);
	  }
     }	  
}

/* Compute the magnetic (H) field in Fourier space from the electric
   field (e) in position space; this amouns to Fourier transforming and
   then taking the curl.  Also, multiply by scale.  Other
   parameters are as in compute_d_from_H. 

   Note: we actually compute (k+G) x E, whereas the actual H field
   is -i/omega i(k+G) x E...so, we are actually computing omega*H, here. */
void maxwell_compute_H_from_e(maxwell_data *d, evectmatrix Hout, 
			      scalar_complex *efield,
			      int cur_band_start, int cur_num_bands,
			      real scale)
{
     scalar *fft_data = (scalar *) efield;
     int i, j, b;

     CHECK(Hout.c == 2, "fields don't have 2 components!");
     CHECK(d, "null maxwell data pointer!");
     CHECK(efield, "null field output data!");
     CHECK(cur_band_start >= 0 && cur_band_start + cur_num_bands <= Hout.p,
	   "invalid range of bands for computing fields");

     /* convert back to Fourier space */
     maxwell_compute_fft(-1, d, fft_data,
			 cur_num_bands*3, cur_num_bands*3, 1);
     
     /* then, compute Hout = curl(fft_data) (* scale factor): */
     
     for (i = 0; i < d->other_dims; ++i)
	  for (j = 0; j < d->last_dim; ++j) {
	       int ij = i * d->last_dim + j;
	       int ij2 = i * d->last_dim_size + j;
	       k_data cur_k = d->k_plus_G[ij];
	       
	       for (b = 0; b < cur_num_bands; ++b)
		    assign_cross_c2t(&Hout.data[ij * 2 * Hout.p + 
					       b + cur_band_start],
				     Hout.p, cur_k, 
				     &fft_data[3 * (ij2*cur_num_bands + b)],
				     scale);
	  }
}


/* Compute H field in position space from Hin.  Parameters and output
   formats are the same as for compute_d_from_H, above. */
void maxwell_compute_h_from_H(maxwell_data *d, evectmatrix Hin, 
			      scalar_complex *hfield,
			      int cur_band_start, int cur_num_bands)
{
     scalar *fft_data = (scalar *) hfield;
     int i, j, b;

     CHECK(Hin.c == 2, "fields don't have 2 components!");
     CHECK(d, "null maxwell data pointer!");
     CHECK(hfield, "null field output data!");
     CHECK(cur_band_start >= 0 && cur_band_start + cur_num_bands <= Hin.p,
	   "invalid range of bands for computing fields");

     /* first, compute fft_data = Hin, with the vector field converted 
	from transverse to cartesian basis: */
     for (i = 0; i < d->other_dims; ++i)
	  for (j = 0; j < d->last_dim; ++j) {
	       int ij = i * d->last_dim + j;
	       int ij2 = i * d->last_dim_size + j;
               k_data cur_k = d->k_plus_G[ij];
	       
	       for (b = 0; b < cur_num_bands; ++b)
		    assign_t2c(&fft_data[3 * (ij2*cur_num_bands 
					      + b)], 
			       cur_k,
			       &Hin.data[ij * 2 * Hin.p + 
					b + cur_band_start],
			       Hin.p);
	  }

     /* now, convert to position space via FFT: */
     maxwell_compute_fft(+1, d, fft_data,
			 cur_num_bands*3, cur_num_bands*3, 1);
}

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

/* The following functions take a complex or real vector field
   in position space, as output by rfftwnd (or rfftwnd_mpi), and
   compute the "other half" of the array.

   That is, rfftwnd outputs only half of the logical FFT output,
   since the other half is redundant (see the FFTW manual).  This
   is fine for computation, but for visualization/output we
   want the whole array, redundant or not.  So, we output the array
   in two stages, first outputting the array we are given, then
   using the functions below to compute the other half and outputting
   that.

   Given an array A(i,j,k), the redundant half is given by the
   following identity for transforms of real data:

          A(nx-i,ny-j,nz-k) = A(i,j,k)*

   where nx-i/ny-j/nz-k are interpreted modulo nx/ny/nz.  (This
   means that zero coordinates are handled specially: nx-0 = 0.)

   Note that actually, the other "half" is actually slightly less
   than half of the array.  Note also that the other half, in the
   case of distributed MPI transforms, is not necessarily contiguous,
   due to special handling of zero coordinates.

   There is an additional complication.  The array with the symmetry
   above may have been multiplied by exp(ikx) phases to get its Bloch
   state.  In this case, we must use the identity:

          A(R-x)*exp(ik(R-x)) = ( A(x) exp(ikx) exp(-ikR) )*

   where R is a lattice vector.  That is, we not only must conjugate
   and reverse the order, but we also may need to multiply by exp(-ikR)
   before conjugating.  Unfortunately, R depends upon where we are
   in the array, because of the fact that the indices are interpreted
   modulo n (i.e. the zero indices aren't reordered).  e.g. for
   the point (nx-i,ny-j,nz-k), R = Rx*(i!=0) + Ry*(j!=0) + Rz*(k!=0).

   Another complication is that, for 2d rfftwnd_mpi transforms, the
   truncated dimension (in the transformed, transposed array) is
   the *first* dimension rather than the last.

   This code is a little too subtle for my tastes; real FFTs are a pain. */

#define TWOPI 6.2831853071795864769252867665590057683943388

/* This function takes a complex vector field and replaces it with its
   other "half."  phase{x,y,z} is the phase k*R{x,y,z}, in "units" of
   2*pi.  (Equivalently, phase{x,y,z} is the k vector in the
   reciprocal lattice basis.) */
void maxwell_vectorfield_otherhalf(maxwell_data *d, scalar_complex *field,
				   real phasex, real phasey, real phasez)
{
#ifndef SCALAR_COMPLEX
     int i, j, jmin = 1;
     int rank, n_other, n_last, n_last_stored, n_last_new, nx, ny, nz, nxmax;
#  ifdef HAVE_MPI
     int local_x_start;
#  endif
     scalar_complex pz, pxz, pyz, pxyz;

     nxmax = nx = d->nx; ny = d->ny; nz = d->nz;
     n_other = d->other_dims;
     n_last = d->last_dim;
     n_last_stored = d->last_dim_size / 2;
     n_last_new = n_last - n_last_stored; /* < n_last_stored always */
     rank = (nz == 1) ? (ny == 1 ? 1 : 2) : 3;

#  ifdef HAVE_MPI
     local_x_start = d->local_y_start;
     CHECK(rank == 2 || rank == 3, "unsupported rfftwnd_mpi rank!");
     if (rank == 2) {
	  n_other = nx;
	  n_last_new = ny = d->local_ny;
	  if (local_x_start == 0)
	       --n_last_new; /* DC frequency should not be in other half */
	  else
	       jmin = 0;
	  if (local_x_start + ny == n_last_stored && n_last % 2 == 0)
	       --n_last_new; /* Nyquist freq. should not be in other half */
	  n_last_stored = ny;
     }
     else { /* rank == 3 */
	  ny = nx;
	  nx = d->local_ny;
	  nxmax = local_x_start ? nx - 1 : nx;
	  n_other = nx * ny;
     }
#  endif /* HAVE_MPI */

     /* compute p = exp(i*phase) factors: */
     phasex *= -TWOPI; phasey *= -TWOPI; phasez *= -TWOPI;
     switch (rank) { /* treat z as the last/truncated dimension always */
	 case 3: break;
#  if defined(HAVE_MPI) && ! defined(SCALAR_COMPLEX)
	 case 2: phasez = phasex; phasex = phasey; phasey = 0; break;
#  else
	 case 2: phasez = phasey; phasey = 0; break;
#  endif
	 case 1: phasez = phasex; phasex = phasey = 0; break;
     }
     CASSIGN_SCALAR(pz, cos(phasez), sin(phasez));
     phasex += phasez;
     CASSIGN_SCALAR(pxz, cos(phasex), sin(phasex));
     phasex += phasey;
     CASSIGN_SCALAR(pxyz, cos(phasex), sin(phasex));
     phasey += phasez;
     CASSIGN_SCALAR(pyz, cos(phasey), sin(phasey));

/* convenience macros to copy vectors, vectors times phases, 
   and conjugated vectors: */
#  define ASSIGN_V(f,k,f2,k2) { f[3*(k)+0] = f2[3*(k2)+0];               \
                                f[3*(k)+1] = f2[3*(k2)+1];               \
                                f[3*(k)+2] = f2[3*(k2)+2]; }
#  define ASSIGN_VP(f,k,f2,k2,p) { CASSIGN_MULT(f[3*(k)+0], f2[3*(k2)+0], p); \
                                   CASSIGN_MULT(f[3*(k)+1], f2[3*(k2)+1], p); \
                                   CASSIGN_MULT(f[3*(k)+2], f2[3*(k2)+2], p); }
#  define ASSIGN_CV(f,k,f2,k2) { CASSIGN_CONJ(f[3*(k)+0], f2[3*(k2)+0]); \
                                 CASSIGN_CONJ(f[3*(k)+1], f2[3*(k2)+1]); \
                                 CASSIGN_CONJ(f[3*(k)+2], f2[3*(k2)+2]); }

     /* First, swap the order of elements and multiply by exp(ikR)
        phase factors.  We have to be careful here not to double-swap
        any element pair; this is prevented by never swapping with a
        "conjugated" point that is earlier in the array.  */

     if (rank == 3) {
	  int ix, iy;
	  for (ix = 0; 2*ix <= nxmax; ++ix) {
	       int xdiff, ixc;
#  ifdef HAVE_MPI
	       if (local_x_start == 0) {
		    xdiff = ix != 0; ixc = (nx - ix) % nx;
	       }
	       else {
		    xdiff = 1; ixc = nx-1 - ix;
	       }
#  else
	       xdiff = ix != 0; ixc = (nx - ix) % nx;
#  endif
	       for (iy = 0; iy < ny; ++iy) {
		    int ydiff = iy != 0;
		    int i = ix * ny + iy, ic = ixc * ny + (ny - iy) % ny, jmax;
		    if (ic < i)
			 continue;
		    jmax = n_last_new;
		    if (ic == i)
			 jmax = (jmax + 1) / 2;
		    for (j = 1; j <= jmax; ++j) {
			 int jc = n_last_new + 1 - j;
			 int ij = i*n_last_stored + j;
			 int ijc = ic*n_last_stored + jc;
			 scalar_complex f_tmp[3];
			 switch (xdiff*2 + ydiff) { /* pick exp(-ikR) phase */
			     case 3: /* xdiff && ydiff */
				  ASSIGN_VP(f_tmp, 0, field, ijc, pxyz);
				  ASSIGN_VP(field, ijc, field, ij, pxyz);
				  ASSIGN_V(field, ij, f_tmp, 0);
				  break;
			     case 2: /* xdiff && !ydiff */
				  ASSIGN_VP(f_tmp, 0, field, ijc, pxz);
				  ASSIGN_VP(field, ijc, field, ij, pxz);
				  ASSIGN_V(field, ij, f_tmp, 0);
				  break;
			     case 1: /* !xdiff && ydiff */
				  ASSIGN_VP(f_tmp, 0, field, ijc, pyz);
				  ASSIGN_VP(field, ijc, field, ij, pyz);
				  ASSIGN_V(field, ij, f_tmp, 0);
				  break;
			     case 0: /* !xdiff && !ydiff */
				  ASSIGN_VP(f_tmp, 0, field, ijc, pz);
				  ASSIGN_VP(field, ijc, field, ij, pz);
				  ASSIGN_V(field, ij, f_tmp, 0);
				  break;
			 }
		    }
	       }
	  }

	  /* Next, conjugate, and remove the holes from the array
	     corresponding to the DC and Nyquist frequencies (which were in
	     the first half already): */
	  for (i = 0; i < n_other; ++i)
	       for (j = 1; j < n_last_new + 1; ++j) {
		    int ij = i*n_last_stored + j, ijnew = i*n_last_new + j-1;
		    ASSIGN_CV(field, ijnew, field, ij);
	       }
     }
     else /* if (rank <= 2) */ {
	  int i;
	  if (rank == 1) /* (note that 1d MPI transforms are not allowed) */
	       nx = 1; /* x dimension is handled by j (last dimension) loop */

#  ifdef HAVE_MPI
	  for (i = 0; i < nx; ++i)
#  else
	  for (i = 0; 2*i <= nx; ++i)
#  endif
	  {
	       int xdiff = i != 0, ic = (nx - i) % nx;
	       int jmax = n_last_new + (jmin - 1);
#  ifndef HAVE_MPI
	       if (ic == i)
		    jmax = (jmax + 1) / 2;
#  endif
	       for (j = jmin; j <= jmax; ++j) {