fields.c

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

{
     int i, j, k, n1, n2, n3, n_other, n_last, rank, last_dim;
#ifdef HAVE_MPI
     int local_n2, local_y_start, local_n3;
#endif
     real s1, s2, s3, c1, c2, c3;
     real *energy = (real *) curfield;
     real energy_sum = 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;
     }

     for (i = 0; i < objects.num_items; ++i)
	  geom_fix_object(objects.items[i]);

     n1 = mdata->nx; n2 = mdata->ny; n3 = mdata->nz;
     n_other = mdata->other_dims;
     n_last = mdata->last_dim_size / (sizeof(scalar_complex)/sizeof(scalar));
     last_dim = mdata->last_dim;
     rank = (n3 == 1) ? (n2 == 1 ? 1 : 2) : 3;

     s1 = geometry_lattice.size.x / n1;
     s2 = geometry_lattice.size.y / n2;
     s3 = geometry_lattice.size.z / n3;
     c1 = n1 <= 1 ? 0 : geometry_lattice.size.x * 0.5;
     c2 = n2 <= 1 ? 0 : geometry_lattice.size.y * 0.5;
     c3 = n3 <= 1 ? 0 : geometry_lattice.size.z * 0.5;

     /* 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 "index" describing the corresponding index in 
	the curfield array.

        This was all stolen from maxwell_eps.c...it would be better
        if we didn't have to cut and paste, sigh. */

#ifdef SCALAR_COMPLEX

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

#  else /* HAVE_MPI */

     local_n2 = mdata->local_ny;
     local_y_start = mdata->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)
     {
	  int i2 = i, j2 = j + local_y_start, k2 = k;
	  int index = ((j * n1 + i) * n3 + k);

#  endif /* HAVE_MPI */

#else /* not SCALAR_COMPLEX */

#  ifndef HAVE_MPI

     for (i = 0; i < n_other; ++i)
	  for (j = 0; j < n_last; ++j)
     {
	  int 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 = mdata->local_ny;
     local_y_start = mdata->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 = mdata->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 index = ((j * n1 + i) * local_n3 + k);

#  endif /* HAVE_MPI */

#endif /* not SCALAR_COMPLEX */

	  {
	       vector3 p;
	       int n;
	       p.x = i2 * s1 - c1; p.y = j2 * s2 - c2; p.z = k2 * s3 - c3;
	       for (n = objects.num_items - 1; n >= 0; --n)
		    if (point_in_periodic_fixed_objectp(p, objects.items[n])) {
			 if (objects.items[n].material.which_subclass
			     == MATERIAL_TYPE_SELF)
			      break; /* treat as a "nothing" object */
			 energy_sum += energy[index];
#ifndef SCALAR_COMPLEX
			 {
			      int last_index;
#  ifdef HAVE_MPI
			      if (n3 == 1)
				   last_index = j + local_y_start;
			      else
				   last_index = k;
#  else
			      last_index = j;
#  endif
			      if (last_index != 0 && 2*last_index != last_dim)
				   energy_sum += energy[index];
			 }
#endif
			 break;
		    }
	  }
     }

     mpi_allreduce_1(&energy_sum, real, SCALAR_MPI_TYPE,
		     MPI_SUM, MPI_COMM_WORLD);
     energy_sum *= Vol / H.N;
     return energy_sum;
}

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

/* Compute the integral of f(energy/field, epsilon, r) over the cell. */
cnumber compute_field_integral(function f)
{
     int i, j, k, n1, n2, n3, n_other, n_last, rank, last_dim;
#ifdef HAVE_MPI
     int local_n2, local_y_start, local_n3;
#endif
     real s1, s2, s3, c1, c2, c3;
     int integrate_energy;
     real *energy = (real *) curfield;
     cnumber integral = {0,0};
     vector3 kvector = {0,0,0};

     if (!curfield || !strchr("dheDHRcv", curfield_type)) {
          mpi_one_fprintf(stderr, "The D or H energy/field must be loaded first.\n");
          return integral;
     }
     if (curfield_type != 'v')
	  kvector = cur_kvector;

     integrate_energy = strchr("DHR", curfield_type) != NULL;

     n1 = mdata->nx; n2 = mdata->ny; n3 = mdata->nz;
     n_other = mdata->other_dims;
     n_last = mdata->last_dim_size / (sizeof(scalar_complex)/sizeof(scalar));
     last_dim = mdata->last_dim;
     rank = (n3 == 1) ? (n2 == 1 ? 1 : 2) : 3;

     s1 = geometry_lattice.size.x / n1;
     s2 = geometry_lattice.size.y / n2;
     s3 = geometry_lattice.size.z / n3;
     c1 = n1 <= 1 ? 0 : geometry_lattice.size.x * 0.5;
     c2 = n2 <= 1 ? 0 : geometry_lattice.size.y * 0.5;
     c3 = n3 <= 1 ? 0 : geometry_lattice.size.z * 0.5;

     /* 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 "index" describing the corresponding index in 
	the curfield array.

        This was all stolen from maxwell_eps.c...it would be better
        if we didn't have to cut and paste, sigh. */

#ifdef SCALAR_COMPLEX

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

#  else /* HAVE_MPI */

     local_n2 = mdata->local_ny;
     local_y_start = mdata->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)
     {
	  int i2 = i, j2 = j + local_y_start, k2 = k;
	  int index = ((j * n1 + i) * n3 + k);

#  endif /* HAVE_MPI */

#else /* not SCALAR_COMPLEX */

#  ifndef HAVE_MPI

     for (i = 0; i < n_other; ++i)
	  for (j = 0; j < n_last; ++j)
     {
	  int 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 = mdata->local_ny;
     local_y_start = mdata->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 = mdata->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 index = ((j * n1 + i) * local_n3 + k);

#  endif /* HAVE_MPI */

#endif /* not SCALAR_COMPLEX */

	  {
	       real epsilon;
	       vector3 p;

	       epsilon = mean_epsilon(mdata->eps_inv + index);
	       
	       p.x = i2 * s1 - c1; p.y = j2 * s2 - c2; p.z = k2 * s3 - c3;
	       if (integrate_energy) {
		    integral.re +=
			 ctl_convert_number_to_c(
			      gh_call3(f,
				     ctl_convert_number_to_scm(energy[index]),
				       ctl_convert_number_to_scm(epsilon),
				       ctl_convert_vector3_to_scm(p)));	  
	       }
	       else {
		    cvector3 F;
		    double phase_phi;
		    scalar_complex phase;
		    cnumber integrand;

		    phase_phi = TWOPI * 
			 (kvector.x * (p.x/geometry_lattice.size.x+0.5) +
			  kvector.y * (p.y/geometry_lattice.size.y+0.5) +
			  kvector.z * (p.z/geometry_lattice.size.z+0.5));
		    CASSIGN_SCALAR(phase, cos(phase_phi), sin(phase_phi));
		    CASSIGN_MULT_RE(F.x.re, curfield[3*index+0], phase);
		    CASSIGN_MULT_IM(F.x.im, curfield[3*index+0], phase);
		    CASSIGN_MULT_RE(F.y.re, curfield[3*index+1], phase);
		    CASSIGN_MULT_IM(F.y.im, curfield[3*index+1], phase);
		    CASSIGN_MULT_RE(F.z.re, curfield[3*index+2], phase);
		    CASSIGN_MULT_IM(F.z.im, curfield[3*index+2], phase);
		    integrand =
			 ctl_convert_cnumber_to_c(
                              gh_call3(f,
				       ctl_convert_cvector3_to_scm(F),
				       ctl_convert_number_to_scm(epsilon),
                                       ctl_convert_vector3_to_scm(p)));
		    integral.re += integrand.re;
		    integral.im += integrand.im;
	       }

#ifndef SCALAR_COMPLEX
	       {
		    int last_index;
#  ifdef HAVE_MPI
		    if (n3 == 1)
			 last_index = j + local_y_start;
		    else
			 last_index = k;
#  else
		    last_index = j;
#  endif
		    
		    if (last_index != 0 && 2*last_index != last_dim) {
			 int i2c, j2c, k2c;
			 i2c = i2 ? (n1 - i2) : 0;
			 j2c = j2 ? (n2 - j2) : 0;
			 k2c = k2 ? (n3 - k2) : 0;
			 p.x = i2c * s1 - c1; 
			 p.y = j2c * s2 - c2; 
			 p.z = k2c * s3 - c3;
			 if (integrate_energy)
			      integral.re += 
				   ctl_convert_number_to_c(
					gh_call3(f,
				      ctl_convert_number_to_scm(energy[index]),
				      ctl_convert_number_to_scm(epsilon),
				      ctl_convert_vector3_to_scm(p)));
			 else {
			      cvector3 F;
			      double phase_phi;
			      scalar_complex phase, Fx, Fy, Fz;
			      cnumber integrand;
			      
			      Fx = curfield[3*index+0];
			      Fy = curfield[3*index+1];
			      Fz = curfield[3*index+2];
			      Fx.im= -Fx.im; Fy.im= -Fy.im; Fz.im= -Fz.im;

			      phase_phi = TWOPI * 
				   (kvector.x 
				    * (p.x/geometry_lattice.size.x+0.5) +
				    kvector.y 
				    * (p.y/geometry_lattice.size.y+0.5) +
				    kvector.z 
				    * (p.z/geometry_lattice.size.z+0.5));
			      CASSIGN_SCALAR(phase, 
					     cos(phase_phi), sin(phase_phi));
			      CASSIGN_MULT_RE(F.x.re, Fx, phase);
			      CASSIGN_MULT_IM(F.x.im, Fx, phase);
			      CASSIGN_MULT_RE(F.y.re, Fy, phase);
			      CASSIGN_MULT_IM(F.y.im, Fy, phase);
			      CASSIGN_MULT_RE(F.z.re, Fz, phase);
			      CASSIGN_MULT_IM(F.z.im, Fz, phase);

			      integrand =
				   ctl_convert_cnumber_to_c(
					gh_call3(f,
					    ctl_convert_cvector3_to_scm(F),
				            ctl_convert_number_to_scm(epsilon),
                                            ctl_convert_vector3_to_scm(p)));
			      integral.re += integrand.re;
			      integral.im += integrand.im;
			 }
		    }
	       }
#endif
	  }
     }

     integral.re *= Vol / H.N;
     integral.im *= Vol / H.N;
     {
	  cnumber integral_sum;
	  mpi_allreduce(&integral, &integral_sum, 2, number, 
			MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
	  return integral_sum;
     }
}

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

     return cnumber_re(compute_field_integral(f));
}

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