symmetry.cpp

采用FDTD时域有限差分法计算电磁波的传播问题等。

/* 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 <signal.h>

#include <meep.hpp>
using namespace meep;

const char *mydirname = "symmetry-out";

double one(const vec &) { return 1.0; }
vec the_center;
double rods_2d(const vec &pp) {
  vec p = pp - the_center;
  while (p.x() > 0.5) p -= vec(1.0,0);
  while (p.x() <-0.5) p += vec(1.0,0);
  while (p.y() > 0.5) p -= vec(0,1.0);
  while (p.y() <-0.5) p += vec(0,1.0);
  if (fabs(p.x()) < 0.314) return 12.0;
  if (fabs(p.y()) < 0.314) return 12.0;
  return 1.0;
}

static const double eps_compare = 7e-15;
static const double thresh_compare = 1e-16;

static inline double max(double a, double b) { return a > b ? a : b; }

int compare(double a, double b, const char *n) {
  if (fabs(a-b) > fabs(b)*eps_compare
      && max(fabs(a),fabs(b)) > thresh_compare) {
    master_printf("%s differs by\t%g out of\t%g\n", n, a-b, b);
    master_printf("This gives a fractional error of %g\n", fabs(a-b)/fabs(b));
    return 0;
  } else {
    return 1;
  }
}

int compare_point(fields &f1, fields &f2, const vec &p) {
  monitor_point m1, m_test;
  f1.get_point(&m_test, p);
  f2.get_point(&m1, p);
  for (int i=0;i<10;i++) {
    component c = (component) i;
    if (f1.v.has_field(c)) {
      complex<double> v1 = m_test.get_component(c), v2 = m1.get_component(c);
      if (abs(v1 - v2) > eps_compare*abs(v2) && abs(v2) > thresh_compare) {
        master_printf("%s differs:  %g %g out of %g %g\n",
               component_name(c), real(v2-v1), imag(v2-v1), real(v2), imag(v2));
        master_printf("This comes out to a fractional error of %g\n",
               abs(v1 - v2)/abs(v2));
        master_printf("Right now I'm looking at ");
        LOOP_OVER_DIRECTIONS(p.dim,d)
          master_printf("%s = %g, ", direction_name(d), p.in_direction(d));
        master_printf("time %g\n", f1.time());
        f1.output_real_imaginary_slices("multi");
        f2.output_real_imaginary_slices("single");
        f1.eps_slices("multi");
        f2.eps_slices("single");
        return 0;
      }
    }
  }
  return 1;
}

void check_unequal_layout(const fields &f1, const fields &f2)
{
  if (f1.equal_layout(f2) ||
      !f1.equal_layout(f1) ||
      !f2.equal_layout(f2))
    abort("fields::equal_layout did not return expected result");
}

int test_cyl_metal_mirror(double eps(const vec &)) {
  master_printf("Testing Z mirror symmetry in Cylindrical...\n");
  double a = 8.0;
  double ttot = 3.0;

  const volume v = volcyl(1.0, 1.0, a);
  the_center = v.center();
  const symmetry S = mirror(Z,v);
  structure s(v, eps, no_pml(), S);
  structure s1(v, eps);

  fields f1(&s1);
  f1.add_point_source(Er, 0.7, 2.5, 0.0, 4.0, veccyl(0.5,0.5));
  f1.add_point_source(Ep, 0.8, 0.6, 0.0, 4.0, veccyl(0.401,0.5));
  fields f(&s);
  f.add_point_source(Er, 0.7, 2.5, 0.0, 4.0, veccyl(0.5,0.5));
  f.add_point_source(Ep, 0.8, 0.6, 0.0, 4.0, veccyl(0.401,0.5));
  check_unequal_layout(f, f1);
  double total_energy_check_time = 1.0;
  while (f.time() < ttot) {
    f.step();
    f1.step();
    if (!compare_point(f, f1, veccyl(0.01,  0.5  ))) return 0;
    if (!compare_point(f, f1, veccyl(0.21,  0.5  ))) return 0;
    if (!compare_point(f, f1, veccyl(0.501, 0.5  ))) return 0;
    if (!compare_point(f, f1, veccyl(0.33,  0.46 ))) return 0;
    if (!compare_point(f, f1, veccyl(0.2,   0.2  ))) return 0;
    if (f.time() >= total_energy_check_time) {
      if (!compare(f.electric_energy_in_box(v.surroundings()),
                   f1.electric_energy_in_box(v.surroundings()),
                   "electric energy")) return 0;
      if (!compare(f.magnetic_energy_in_box(v.surroundings()),
                   f1.magnetic_energy_in_box(v.surroundings()),
                   "magnetic energy")) return 0;
      if (!compare(f.total_energy(), f1.total_energy(),
                   "   total energy")) return 0;
      total_energy_check_time += 1.0;
    }
  }
  return 1;
}

int test_cyl_metal_mirror_nonlinear(double eps(const vec &)) {
  master_printf("Testing Z mirror symmetry in Cylindrical...\n");
  double a = 16.0;
  double ttot = 3.0;

  const volume v = volcyl(1.0, 1.0, a);
  the_center = v.center();
  const symmetry S = mirror(Z,v);
  structure s(v, eps, no_pml(), S);
  structure s1(v, eps);
  s.set_kerr(one);
  s1.set_kerr(one);

  fields f1(&s1);
  f1.add_point_source(Er, 0.7, 2.5, 0.0, 4.0, veccyl(0.5,0.5));
  f1.add_point_source(Ep, 0.8, 0.6, 0.0, 4.0, veccyl(0.401,0.5));
  fields f(&s);
  f.add_point_source(Er, 0.7, 2.5, 0.0, 4.0, veccyl(0.5,0.5));
  f.add_point_source(Ep, 0.8, 0.6, 0.0, 4.0, veccyl(0.401,0.5));
  check_unequal_layout(f, f1);
  double total_energy_check_time = 1.0;
  while (f.time() < ttot) {
    f.step();
    f1.step();
    if (!compare_point(f, f1, veccyl(0.01,  0.5  ))) return 0;
    if (!compare_point(f, f1, veccyl(0.21,  0.5  ))) return 0;
    if (!compare_point(f, f1, veccyl(0.501, 0.5  ))) return 0;
    if (!compare_point(f, f1, veccyl(0.33,  0.46 ))) return 0;
    if (!compare_point(f, f1, veccyl(0.2,   0.2  ))) return 0;
    if (f.time() >= total_energy_check_time) {
      if (!compare(f.electric_energy_in_box(v.surroundings()),
                   f1.electric_energy_in_box(v.surroundings()),
                   "electric energy")) return 0;
      if (!compare(f.magnetic_energy_in_box(v.surroundings()),
                   f1.magnetic_energy_in_box(v.surroundings()),
                   "magnetic energy")) return 0;
      if (!compare(f.total_energy(), f1.total_energy(),
                   "   total energy")) return 0;
      total_energy_check_time += 1.0;
    }
  }
  return 1;
}

int test_1d_periodic_mirror(double eps(const vec &)) {
  master_printf("Testing Z mirror symmetry in 1D...\n");
  double a = 16.0;
  double ttot = 3.0;

  const volume v = volone(1.0, a);
  the_center = v.center();
  const symmetry S = mirror(Z,v);
  structure s(v, eps, no_pml(), S);
  structure s1(v, eps);

  fields f1(&s1);
  f1.use_bloch(0.0);
  f1.add_point_source(Ex, 0.7, 2.5, 0.0, 4.0, vec(0.5));
  fields f(&s);
  f.use_bloch(0.0);
  f.add_point_source(Ex, 0.7, 2.5, 0.0, 4.0, vec(0.5));
  check_unequal_layout(f, f1);
  double total_energy_check_time = 1.0;
  while (f.time() < ttot) {
    f.step();
    f1.step();
    if (!compare_point(f, f1, vec(0.01))) return 0;
    if (!compare_point(f, f1, vec(0.33))) return 0;
    if (!compare_point(f, f1, vec(0.50))) return 0;
    if (f.time() >= total_energy_check_time) {
      if (!compare(f.electric_energy_in_box(v.surroundings()),
                   f1.electric_energy_in_box(v.surroundings()),
                   "electric energy")) return 0;
      if (!compare(f.magnetic_energy_in_box(v.surroundings()),
                   f1.magnetic_energy_in_box(v.surroundings()),
                   "magnetic energy")) return 0;
      if (!compare(f.total_energy(), f1.total_energy(),
                   "   total energy")) return 0;
      total_energy_check_time += 1.0;
    }
  }
  return 1;
}

int test_origin_shift(const char *mydirname) {
  master_printf("Testing origin shift in 2D...\n");
  double a = 8.0;
  double ttot = 3.0;

  const volume v = voltwo(1.0, 1.0, a);
  volume vcentered = v;
  vcentered.shift_origin(-v.center());
  structure s(vcentered, one);
  structure s1(v, one);
  s.set_output_directory(mydirname);
  s1.set_output_directory(mydirname);

  fields f1(&s1);
  fields f(&s);
  f1.add_point_source(Ey, 0.7, 2.5, 0.0, 4.0, v.center());
  f1.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, v.center());
  f.add_point_source(Ey, 0.7, 2.5, 0.0, 4.0, vec(0.0,0.0));
  f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(0.0,0.0));
  check_unequal_layout(f, f1);
  while (f.time() < ttot) {
    f.step();
    f1.step();
    if (!compare(f.total_energy(), f1.total_energy(), "   total energy")) {
      master_printf("Time is %g\n", f.time());
      f1.output_real_imaginary_slices("unshifted");
      f.output_real_imaginary_slices("shifted");
      f1.eps_slices("unshifted");
      f.eps_slices("shifted");
      return 0;
    }
  }
  return 1;
}

int test_metal_xmirror(double eps(const vec &)) {
  master_printf("Testing X mirror symmetry...\n");
  double a = 8.0;
  double ttot = 3.0;

  const volume v = voltwo(1.0, 1.0, a);
  the_center = v.center();
  const symmetry S = mirror(X,v);
  structure s(v, eps, no_pml(), S);
  structure s1(v, eps);

  fields f1(&s1);
  f1.add_point_source(Ey, 0.7, 2.5, 0.0, 4.0, vec(0.5,0.5));
  f1.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(0.5,0.401));
  fields f(&s);
  f.add_point_source(Ey, 0.7, 2.5, 0.0, 4.0, vec(0.5,0.5));
  f.add_point_source(Ez, 0.8, 0.6, 0.0, 4.0, vec(0.5,0.401));
  check_unequal_layout(f, f1);
  double total_energy_check_time = 1.0;
  while (f.time() < ttot) {
    f.step();
    f1.step();
    if (!compare_point(f, f1, vec(0.5  , 0.01))) return 0;
    if (!compare_point(f, f1, vec(0.5  , 0.21))) return 0;
    if (!compare_point(f, f1, vec(0.5  , 0.501))) return 0;
    if (!compare_point(f, f1, vec(0.46 , 0.33))) return 0;
    if (!compare_point(f, f1, vec(0.2  , 0.2 ))) return 0;
    if (f.time() >= total_energy_check_time) {
      if (!compare(f.electric_energy_in_box(v.surroundings()),
                   f1.electric_energy_in_box(v.surroundings()),
                   "electric energy")) return 0;
      if (!compare(f.magnetic_energy_in_box(v.surroundings()),
                   f1.magnetic_energy_in_box(v.surroundings()),
                   "magnetic energy")) return 0;
      if (!compare(f.total_energy(), f1.total_energy(),
                   "   total energy")) return 0;
      total_energy_check_time += 1.0;
    }
  }
  return 1;
}

int test_3D_metal_xmirror(double eps(const vec &)) {
  double a = 8.0;
  double ttot = 3.0;

  const volume v = vol3d(1.0, 1.0, 1.0, a);
  const symmetry S = mirror(X,v);
  structure s(v, eps, no_pml(), S);
  structure s1(v, eps);
  master_printf("Testing X mirror symmetry in 3D...\n");

  fields f1(&s1);
  f1.add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(0.5,0.51,0.55));
  f1.add_point_source(Hx, 0.8, 0.6, 0.0, 4.0, vec(0.5,0.401,0.43));
  fields f(&s);
  f.add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(0.5,0.51,0.55));
  f.add_point_source(Hx, 0.8, 0.6, 0.0, 4.0, vec(0.5,0.401,0.43));
  check_unequal_layout(f, f1);
  double total_energy_check_time = 1.0;
  while (f.time() < ttot) {
    f.step();
    f1.step();
    if (!compare_point(f, f1, vec(0.5  , 0.01 , 0.5))) return 0;
    if (!compare_point(f, f1, vec(0.5  , 0.21 , 0.5))) return 0;
    if (!compare_point(f, f1, vec(0.5  , 0.501, 0.5))) return 0;
    if (!compare_point(f, f1, vec(0.46 , 0.33 , 0.5))) return 0;
    if (!compare_point(f, f1, vec(0.2  , 0.2  , 0.5))) return 0;
    if (f.time() >= total_energy_check_time) {
      if (!compare(f.electric_energy_in_box(v.surroundings()),
                   f1.electric_energy_in_box(v.surroundings()),
                   "electric energy")) return 0;
      if (!compare(f.magnetic_energy_in_box(v.surroundings()),
                   f1.magnetic_energy_in_box(v.surroundings()),
                   "magnetic energy")) return 0;
      if (!compare(f.total_energy(), f1.total_energy(),
                   "   total energy")) return 0;
      total_energy_check_time += 1.0;
    }
  }
  return 1;
}

int test_3D_metal_zmirror(double eps(const vec &)) {
  double a = 8.0;
  double ttot = 3.0;

  const volume v = vol3d(1.1, 0.6, 1.0, a);
  const symmetry S = mirror(Z,v);
  structure s(v, eps, no_pml(), S);
  structure s1(v, eps);
  master_printf("Testing Z mirror symmetry in 3D...\n");

  fields f1(&s1);
  f1.add_point_source(Ex, 0.7, 2.5, 0.0, 4.0, vec(0.55,0.51,0.5));
  f1.add_point_source(Ey, 0.8, 0.6, 0.0, 4.0, vec(0.43,0.401,0.5));
  fields f(&s);
  f.add_point_source(Ex, 0.7, 2.5, 0.0, 4.0, vec(0.55,0.51,0.5));
  f.add_point_source(Ey, 0.8, 0.6, 0.0, 4.0, vec(0.43,0.401,0.5));
  check_unequal_layout(f, f1);
  double total_energy_check_time = 1.0;
  while (f.time() < ttot) {
    f.step();
    f1.step();
    if (!compare_point(f, f1, vec(0.5  , 0.01 , 0.75))) return 0;
    if (!compare_point(f, f1, vec(0.5  , 0.21 , 0.15))) return 0;
    if (!compare_point(f, f1, vec(0.5  , 0.501, 0.5))) return 0;
    if (!compare_point(f, f1, vec(0.46 , 0.33 , 0.51))) return 0;
    if (!compare_point(f, f1, vec(0.2  , 0.2  , 0.05))) return 0;
    if (f.time() >= total_energy_check_time) {
      if (!compare(f.electric_energy_in_box(v.surroundings()),
                   f1.electric_energy_in_box(v.surroundings()),
                   "electric energy")) return 0;
      if (!compare(f.magnetic_energy_in_box(v.surroundings()),
                   f1.magnetic_energy_in_box(v.surroundings()),
                   "magnetic energy")) return 0;
      if (!compare(f.total_energy(), f1.total_energy(),
                   "   total energy")) return 0;
      total_energy_check_time += 1.0;
    }