Function test_pml

tests/two_dimensional.cpp:184–227

}

int test_pml(double eps(const vec &), int splitting) {
  double a = 10.0;

  grid_volume gv = voltwo(3.0, 2.0, a);
  structure s1(gv, eps, pml(1.0, X) + pml(1.0, Y, High));
  structure s(gv, eps, pml(1.0, X) + pml(1.0, Y, High), identity(), splitting);

  master_printf("Testing pml while splitting into %d chunks...\n", splitting);
  fields f(&s);
  f.add_point_source(Hz, 0.7, 1.5, 0.0, 4.0, vec(1.5, 0.5), 1.0);
  f.add_point_source(Ez, 0.8, 1.6, 0.0, 4.0, vec(1.299, 0.401), 1.0);
  fields f1(&s1);
  f1.add_point_source(Hz, 0.7, 1.5, 0.0, 4.0, vec(1.5, 0.5), 1.0);
  f1.add_point_source(Ez, 0.8, 1.6, 0.0, 4.0, vec(1.299, 0.401), 1.0);
  const double deltaT = 100.0;
  const double ttot = 3.1 * deltaT;
  double field_energy_check_time = deltaT;

  while (f.time() < f.last_source_time())
    f.step();
  while (f1.time() < f1.last_source_time())
    f1.step();

  double last_energy = f.field_energy();
  while (f.time() < ttot) {
    f.step();
    f1.step();
    if (f.time() >= field_energy_check_time) {
      if (!compare_point(f, f1, vec(0.5, 0.01))) return 0;
      if (!compare_point(f, f1, vec(0.46, 0.33))) return 0;
      if (!compare_point(f, f1, vec(1.0, 1.0))) return 0;
      const double new_energy = f.field_energy();
      if (!compare(new_energy, f1.field_energy(), "   total energy")) return 0;
      if (new_energy > last_energy * 1e-6) {
        master_printf("Energy decaying too slowly: from %g to %g (%g)\n", last_energy, new_energy,
                      new_energy / last_energy);
        return 0;
      }
      else { master_printf("Got newE/oldE of %g\n", new_energy / last_energy); }
      field_energy_check_time += deltaT;
    }
  }
  return 1;
}

int test_pml_tm(double eps(const vec &), int splitting) {
  double a = 10.0;