symmetry.cpp

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

    if (f_mm.time() >= total_energy_check_time) {
      if (!compare(f_mm.electric_energy_in_box(v.surroundings()),
                   f1.electric_energy_in_box(v.surroundings()),
                   "electric energy")) return 0;
      total_energy_check_time += 3.0;
    }
  }
  return 1;
}

int exact_metal_rot4z(double eps(const vec &)) {
  double a = 8.0;
  double ttot = 5.0;

  const volume v = voltwo(1.0, 1.0, a);
  the_center = v.center();
  const symmetry S = rotate4(Z,v);

  structure s(v, eps, no_pml(), S);
  structure s1(v, eps);
  master_printf("Testing Z fourfold rotational symmetry...\n");

  fields f1(&s1);
  f1.add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(0.5,0.5));
  f1.add_point_source(Hz, 0.8, 0.6, 0.0, 4.0, vec(0.5,0.5));
  fields f(&s);
  f.add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(0.5,0.5));
  f.add_point_source(Hz, 0.8, 0.6, 0.0, 4.0, vec(0.5,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 ,  0.5))) return 0;
    if (!compare_point(f, f1, vec(0.21 ,  0.5))) 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 exact_metal_rot4z_nonlinear(double eps(const vec &)) {
  double a = 8.0;
  double ttot = 5.0;

  const volume v = voltwo(1.0, 1.0, a);
  the_center = v.center();
  const symmetry S = rotate4(Z,v);

  structure s(v, eps, no_pml(), S);
  structure s1(v, eps);
  s.set_chi3(one);
  s1.set_chi3(one);
  master_printf("Testing nonlinear Z fourfold rotational symmetry...\n");

  fields f1(&s1);
  //f1.add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(0.5,0.5));
  f1.add_point_source(Hz, 0.8, 0.6, 0.0, 4.0, vec(0.5,0.5));
  fields f(&s);
  //f.add_point_source(Ez, 0.7, 2.5, 0.0, 4.0, vec(0.5,0.5));
  f.add_point_source(Hz, 0.8, 0.6, 0.0, 4.0, vec(0.5,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 ,  0.5))) return 0;
    if (!compare_point(f, f1, vec(0.21 ,  0.5))) 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 exact_pml_rot2x_tm(double eps(const vec &)) {
  double a = 8.0;
  double ttot = 30.0;

  const volume v = voltwo(3.0, 3.0, a);
  the_center = v.center();
  const symmetry S = rotate2(X,v);

  structure s(v, eps, pml(1.0), S);
  structure s1(v, eps, pml(1.0), identity());
  s.set_output_directory(mydirname);
  s1.set_output_directory(mydirname);
  master_printf("Testing X twofold rotational symmetry with PML...\n");

  fields f1(&s1);
  f1.add_point_source(Hx, 0.7, 2.5, 0.0, 4.0, vec(1.3,1.5));
  fields f(&s);
  f.add_point_source(Hx, 0.7, 2.5, 0.0, 4.0, vec(1.3,1.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 ,  1.5))) return 0;
    if (!compare_point(f, f1, vec(1.21 ,  1.5))) return 0;
    if (!compare_point(f, f1, vec(1.46 , 0.33))) return 0;
    if (!compare_point(f, f1, vec(1.2  , 1.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;
}

double polariton_ex(const volume &v, double eps(const vec &)) {
  const double ttot = 10.0;
  master_printf("Testing polariton in %s...\n", dimension_name(v.dim));
  the_center = v.center();
  const symmetry S = mirror(Z,v);
  structure s(v, eps);
  structure sS(v, eps, no_pml(), S);
  s.add_polarizability(one, 0.3, 0.1, 7.63);
  sS.add_polarizability(one, 0.3, 0.1, 7.63);
  fields f(&s);
  f.use_real_fields();
  f.add_point_source(Ex, 0.2, 3.0, 0.0, 2.0, v.center());
  fields fS(&sS);
  fS.use_real_fields();
  fS.add_point_source(Ex, 0.2, 3.0, 0.0, 2.0, v.center());
  f.use_bloch(zero_vec(v.dim));
  fS.use_bloch(zero_vec(v.dim));
  check_unequal_layout(f, fS);
  while (f.time() < ttot) {
    f.step();
    fS.step();
    if (!compare_point(fS, f, v.center())) return 0;
    if (!compare_point(fS, f, zero_vec(v.dim))) return 0;
    if (!compare_point(fS, f, v.center()*0.3)) return 0;
  }
  return 1;
}

double nonlinear_ex(const volume &v, double eps(const vec &)) {
  const double ttot = 10.0;
  master_printf("Testing nonlinear in %s...\n", dimension_name(v.dim));
  the_center = v.center();
  const symmetry S = mirror(Z,v);
  structure s(v, eps);
  structure sS(v, eps, no_pml(), S);
  s.set_chi3(one);
  sS.set_chi3(one);
  fields f(&s);
  f.use_real_fields();
  f.add_point_source(Ex, 0.2, 3.0, 0.0, 2.0, v.center());
  fields fS(&sS);
  fS.use_real_fields();
  fS.add_point_source(Ex, 0.2, 3.0, 0.0, 2.0, v.center());
  f.use_bloch(zero_vec(v.dim));
  fS.use_bloch(zero_vec(v.dim));
  check_unequal_layout(f, fS);
  while (f.time() < ttot) {
    f.step();
    fS.step();
    if (!compare_point(fS, f, v.center())) return 0;
    if (!compare_point(fS, f, zero_vec(v.dim))) return 0;
    if (!compare_point(fS, f, v.center()*0.3)) return 0;
  }
  return 1;
}

double saturated_gain_ez(const volume &v, double eps(const vec &)) {
  const double ttot = 10.0;
  master_printf("Testing saturated gain in %s...\n", dimension_name(v.dim));
  the_center = v.center();
  const symmetry S = -mirror(Z,v);
  structure s(v, eps);
  structure sS(v, eps, no_pml(), S);
  s.add_polarizability(one, 0.3, -0.1, 7.63, 0.5);
  sS.add_polarizability(one, 0.3, -0.1, 7.63, 0.5);
  fields f(&s);
  f.use_real_fields();
  f.add_point_source(Ez, 0.2, 3.0, 0.0, 2.0, v.center());
  fields fS(&sS);
  fS.use_real_fields();
  fS.add_point_source(Ez, 0.2, 3.0, 0.0, 2.0, v.center());
  f.use_bloch(zero_vec(v.dim));
  fS.use_bloch(zero_vec(v.dim));
  check_unequal_layout(f, fS);
  while (f.time() < ttot) {
    f.step();
    fS.step();
    if (!compare_point(fS, f, v.center())) return 0;
    if (!compare_point(fS, f, zero_vec(v.dim))) return 0;
    if (!compare_point(fS, f, v.center()*0.3)) return 0;
  }
  return 1;
}

double saturated_gain_te(const volume &v, double eps(const vec &)) {
  const double ttot = 10.0;
  master_printf("Testing saturated gain in %s...\n", dimension_name(v.dim));
  the_center = v.center();
  const symmetry S = -mirror(X,v);
  structure s(v, eps);
  structure sS(v, eps, no_pml(), S);
  s.add_polarizability(one, 0.3, -0.1, 7.63, 0.5);
  sS.add_polarizability(one, 0.3, -0.1, 7.63, 0.5);
  fields f(&s);
  f.use_real_fields();
  f.add_point_source(Ex, 0.2, 3.0, 0.0, 2.0, v.center());
  fields fS(&sS);
  fS.use_real_fields();
  fS.add_point_source(Ex, 0.2, 3.0, 0.0, 2.0, v.center());
  f.use_bloch(zero_vec(v.dim));
  fS.use_bloch(zero_vec(v.dim));
  check_unequal_layout(f, fS);
  while (f.time() < ttot) {
    f.step();
    fS.step();
    if (!compare_point(fS, f, v.center())) return 0;
    if (!compare_point(fS, f, zero_vec(v.dim))) return 0;
    if (!compare_point(fS, f, v.center()*0.3)) return 0;
  }
  return 1;
}

int main(int argc, char **argv) {
  initialize mpi(argc, argv);
  quiet = true;
  trash_output_directory(mydirname);
  master_printf("Testing with various kinds of symmetry...\n");

  if (0) // FIXME: disable until divergence bug is fixed
  if (!saturated_gain_ez(vol3d(0.5, 1.2, 0.8, 10.0), one))
    abort("error in 3D saturated gain\n");

  if (0) // FIXME: disable until divergence bug is fixed
  if (!saturated_gain_ez(volcyl(0.5, 1.2, 10.0), one))
    abort("error in cylindrical saturated gain\n");

  if (0) // FIXME: disable until divergence bug is fixed
  if (!saturated_gain_te(vol2d(0.6, 1.2, 10.0), one))
    abort("error in 2D TE saturated gain\n");

  if (!test_1d_periodic_mirror(one))
    abort("error in test_1d_periodic_mirror vacuum\n");

  if (!test_cyl_metal_mirror(one))
    abort("error in test_cyl_metal_mirror vacuum\n");

  if (!test_yperiodic_ymirror(one))
    abort("error in test_yperiodic_ymirror vacuum\n");
  if (!test_yperiodic_ymirror(rods_2d))
    abort("error in test_yperiodic_ymirror rods2d\n");

  if (!pml_twomirrors(one))
    abort("error in pml_twomirrors vacuum\n");

  if (!test_origin_shift())
    abort("error in test_origin_shift\n");

  if (!exact_pml_rot2x_tm(one))
    abort("error in exact_pml_rot2x_tm vacuum\n");

  if (!test_metal_xmirror(rods_2d))
    abort("error in test_metal_xmirror rods_2d\n");

  if (!test_metal_xmirror(one))
    abort("error in test_metal_xmirror vacuum\n");
  if (!test_metal_ymirror(one))
    abort("error in test_metal_ymirror vacuum\n");
  if (!test_metal_ymirror(rods_2d))
    abort("error in test_metal_ymirror rods_2d\n");

  if (!test_metal_rot2y(one))
    abort("error in test_metal_rot2y vacuum\n");
  if (!test_metal_rot2y(rods_2d))
    abort("error in test_metal_rot2y rods_2d\n");

  if (!exact_metal_rot2y(one))
    abort("error in exact_metal_rot2y vacuum\n");
  if (!exact_metal_rot2y(rods_2d))
    abort("error in exact_metal_rot2y rods_2d\n");

  if (!exact_metal_rot4z(one))
    abort("error in exact_metal_rot4z vacuum\n");
  if (!exact_metal_rot4z(rods_2d))
    abort("error in exact_metal_rot4z rods_2d\n");


  if (!test_3D_metal_xmirror(one))
    abort("error in test_3D_metal_xmirror vacuum\n");

  if (!test_3D_metal_zmirror(one))
    abort("error in test_3D_metal_zmirror vacuum\n");

  if (!test_3D_metal_odd_zmirror(one))
    abort("error in test_3D_metal_odd_zmirror vacuum\n");

  if (!test_3D_metal_rot4z(one)) {
    all_wait();
    abort("error in test_3D_metal_rot4z vacuum\n");
  }

  if (!test_3D_metal_rot4z_mirror(one))
    abort("error in test_3D_metal_rot4z_mirror vacuum\n");

  if (!test_3D_metal_3mirror(one))
    abort("error in test_3D_metal_3mirror\n");

  /**************************************************************************/
  /* For the following tests, we increase the check tolerance slightly.
     Floating-point errors can cause these tests to have slightly different
     results with and without symmetry. 

     Note also that symmetry is tricky with nonlinearity, since in
     general a nonlinear system does *not* conserve the irreducible
     representation of the symmetry group (i.e. symmetry doesn't work).
     The simulations here are chosen to preserve the symmetry, however. */
  thresh_compare = 1e-10;

  if (!nonlinear_ex(vol1d(1.0, 30.0), one))
    abort("error in 1D nonlinear vacuum\n");
  if (!nonlinear_ex(vol3d(1.0, 1.2, 0.8, 10.0), one))
    abort("error in 3D nonlinear vacuum\n");

  if (!test_cyl_metal_mirror_nonlinear(one))
    abort("error in test_cyl_metal_mirror nonlinear vacuum\n");

  if (!exact_metal_rot4z_nonlinear(one))
    abort("error in exact_metal_rot4z nonlinear vacuum\n");
  if (!exact_metal_rot4z_nonlinear(rods_2d))
    abort("error in exact_metal_rot4z nonlinear rods_2d\n");

  // I'm not sure why the polariton tests require increased tolerances...?
  if (!polariton_ex(vol1d(1.0, 30.0), one))
    abort("error in 1D polariton vacuum\n");
  if (!polariton_ex(vol3d(1.0, 1.2, 0.8, 10.0), one))
    abort("error in 3D polariton vacuum\n");

  return 0;
}