fields.c

麻省理工的计算光子晶体的程序

	  CASSIGN_MULT(field[2], field[2], phase);
     }

     F.x = cscalar2cnumber(field[0]);
     F.y = cscalar2cnumber(field[1]);
     F.z = cscalar2cnumber(field[2]);
     return F;
}

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

/* compute the fraction of the field energy that is located in the
   given range of dielectric constants: */
number compute_energy_in_dielectric(number eps_low, number eps_high)
{
     int N, i, last_dim, last_dim_stored, nx, nz, local_y_start;
     real *energy = (real *) curfield;
     real epsilon, energy_sum = 0.0;

     if (!curfield || !strchr("DHR", curfield_type)) {
          mpi_one_fprintf(stderr, "The D or H energy density must be loaded first.\n");
          return 0.0;
     }

     N = mdata->fft_output_size;
     last_dim = mdata->last_dim;
     last_dim_stored =
	  mdata->last_dim_size / (sizeof(scalar_complex)/sizeof(scalar));
     nx = mdata->nx; nz = mdata->nz; local_y_start = mdata->local_y_start;

     for (i = 0; i < N; ++i) {
	  epsilon = mean_epsilon(mdata->eps_inv +i);
	  if (epsilon >= eps_low && epsilon <= eps_high) {
	       energy_sum += energy[i];
#ifndef SCALAR_COMPLEX
	       /* most points are counted twice, by rfftw output symmetry: */
	       {
		    int last_index;
#  ifdef HAVE_MPI
		    if (nz == 1) /* 2d: 1st dim. is truncated one */
			 last_index = i / nx + local_y_start;
		    else
			 last_index = i % last_dim_stored;
#  else
		    last_index = i % last_dim_stored;
#  endif
		    if (last_index != 0 && 2*last_index != last_dim)
			 energy_sum += energy[i];
	       }
#endif
	  }
     }
     mpi_allreduce_1(&energy_sum, real, SCALAR_MPI_TYPE,
		     MPI_SUM, MPI_COMM_WORLD);
     energy_sum *= Vol / H.N;
     return energy_sum;
}

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

/* Prepend the prefix to the fname, and (if parity_suffix is true)
   append a parity specifier (if any) (e.g. ".te"), returning a new
   string, which should be deallocated with free().  fname or prefix
   may be NULL, in which case they are treated as the empty string. */
static char *fix_fname(const char *fname, const char *prefix,
		       maxwell_data *d, int parity_suffix)
{
     char *s;
     CHK_MALLOC(s, char,
		(fname ? strlen(fname) : 0) + 
		(prefix ? strlen(prefix) : 0) + 20);
     strcpy(s, prefix ? prefix : "");
     strcat(s, fname ? fname : "");
     if (parity_suffix && d->parity != NO_PARITY) {
	  /* assumes parity suffix is less than 20 characters;
	     currently it is less than 12 */
	  strcat(s, ".");
	  strcat(s, parity_string(d));
     }
     return s;
}

static void output_scalarfield(real *vals,
			       const int dims[3],
			       const int local_dims[3],
			       const int start[3],
			       matrixio_id file_id,
			       const char *dataname,
			       int last_dim_index,
			       int last_dim_start, int last_dim_size,
			       int first_dim_start, int first_dim_size,
			       int write_start0_special)
{
     matrixio_id data_id = -1;

     fieldio_write_real_vals(vals, 3, dims, 
			     local_dims, start, file_id, 0, 
			     dataname, &data_id);
     
#ifndef SCALAR_COMPLEX
     {
	  int start_new[3], local_dims_new[3];

	  start_new[0] = start[0];
	  start_new[1] = start[1];
	  start_new[2] = start[2];
	  local_dims_new[0] = local_dims[0];
	  local_dims_new[1] = local_dims[1];
	  local_dims_new[2] = local_dims[2];

	  maxwell_scalarfield_otherhalf(mdata, vals);
	  start_new[last_dim_index] = last_dim_start;
	  local_dims_new[last_dim_index] = last_dim_size;
	  start_new[0] = first_dim_start;
	  local_dims_new[0] = first_dim_size;
	  if (write_start0_special) {
	       /* The conjugated array half may be discontiguous.
		  First, write the part not containing start_new[0], and
		  then write the start_new[0] slab. */
	       fieldio_write_real_vals(vals +
				       local_dims_new[1] * local_dims_new[2],
				       3, dims, local_dims_new, start_new, 
				       file_id, 1, dataname, &data_id);
	       local_dims_new[0] = 1;
	       start_new[0] = 0;
	       fieldio_write_real_vals(vals, 3, dims, 
				       local_dims_new, start_new,
				       file_id, 1, dataname, &data_id);
	  }
	  else {
	       fieldio_write_real_vals(vals, 3, dims, 
				       local_dims_new, start_new, 
				       file_id, 1, dataname, &data_id);
	  }
     }
#endif

     if (data_id >= 0)
	  matrixio_close_dataset(data_id);
}

/* given the field in curfield, store it to HDF (or whatever) using
   the matrixio (fieldio) routines.  Allow the component to be specified
   (which_component 0/1/2 = x/y/z, -1 = all) for vector fields. 
   Also allow the user to specify a prefix string for the filename. */
void output_field_to_file(integer which_component, string filename_prefix)
{
     char fname[100], *fname2, description[100];
     int dims[3], local_dims[3], start[3] = {0,0,0};
     matrixio_id file_id = -1;
     int attr_dims[2] = {3, 3};
     real output_k[3]; /* kvector in reciprocal lattice basis */
     real output_R[3][3];

     /* where to put "otherhalf" block of output, only used for real scalars */
     int last_dim_index = 0;
     int last_dim_start = 0, last_dim_size = 0;
     int first_dim_start = 0, first_dim_size = 0;
     int write_start0_special = 0;

     if (!curfield) {
	  mpi_one_fprintf(stderr, 
		  "fields, energy dens., or epsilon must be loaded first.\n");
	  return;
     }
     
#ifdef HAVE_MPI
     /* The first two dimensions (x and y) of the position-space fields
	are transposed when we use MPI, so we need to transpose everything. */
     dims[0] = mdata->ny;
     local_dims[1] = dims[1] = mdata->nx;
     local_dims[2] = dims[2] = mdata->nz;
     local_dims[0] = mdata->local_ny;
     start[0] = mdata->local_y_start;
#  ifndef SCALAR_COMPLEX
     /* Ugh, hairy.  See also maxwell_vectorfield_otherhalf. */
     if (dims[2] == 1) {
	  last_dim_index = 0;
	  first_dim_size = local_dims[0];
	  first_dim_start = dims[0] - (start[0] + local_dims[0] - 1);

	  if (start[0] == 0)
	       --first_dim_size; /* DC frequency is not in other half */
	  if (start[0] + local_dims[0] == mdata->last_dim_size / 2 &&
	      dims[0] % 2 == 0) {
	       --first_dim_size; /* Nyquist frequency is not in other half */
	       ++first_dim_start;
	  }

	  last_dim_start = first_dim_start;
	  last_dim_size = first_dim_size;
     }
     else {
	  last_dim_index = 2;
	  local_dims[last_dim_index] = mdata->last_dim_size / 2;
	  if (start[0] == 0) {
	       first_dim_size = local_dims[0] - 1;
	       first_dim_start = dims[0] - first_dim_size;
	       write_start0_special = 1;
	  }
	  else {
	       first_dim_start = dims[0] - (start[0] + local_dims[0] - 1);
	       first_dim_size = local_dims[0];
	  }
	  last_dim_start = local_dims[last_dim_index];
	  last_dim_size = dims[last_dim_index] - local_dims[last_dim_index];
     }
#  endif /* ! SCALAR_COMPLEX */
     output_k[0] = R[1][0]*mdata->current_k[0] + R[1][1]*mdata->current_k[1]
	  + R[1][2]*mdata->current_k[2];
     output_k[1] = R[0][0]*mdata->current_k[0] + R[0][1]*mdata->current_k[1]
	  + R[0][2]*mdata->current_k[2];
     output_k[2] = R[2][0]*mdata->current_k[0] + R[2][1]*mdata->current_k[1]
	  + R[2][2]*mdata->current_k[2];
     output_R[0][0]=R[1][0]; output_R[0][1]=R[1][1]; output_R[0][2]=R[1][2];
     output_R[1][0]=R[0][0]; output_R[1][1]=R[0][1]; output_R[1][2]=R[0][2];
     output_R[2][0]=R[2][0]; output_R[2][1]=R[2][1]; output_R[2][2]=R[2][2];
#else /* ! HAVE_MPI */
     dims[0] = mdata->nx;
     local_dims[1] = dims[1] = mdata->ny;
     local_dims[2] = dims[2] = mdata->nz;
     local_dims[0] = mdata->local_nx;
#  ifndef SCALAR_COMPLEX
     last_dim_index = dims[2] == 1 ? (dims[1] == 1 ? 0 : 1) : 2;
     local_dims[last_dim_index] = mdata->last_dim_size / 2;
     last_dim_start = local_dims[last_dim_index];
     last_dim_size = dims[last_dim_index] - local_dims[last_dim_index];
     first_dim_start = last_dim_index ? 0 : last_dim_start;
     first_dim_size = last_dim_index ? local_dims[0] : last_dim_size;
#  endif
     start[0] = mdata->local_x_start;
     output_k[0] = R[0][0]*mdata->current_k[0] + R[0][1]*mdata->current_k[1]
	  + R[0][2]*mdata->current_k[2];
     output_k[1] = R[1][0]*mdata->current_k[0] + R[1][1]*mdata->current_k[1]
	  + R[1][2]*mdata->current_k[2];
     output_k[2] = R[2][0]*mdata->current_k[0] + R[2][1]*mdata->current_k[1]
	  + R[2][2]*mdata->current_k[2];
     output_R[0][0]=R[0][0]; output_R[0][1]=R[0][1]; output_R[0][2]=R[0][2];
     output_R[1][0]=R[1][0]; output_R[1][1]=R[1][1]; output_R[1][2]=R[1][2];
     output_R[2][0]=R[2][0]; output_R[2][1]=R[2][1]; output_R[2][2]=R[2][2];
#endif /* ! HAVE_MPI */

     if (strchr("Rv", curfield_type)) /* generic scalar/vector field */
	  output_k[0] = output_k[1] = output_k[2] = 0.0; /* don't know k */
     
     if (strchr("dhecv", curfield_type)) { /* outputting vector field */
	  matrixio_id data_id[6] = {-1,-1,-1,-1,-1,-1};
	  int i;

	  sprintf(fname, "%c.k%02d.b%02d",
		  curfield_type, kpoint_index, curfield_band);
	  if (which_component >= 0) {
	       char comp_str[] = ".x";
	       comp_str[1] = 'x' + which_component;
	       strcat(fname, comp_str);
	  }
	  sprintf(description, "%c field, kpoint %d, band %d, freq=%g",
		  curfield_type, kpoint_index, curfield_band, 
		  freqs.items[curfield_band - 1]);
	  fname2 = fix_fname(fname, filename_prefix, mdata, 1);
	  mpi_one_printf("Outputting fields to %s...\n", fname2);
	  file_id = matrixio_create(fname2);
	  free(fname2);
	  fieldio_write_complex_field(curfield, 3, dims, local_dims, start,
				      which_component, output_k,
				      file_id, 0, data_id);

#ifndef SCALAR_COMPLEX
	  /* Here's where it gets hairy. */
	  maxwell_vectorfield_otherhalf(mdata, curfield,
					output_k[0], output_k[1], output_k[2]);
	  start[last_dim_index] = last_dim_start;
	  local_dims[last_dim_index] = last_dim_size;
	  start[0] = first_dim_start;
	  local_dims[0] = first_dim_size;
	  if (write_start0_special) {
	       /* The conjugated array half may be discontiguous.
		  First, write the part not containing start[0], and
		  then write the start[0] slab. */
	       fieldio_write_complex_field(curfield + 
					   3 * local_dims[1] * local_dims[2],
					   3, dims, local_dims, start,
					   which_component, NULL, 
					   file_id, 1, data_id);
	       local_dims[0] = 1;
	       start[0] = 0;
	       fieldio_write_complex_field(curfield, 3, dims,local_dims,start,
					   which_component, NULL,
					   file_id, 1, data_id);
	  }
	  else {
	       fieldio_write_complex_field(curfield, 3, dims,local_dims,start,
					   which_component, NULL,
					   file_id, 1, data_id);
	  }
#endif

	  for (i = 0; i < 6; ++i)
	       if (data_id[i] >= 0)
		    matrixio_close_dataset(data_id[i]);
	  matrixio_write_data_attr(file_id, "Bloch wavevector",
				   output_k, 1, attr_dims);
     }
     else if (strchr("DHnR", curfield_type)) { /* scalar field */
	  if (curfield_type == 'n') {
	       sprintf(fname, "epsilon");
	       sprintf(description, "dielectric function, epsilon");
	  }
	  else {
	       sprintf(fname, "%cpwr.k%02d.b%02d",
		       tolower(curfield_type), kpoint_index, curfield_band);
	       sprintf(description,
		       "%c field energy density, kpoint %d, band %d, freq=%g",
		       curfield_type, kpoint_index, curfield_band, 
		       freqs.items[curfield_band - 1]);
	  }
	  fname2 = fix_fname(fname, filename_prefix, mdata, 
			     /* no parity suffix for epsilon: */
			     curfield_type != 'n');
	  mpi_one_printf("Outputting %s...\n", fname2);
	  file_id = matrixio_create(fname2);
	  free(fname2);

	  output_scalarfield((real *) curfield, dims, 
			     local_dims, start, file_id, "data",
			     last_dim_index, last_dim_start, last_dim_size,
			     first_dim_start, first_dim_size,
			     write_start0_special);

	  if (curfield_type == 'n') {
	       int c1, c2, inv;
	       char dataname[100];

	       for (inv = 0; inv < 2; ++inv)
		    for (c1 = 0; c1 < 3; ++c1)
			 for (c2 = c1; c2 < 3; ++c2) {
			      get_epsilon_tensor(c1,c2, 0, inv);
			      sprintf(dataname, "%s.%c%c",
				      inv ? "epsilon_inverse" : "epsilon",
				      c1 + 'x', c2 + 'x');
			      output_scalarfield((real *) curfield, dims,
						 local_dims, start,
						 file_id, dataname,
						 last_dim_index,
						 last_dim_start, last_dim_size,
						 first_dim_start,
						 first_dim_size,
						 write_start0_special);
#if defined(WITH_HERMITIAN_EPSILON)
			      if (c1 != c2) {
				   get_epsilon_tensor(c1,c2, 1, inv);
				   strcat(dataname, ".i");
#ifndef SCALAR_COMPLEX /* scalarfield_otherhalf isn't right */
				   strcat(dataname, ".screwy");
#endif
				   output_scalarfield((real *) curfield, dims,
						      local_dims, start,
						      file_id, dataname,
						      last_dim_index,
						      last_dim_start,
						      last_dim_size,
						      first_dim_start,
						      first_dim_size,
						      write_start0_special);
			      }
#endif
			 }
	  }

     }
     else
	  mpi_one_fprintf(stderr, "unknown field type!\n");

     if (file_id >= 0) {
	  matrixio_write_data_attr(file_id, "lattice vectors",
				   &output_R[0][0], 2, attr_dims);
	  matrixio_write_string_attr(file_id, "description", description);

	  matrixio_close(file_id);
     }

     /* We have destroyed curfield (by multiplying it by phases,
	and/or reorganizing in the case of real-amplitude fields). */
     curfield_reset();
}

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

/* For curfield an energy density, compute the fraction of the energy
   that resides inside the given list of geometric objects.   Later
   objects in the list have precedence, just like the ordinary
   geometry list. */
number compute_energy_in_object_list(geometric_object_list objects)