monitor.cpp

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

/* Copyright (C) 2006 Massachusetts Institute of Technology
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2, or (at your option)
%  any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software Foundation,
%  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "meep.hpp"
#include "meep_internals.hpp"

#include "config.h"
#if defined(HAVE_LIBFFTW3)
#  include <fftw3.h>
#elif defined(HAVE_LIBDFFTW)
#  include <dfftw.h>
#elif defined(HAVE_LIBFFTW)
#  include <fftw.h>
#endif
#define HAVE_SOME_FFTW (defined(HAVE_LIBFFTW3) || defined(HAVE_LIBFFTW) || defined(HAVE_LIBDFFTW))

/* Below are the monitor point routines. */

namespace meep {

monitor_point::monitor_point() {
  next = NULL;
}

monitor_point::~monitor_point() {
  if (next) delete next;
}

inline complex<double> getcm(const double * const f[2], int i) {
  return complex<double>(f[0][i],f[1][i]);
}

static void dumbsort(complex<double> val[8]) {
  for (int i=0;i<7;i++) {
    int lowest = i;
    for (int j=i+1;j<8;j++) if (abs(val[j]) < abs(val[lowest])) lowest = j;
    complex<double> tmp = val[i];
    val[i] = val[lowest];
    val[lowest] = tmp;
  }
}

static void dumbsort(double val[8]) {
  for (int i=0;i<7;i++) {
    int lowest = i;
    for (int j=i+1;j<8;j++) if (abs(val[j]) < abs(val[lowest])) lowest = j;
    double tmp = val[i];
    val[i] = val[lowest];
    val[lowest] = tmp;
  }
}

void fields::get_point(monitor_point *pt, const vec &loc) const {
  if (pt == NULL) abort("Error:  get_point passed a null pointer!\n");
  for (int i=0;i<10;i++) pt->f[i] = 0.0;
  pt->loc = loc;
  pt->t = time();
  FOR_COMPONENTS(c)
    if (v.has_field(c))
      pt->f[c] = get_field(c,loc);
}

complex<double> fields::get_field(int c, const vec &loc) const {
  return (is_derived(c) ? get_field(derived_component(c), loc)
	  : get_field(component(c), loc));
}

double fields::get_field(derived_component c, const vec &loc) const {
  component c1, c2;
  double sum = 0;
  switch (c) {
  case Sx: case Sy: case Sz: case Sr: case Sp:
    switch (c) {
    case Sx: c1 = Ey; c2 = Hz; break;
    case Sy: c1 = Ez; c2 = Hx; break;
    case Sz: c1 = Ex; c2 = Hy; break;
    case Sr: c1 = Ep; c2 = Hz; break;
    case Sp: c1 = Ez; c2 = Hr; break;
    default: break; // never
    }
    sum += real(conj(get_field(c1, loc)) * get_field(c2, loc));
    sum -= real(conj(get_field(direction_component(Ex, component_direction(c2)), loc)) * get_field(direction_component(Hx, component_direction(c1)), loc));
    return sum;
  case EnergyDensity: case D_EnergyDensity: case H_EnergyDensity:
    if (c != H_EnergyDensity)
      FOR_ELECTRIC_COMPONENTS(c1) if (v.has_field(c1)) {
	c2 = direction_component(Dx, component_direction(c1));
	sum += real(conj(get_field(c1, loc)) * get_field(c2, loc));
      }
    if (c != D_EnergyDensity)
      FOR_MAGNETIC_COMPONENTS(c1) if (v.has_field(c1)) {
	complex<double> f = get_field(c1, loc);
	sum += real(conj(f) * f);
      }
    return sum * 0.5;
  default: abort("unknown derived_component in get_field");
  }
}

complex<double> fields::get_field(component c, const vec &loc) const {
  switch (c) {
  case Dielectric:
    return get_eps(loc);
  default:
    ivec ilocs[8];
    double w[8];
    complex<double> val[8];
    for (int i=0;i<8;i++) val[i] = 0.0;
    v.interpolate(c, loc, ilocs, w);
    for (int argh=0;argh<8&&w[argh];argh++)
      val[argh] = w[argh]*get_field(c,ilocs[argh]);
    dumbsort(val);
    complex<double> res = 0.0;
    for (int i=0;i<8;i++) res += val[i];
    return res;
  }
}

complex<double> fields::get_field(component c, const ivec &origloc) const {
  ivec iloc = origloc;
  complex<double> kphase = 1.0;
  locate_point_in_user_volume(&iloc, &kphase);
  for (int sn=0;sn<S.multiplicity();sn++)
    for (int i=0;i<num_chunks;i++)
      if (chunks[i]->v.contains(S.transform(iloc,sn)))
        return S.phase_shift(c,sn)*kphase*
          chunks[i]->get_field(S.transform(c,sn),S.transform(iloc,sn));
  return 0.0;
}

complex<double> fields_chunk::get_field(component c, const ivec &iloc) const {
  complex<double> res = 0.0;
  if (f[c][0] && f[c][1]) res = getcm(f[c], v.index(c, iloc));
  else if (f[c][0]) res = f[c][0][v.index(c,iloc)];
  return broadcast(n_proc(), res);
}

/* Bounding box for zero-communication get_field, below.  This is the
   largest box in which you can interpolate the fields without communication.
   It is *not* necessarily non-overlapping with other chunks. */
geometric_volume fields_chunk::get_field_gv(component c) const {
  switch (c) {
  case Dielectric:
    c = v.eps_component();
  default:
    return geometric_volume(v.loc(c, 0), v.loc(c, v.ntot() - 1));
  }
}

/* Non-collective, zero-communication get_field... loc *must*
   be in get_field_gv(c). */
complex<double> fields_chunk::get_field(component c, const vec &loc) const {
  ivec ilocs[8];
  double w[8];
  switch (c) {
  case Dielectric: {
    c = v.eps_component();
    v.interpolate(c, loc, ilocs, w);
    double res = 0.0;
    for (int i = 0; i < 8 && w[i] != 0.0; ++i) {
      if (!v.contains(ilocs[i]))
	abort("invalid loc in chunk get_field, weight = %g", w[i]);
      res += s->eps[v.index(c,ilocs[i])] * w[i];
    }
    return res;
  }
  default: {
    v.interpolate(c, loc, ilocs, w);
    complex<double> res = 0.0;
    for (int i = 0; i < 8 && w[i] != 0.0; ++i) {
      if (!v.contains(ilocs[i]))
	abort("invalid loc in chunk get_field, weight = %g", w[i]);
      if (f[c][0] && f[c][1]) res += getcm(f[c], v.index(c, ilocs[i])) * w[i];
      else if (f[c][0]) res += f[c][0][v.index(c,ilocs[i])] * w[i];
    }
    return res;
  }
  }
}

complex<double> fields_chunk::get_polarization_field(const polarizability_identifier &p,
                                                     component c, const ivec &iloc) const {
  complex<double> res = 0.0;
  for (polarization *P = olpol; P; P = P->next)
    if (P->pb->get_identifier() == p) {
      if (P->P[c][0] && P->P[c][1]) res = getcm(P->P[c], v.index(c, iloc));
      else if (P->P[c][0]) res = P->P[c][0][v.index(c,iloc)];
    }
  return broadcast(n_proc(), res);
}

double fields::get_polarization_energy(const vec &loc) const {
  ivec ilocs[8];
  double w[8], val[8];
  for (int i=0;i<8;i++) val[i] = 0.0;
  v.interpolate(v.eps_component(), loc, ilocs, w);
  for (int argh=0;argh<8&&w[argh];argh++)
    val[argh] = w[argh]*get_polarization_energy(ilocs[argh]);
  dumbsort(val);
  double res = 0.0;
  for (int i=0;i<8;i++) res += val[i];
  return res;
}

double fields::get_polarization_energy(const ivec &origloc) const {
  ivec iloc = origloc;
  complex<double> aaack = 1.0;
  locate_point_in_user_volume(&iloc, &aaack);
  for (int sn=0;sn<S.multiplicity();sn++)
    for (int i=0;i<num_chunks;i++)
      if (chunks[i]->v.contains(S.transform(iloc,sn)))
        return chunks[i]->get_polarization_energy(S.transform(iloc,sn));
  return 0.0;
}

double fields_chunk::get_polarization_energy(const ivec &iloc) const {
  double res = 0.0;
  polarization *p = pol;
  while (is_mine() && p) {
    res += p->local_energy(iloc);
    p = p->next;
  }
  return broadcast(n_proc(), res);
}

double fields_chunk::my_polarization_energy(const ivec &iloc) const {
  if (!is_mine())
    abort("Can't call my_polarization_energy on someone else's chunk!\n");
  double res = 0.0;
  polarization *p = pol;
  while (p) {
    res += p->local_energy(iloc);
    p = p->next;
  }
  return res;
}

double fields::get_polarization_energy(const polarizability_identifier &p,
                                       const vec &loc) const {
  ivec ilocs[8];
  double w[8], val[8];
  for (int i=0;i<8;i++) val[i] = 0.0;
  v.interpolate(v.eps_component(), loc, ilocs, w);
  for (int argh=0;argh<8&&w[argh];argh++)
    val[argh] = w[argh]*get_polarization_energy(p, ilocs[argh]);
  dumbsort(val);
  double res = 0.0;
  for (int i=0;i<8;i++) res += val[i];
  return res;
}

double fields::get_polarization_energy(const polarizability_identifier &p,
                                       const ivec &origloc) const {
  ivec iloc = origloc;
  complex<double> aaack = 1.0;
  locate_point_in_user_volume(&iloc, &aaack);
  for (int sn=0;sn<S.multiplicity();sn++)
    for (int i=0;i<num_chunks;i++)
      if (chunks[i]->v.contains(S.transform(iloc,sn)))
        return chunks[i]->get_polarization_energy(p, S.transform(iloc,sn));
  return 0.0;
}

double fields_chunk::get_polarization_energy(const polarizability_identifier &pi,