known_results.cpp

来自mit的fdtd开放性源代码meep

/* Copyright (C) 2005-2008 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 <meep.hpp>
using namespace meep;

#include "config.h"

double one(const vec &) { return 1.0; }

double rods(const vec &r) {
  vec p = r;
  while (p.x() < -0.5) p.set_direction(X, p.x() + 1.0);
  while (p.x() >  0.5) p.set_direction(X, p.x() - 1.0);
  while (p.y() < -0.5) p.set_direction(Y, p.y() + 1.0);
  while (p.y() >  0.5) p.set_direction(Y, p.y() - 1.0);
  if (p.x()*p.x() + p.y()*p.y() < 0.3) return 12.0;
  return 1.0;
}

void compare(double b, double a, const char *n) {
  if (fabs(a-b) > fabs(b)*1e-5 || b != b) {
    abort("Failed %s (%g instead of %g, relerr %0.2g)\n", n,
	  a, b, fabs(a-b)/fabs(b));
  } else {
    master_printf("Passed %s\n", n);
  }
}

static double dpml = 1.0;

double using_pml_ez(const volume &v, double eps(const vec &)) {
  const double ttot = 30.0;
  structure s(v, eps, pml(dpml));
  fields f(&s);
  f.add_point_source(Ez, 0.2, 3.0, 0.0, 2.0, v.center(),
		     complex<double>(0,-2*pi*0.2));
  while (f.time() < ttot) f.step();
  monitor_point p;
  f.get_point(&p, v.center());
  return real(p.get_component(Ez));
}

double x_periodic_y_pml(const volume &v, double eps(const vec &)) {
  const double ttot = 30.0;
  structure s(v, eps, pml(dpml, Y));
  fields f(&s);
  f.add_point_source(Ez, 0.2, 3.0, 0.0, 2.0, v.center(),
		     complex<double>(0,-2*pi*0.2));
  f.use_bloch(X, 0.1);
  while (f.time() < ttot) f.step();
  monitor_point p;
  f.get_point(&p, v.center());
  return real(p.get_component(Ez));
}

double x_periodic(const volume &v, double eps(const vec &)) {
  const double ttot = 30.0;
  structure s(v, eps);
  fields f(&s);
  f.add_point_source(Ez, 0.2, 3.0, 0.0, 2.0, v.center(),
		     complex<double>(0,-2*pi*0.2));
  f.use_bloch(X, 0.1);
  while (f.time() < ttot) f.step();
  monitor_point p;
  f.get_point(&p, v.center());
  return real(p.get_component(Ez));
}

double periodic_ez(const volume &v, double eps(const vec &)) {
  const double ttot = 30.0;
  structure s(v, eps);
  fields f(&s);
  f.add_point_source(Ez, 0.2, 3.0, 0.0, 2.0, v.center(),
		     complex<double>(0,-2*pi*0.2));
  vec k;
  switch (v.dim) {
  case D1: k = vec(0.3); break;
  case D2: k = vec(0.3,0.4); break;
  case D3: k = vec(0.3,0.5,0.8); break;
  case Dcyl: k = veccyl(0.3,0.2); break;
  }
  f.use_bloch(k);
  while (f.time() < ttot) f.step();
  monitor_point p;
  f.get_point(&p, v.center());
  return real(p.get_component(Ez));
}

double metallic_ez(const volume &v, double eps(const vec &)) {
  const double ttot = 10.0;
  structure s(v, eps);
  fields f(&s);
  f.add_point_source(Ez, 0.2, 3.0, 0.0, 2.0, v.center(),
		     complex<double>(0,-2*pi*0.2));
  while (f.time() < ttot) f.step();
  monitor_point p;
  f.get_point(&p, v.center());
  return real(p.get_component(Ez));
}

double polariton_ex(const volume &v, double eps(const vec &)) {
  const double ttot = 10.0;
  structure s(v, eps);
  s.add_polarizability(one, 0.3, 0.1, 7.63);
  fields f(&s);
  f.add_point_source(Ex, 0.2, 3.0, 0.0, 2.0, v.center(),
		     complex<double>(0,-2*pi*0.2));
  while (f.time() < ttot) f.step();
  monitor_point p;
  f.get_point(&p, v.center());
  return real(p.get_component(Ex));
}

double polariton_energy(const volume &v, double eps(const vec &)) {
  const double ttot = 10.0;
  structure s(v, eps);
  s.add_polarizability(one, 0.3, 0.1, 7.63);
  fields f(&s, 0, 1);
  f.add_point_source(Ex, 0.2, 3.0, 0.0, 2.0, v.center(),
		     complex<double>(0,-2*pi*0.2));
  while (f.time() < ttot) f.step();
  return f.total_energy();
}

double saturated_polariton_ex(const volume &v, double eps(const vec &)) {
  const double ttot = 10.0;
  structure s(v, eps);
  polarizability_identifier thep = s.add_polarizability(one, 0.3, 0.1, -0.063, 0.1);
  fields f(&s);
  f.use_real_fields();
  f.add_point_source(Ex, 0.2, 3.0, 0.0, 2.0, v.center(), 
		     1/sqrt(4*pi) * complex<double>(0,-2*pi*0.2));
  while (f.time() < ttot) f.step();
  monitor_point p;
  f.get_point(&p, v.center());
  return real(p.get_component(Ex)) * sqrt(4*pi);
}

int main(int argc, char **argv) {
  initialize mpi(argc, argv);
  quiet = true;
  const char *mydirname = "known_results-out";
  trash_output_directory(mydirname);
  master_printf("Testing with some known results...\n");
  const double a = 10.0;

  compare(-0.0894851, polariton_ex(volone(1.0, a), one),
          "1D polariton");
#ifdef WITH_SATURABLE_ABSORBERS
  compare(-0.0384049, saturated_polariton_ex(volone(1.0, a), one),
          "1D saturated polariton");
  if (count_processors() <= 2) // FIXME: broken for NCPUS > 2 ???
    compare(-7.57915, saturated_polariton_ex(vol2d(1.0,1.0, a), one),
	    "2D saturated polariton");
  compare(-23.8506, saturated_polariton_ex(vol3d(1.0,1.0,0.5, a), one),
          "3D saturated polariton");
#endif
  compare(0.32617294, polariton_energy(volone(1.0, a), one),
          "1D polariton energy");
  compare(5.20605, metallic_ez(voltwo(1.0, 1.0, a), one),
          "1x1 metallic 2D TM");
  compare(0.883776, using_pml_ez(voltwo(1.0+2*dpml, 1.0+2*dpml, a), one),
          "1x1 PML 2D TM");
  compare(0.110425, x_periodic(voltwo(1.0, 1.0, a), one),
          "1x1 X periodic 2D TM");
  compare(-4.66027, periodic_ez(voltwo(1.0, 3.0, a), rods),
          "1x1 fully periodic 2D TM rods");
  compare(1.12502, periodic_ez(voltwo(1.0, 3.0, a), one),
          "1x1 fully periodic 2D TM");
  compare(0.608815, x_periodic_y_pml(voltwo(1.0, 1.0+2*dpml, a), one),
          "1x1 X periodic Y PML 2D TM");
  compare(-41.8057, metallic_ez(vol3d(1.0, 1.0, 1.0, a), one),
          "1x1x1 metallic 3D");
  compare(-100.758, x_periodic(vol3d(1.0, 1.0, 1.0, a), one),
          "1x1x1 X periodic 3D");
  compare(-101.398, x_periodic_y_pml(vol3d(1.0, 1.0+2*dpml, 1.0, a), one),
          "1x1x1 X periodic Y PML 3D");
  compare(-97.7989, periodic_ez(vol3d(1.0, 1.0, 1.0, a), rods),
          "1x1x1 fully periodic 3D rods");
  compare(-99.1618, periodic_ez(vol3d(1.0, 1.0, 1.0, a), one),
          "1x1x1 fully periodic 3D");

  return 0;
}