maxwell_op.c

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

		    scalar_complex f_tmp[3];
#  ifdef HAVE_MPI
		    int jc = jmax + jmin - j;
		    int ij = j * nx + i;
		    int ijc = jc * nx + ic;
		    if (ijc < ij)
			 break;
#  else  /* ! HAVE_MPI */
		    int jc = n_last_new + 1 - j;
		    int ij = i*n_last_stored + j;
		    int ijc = ic*n_last_stored + jc;
#  endif /* ! HAVE_MPI */
		    if (xdiff) {
			 ASSIGN_VP(f_tmp, 0, field, ijc, pxz);
			 ASSIGN_VP(field, ijc, field, ij, pxz);
			 ASSIGN_V(field, ij, f_tmp, 0);
		    }
		    else {
			 ASSIGN_VP(f_tmp, 0, field, ijc, pz);
			 ASSIGN_VP(field, ijc, field, ij, pz);
			 ASSIGN_V(field, ij, f_tmp, 0);
		    }
	       }
	  }

	  /* 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 < nx; ++i)
	       for (j = jmin; j < n_last_new + jmin; ++j) {
#  ifdef HAVE_MPI
		    int ij = j*nx + i, ijnew = (j-jmin)*nx + i;
#  else
		    int ij = i*n_last_stored + j, ijnew = i*n_last_new + j-1;
#  endif
		    ASSIGN_CV(field, ijnew, field, ij);
	       }
     }

#  undef ASSIGN_V
#  undef ASSIGN_VP
#  undef ASSIGN_CV

#endif /* ! SCALAR_COMPLEX */
}

/* Similar to vectorfield_otherhalf, above, except that it operates on
   a real scalar field, which is assumed to have come from e.g. the
   absolute values of a complex field (and thus no phase factors or
   conjugations are necessary). */
void maxwell_scalarfield_otherhalf(maxwell_data *d, real *field)
{
#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

     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 */

     /* 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 ixc;
#  ifdef HAVE_MPI
	       if (local_x_start == 0)
		    ixc = (nx - ix) % nx;
	       else
		    ixc = nx-1 - ix;
#  else
	       ixc = (nx - ix) % nx;
#  endif
	       for (iy = 0; iy < ny; ++iy) {
		    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;
			 real f_tmp;
			 f_tmp = field[ijc];
			 field[ijc] = field[ij];
			 field[ij] = f_tmp;
		    }
	       }
	  }

	  /* 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;
		    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 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) {
		    real f_tmp;
#  ifdef HAVE_MPI
		    int jc = jmax + jmin - j;
		    int ij = j * nx + i;
		    int ijc = jc * nx + ic;
		    if (ijc < ij)
			 break;
#  else  /* ! HAVE_MPI */
		    int jc = n_last_new + 1 - j;
		    int ij = i*n_last_stored + j;
		    int ijc = ic*n_last_stored + jc;
#  endif /* ! HAVE_MPI */
		    f_tmp = field[ijc];
		    field[ijc] = field[ij];
		    field[ij] = f_tmp;
	       }
	  }

	  /* Next, remove the holes from the array corresponding to
	     the DC and Nyquist frequencies (which were in the first
	     half already): */
	  for (i = 0; i < nx; ++i)
	       for (j = jmin; j < n_last_new + jmin; ++j) {
#  ifdef HAVE_MPI
		    int ij = j*nx + i, ijnew = (j-jmin)*nx + i;
#  else
		    int ij = i*n_last_stored + j, ijnew = i*n_last_new + j-1;
#  endif
		    field[ijnew] = field[ij];
	       }
     }
#endif /* ! SCALAR_COMPLEX */
}

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

#define MIN2(a,b) ((a) < (b) ? (a) : (b))

/* Compute Xout = curl(1/epsilon * curl(Xin)) */
void maxwell_operator(evectmatrix Xin, evectmatrix Xout, void *data,
		      int is_current_eigenvector, evectmatrix Work)
{
     maxwell_data *d = (maxwell_data *) data;
     int cur_band_start;
     scalar_complex *cdata;
     real scale;
     
     CHECK(d, "null maxwell data pointer!");
     CHECK(Xin.c == 2, "fields don't have 2 components!");

     (void) is_current_eigenvector;  /* unused */
     (void) Work;

     cdata = (scalar_complex *) d->fft_data;
     scale = -1.0 / Xout.N;  /* scale factor to normalize FFT; 
				negative sign comes from 2 i's from curls */

     /* compute the operator, num_fft_bands at a time: */
     for (cur_band_start = 0; cur_band_start < Xin.p; 
	  cur_band_start += d->num_fft_bands) {
	  int cur_num_bands = MIN2(d->num_fft_bands, Xin.p - cur_band_start);

	  maxwell_compute_d_from_H(d, Xin, cdata,
				   cur_band_start, cur_num_bands);
	  maxwell_compute_e_from_d(d, cdata, cur_num_bands);
	  maxwell_compute_H_from_e(d, Xout, cdata,
				   cur_band_start, cur_num_bands, scale);
     }
}

/* Compute the operation Xout = (M - w^2) Xin, where M is the Maxwell
   operator and w is the target frequency.  This shifts the eigenvalue
   spectrum so that the smallest magnitude eigenvalues are those
   nearest to w.  However, there are negative eigenvalues (the
   operator is not positive-definite), and the smallest eigenvectors
   (not taking the absolute value) are the same as those of M. */
void maxwell_target_operator1(evectmatrix Xin, evectmatrix Xout, void *data,
			     int is_current_eigenvector, evectmatrix Work)
{
     maxwell_target_data *d = (maxwell_target_data *) data;
     real omega_sqr = d->target_frequency * d->target_frequency;

     maxwell_operator(Xin, Xout, d->d, is_current_eigenvector, Work);
     evectmatrix_aXpbY(1.0, Xout, -omega_sqr, Xin);
}

/* Compute the operation Xout = (M - w^2)^2 Xin, where M is the
   Maxwell operator and w is the target frequency.  This shifts the
   eigenvalue spectrum so that the smallest eigenvalues are those
   nearest to w. */
void maxwell_target_operator(evectmatrix Xin, evectmatrix Xout, void *data,
			     int is_current_eigenvector, evectmatrix Work)
{
     if (Xin.n != 0)
	  CHECK(Work.data && Work.data != Xin.data && Work.data != Xout.data,
		"maxwell_target_operator must have distinct workspace!");

     maxwell_target_operator1(Xin, Work, data, is_current_eigenvector, Xout);

     /* N.B. maxwell_operator(), and thus maxwell_target_operator1(),
	doesn't actually need the workspace, so we can safely pass
	Work here for the scratch parameter: */
     maxwell_target_operator1(Work, Xout, data, is_current_eigenvector, Work);
}

/* Compute the operation Xout = curl 1/epsilon * i u x Xin, which 
   is useful operation in computing the group velocity (derivative
   of the maxwell operator).  u is a vector in cartesian coordinates. */
void maxwell_ucross_op(evectmatrix Xin, evectmatrix Xout,
		       maxwell_data *d, const real u[3])
{
     scalar *fft_data;
     scalar_complex *cdata;
     real scale;
     int cur_band_start;
     int i, j, b;

     CHECK(d, "null maxwell data pointer!");
     CHECK(Xin.c == 2, "fields don't have 2 components!");

     cdata = (scalar_complex *) (fft_data = d->fft_data);
     scale = -1.0 / Xout.N;  /* scale factor to normalize FFT;
                                negative sign comes from 2 i's from curls */

     /* compute the operator, num_fft_bands at a time: */
     for (cur_band_start = 0; cur_band_start < Xin.p;
          cur_band_start += d->num_fft_bands) {
          int cur_num_bands = MIN2(d->num_fft_bands, Xin.p - cur_band_start);
	  
	  /* first, compute fft_data = u x Xin: */
	  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_ucross_t2c(&fft_data[3 * (ij2*cur_num_bands
							  + b)], 
					   u, cur_k, 
					   &Xin.data[ij * 2 * Xin.p + 
						    b + cur_band_start],
					   Xin.p);
	       }
	  
	  /* now, convert to position space via FFT: */
	  maxwell_compute_fft(+1, d, fft_data,
			      cur_num_bands*3, cur_num_bands*3, 1);
	  
          maxwell_compute_e_from_d(d, cdata, cur_num_bands);
          maxwell_compute_H_from_e(d, Xout, cdata,
                                   cur_band_start, cur_num_bands, scale);
     }
}